Bivariate Copula Regression with gamboostLSS

Introduction

gamboostLSS fits one additive predictor per distribution parameter. With CopulaFamilies this is extended to a bivariate response: each margin keeps its own parameters, and an additional parameter theta describes how the two responses depend on each other. Every one of these parameters gets its own predictor and its own variable selection.

By Sklar’s theorem the joint density of two continuous responses factorises into a copula density evaluated at the marginal CDFs and the two marginal densities,

\[ f(y_1, y_2) = c(F_1(y_1), F_2(y_2);\theta)f_1(y_1)f_2(y_2) \]

so the joint log-likelihood is a sum of a dependence term and two marginal terms. This is what CopulaFamilies() boosts.

A first model

The response is a two-column matrix. Everything else works exactly as in a normal gamboostLSS() call.

set.seed(1)
n <- 150
x <- rnorm(n)
eps <- matrix(rnorm(2 * n), n, 2) %*% chol(matrix(c(1, 0.3, 0.3, 1), 2, 2))
y <- cbind(x + eps[, 1], 0.5 * x + eps[, 2])

cf <- CopulaFamilies(GaussianLSS(), GaussianLSS(), copula = "clayton")
fit <- gamboostLSS(y ~ bbs(x), families = cf, data = data.frame(x = x),
                   control = boost_control(mstop = 50, nu = 0.1))

names(fit)
#> [1] "mu1"    "sigma1" "mu2"    "sigma2" "theta"

Five sub-models: two parameters per margin and the copula parameter. The suffixes 1 and 2 indicate which margin a parameter belongs to.

predict() with type = "response" returns every parameter on its natural scale, including theta on the copula’s own scale.

pr <- predict(fit, newdata = data.frame(x = c(-1, 0, 1)), type = "response")
round(cbind(mu1 = pr$mu1, sigma1 = pr$sigma1, theta = pr$theta), 3)
#>        [,1]  [,2]  [,3]
#> [1,] -0.952 0.850 0.577
#> [2,]  0.049 1.166 0.578
#> [3,]  0.789 1.049 0.773

Choosing a copula

Four copulas are available. They differ in how they distribute dependence across the range of the data, in particular the tails.

data.frame(
  copula = c("gaussian", "clayton", "gumbel", "frank"),
  range = c("(-1,1)", "(0, Inf)", "[1, Inf)", "R \\ {0}"),
  tail = c("none", "lower", "upper", "none")
)
#>     copula    range  tail
#> 1 gaussian   (-1,1)  none
#> 2  clayton (0, Inf) lower
#> 3   gumbel [1, Inf) upper
#> 4    frank R \\ {0}  none

Clayton captures joint low values, Gumbel joint high values. Gaussian and Frank are symmetric. Note that Clayton and Gumbel cannot represent negative dependence at all.

Because theta lives on a different scale in each family, the fitted values are not comparable across copulas. copula_tau() converts them to Kendall’s tau, which is

taus <- sapply(c("clayton", "gumbel", "frank"), function(cop){
  cf <- CopulaFamilies(GaussianLSS(), GaussianLSS(), copula = cop)
  fit <- gamboostLSS(y ~ bbs(x), families = cf, data = data.frame(x = x),
                     control = boost_control(mstop = 50, nu = 0.1))
  c(theta = mean(predict(fit, type = "response")$theta),
    tau = mean(copula_tau(fit)))
})
round(taus, 3)
#>       clayton gumbel frank
#> theta   0.641  1.372 2.740
#> tau     0.242  0.271 0.284

The theta row is not comparable across columns, while the tau row is.

Discrete and binary margins

The marginal combination is detected automatically, and the likelihood is built accordingly - from copula CDF differences for discrete margins and from the h-function (conditional copula CDF) for mixed cases.

set.seed(42)
y_mix <- cbind(rbinom(n, 1, plogis(x)), rnorm(n))
cf_mix <- CopulaFamilies(BernoulliLSS(), GaussianLSS(), copula = "frank")
fit_mix <- gamboostLSS(y_mix ~ bbs(x), families = cf_mix, 
                       data = data.frame(x = x),
                       control = boost_control(mstop = 30, nu = 0.1))
names(fit_mix)
#> [1] "mu1"    "mu2"    "sigma2" "theta"

The binary margin contributes a single parameter, so this model has four sub-models rather than five.

Selecting mstop

mstop is the main tuning parameter. make.grid() builds the full factorial grid over all parameters, which is infeasible - five parameters at ten candidate values each is \(10^5 = 100000\) combinations. make_copula_grid() instead searches a single shared stopping iteration.

grid <- make_copula_grid(fit, from = 10, by = 10)
dim(grid)
#> [1] 5 5
cvr <- cvrisk(fit, folds = cv(model.weights(fit), type = "subsampling", B = 3),
              grid = grid, papply = lapply, trace = FALSE)
mstop(cvr)
#>    mu1 sigma1    mu2 sigma2  theta 
#>     20     20     20     20     20

The selected values is applied with:

invisible(fit[mstop(cvr)])
mstop(fit)
#>    mu1 sigma1    mu2 sigma2  theta 
#>     20     20     20     20     20

Note that [ modifies the model in place.

A worked example

The india dataset records childhood malnutrition in India (see gamboostLSS_Tutorial for more information). We model two child health outcomes jointly: stunting (height-for-age z-score) and cbmi (BMI of the child), as smooth functions of the child’s age cage.

data(india)
y_india <- cbind(india$stunting, india$cbmi)

cf_india <- CopulaFamilies(GaussianLSS(), GaussianLSS(), copula = "gaussian")
fit_india <- gamboostLSS(y_india ~ bbs(cage), families = cf_india, data = india,
                         control = boost_control(mstop = 100, nu = 0.1))

fv <- lapply(fit_india, fitted)
round(c(mu1 = mean(fv$mu1), mu2 = mean(fv$mu2)), 3)
#>    mu1    mu2 
#> -1.749 15.526
round(c(observed1 = mean(india$stunting), observed2 = mean(india$cbmi)), 3)
#> observed1 observed2 
#>    -1.754    15.525

The fitted means recover the observed ones. The interesting part is the dependence, which is allowed to vary with age:

tau_india <- copula_tau(fit_india)
plot(india$cage, tau_india, pch = 16, cex = 0.4, xlab = "child age (months)",
     ylab = expression(hat(tau)))
abline(h = cor(india$stunting, india$cbmi, method = "kendall"), col = "red", 
       lty = 2)

The dashed line is the unconditional Kendall’s tau of the two responses. The fitted curve shows how the association changes with age - something a single correlation coefficient cannot express.

Practical Notes

Starting values for theta The offset for theta is derived from the empirical Kendall’s tau of the responses, inverted through the copula’s tau-to-theta relationship. When that value lies outside what a copula can represent (e.g. negative dependence with Clayton) a static default is used instead.

Numerical stability Very strong dependence combined with a large step size can destabilise the fit. If a fit aborts with “Infinite residuals” or “could not fit base-learner”, reduce nu! The same data usually fits at nu = 0.05 or 0.01.

Session info

sessionInfo()
#> R version 4.6.1 (2026-06-24)
#> Platform: x86_64-pc-linux-gnu
#> Running under: Ubuntu 26.04.1 LTS
#> 
#> Matrix products: default
#> BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3 
#> LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.32.so;  LAPACK version 3.12.0
#> 
#> locale:
#>  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
#>  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
#>  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
#>  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
#>  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
#> [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
#> 
#> time zone: Etc/UTC
#> tzcode source: system (glibc)
#> 
#> attached base packages:
#> [1] parallel  stats     graphics  grDevices utils     datasets  methods  
#> [8] base     
#> 
#> other attached packages:
#> [1] gamboostLSS_2.2-0  mboost_2.9-14      stabs_0.7-1        gamlss.dist_6.1-11
#> [5] rmarkdown_2.32    
#> 
#> loaded via a namespace (and not attached):
#>  [1] Matrix_1.7-6      strucchange_1.6-0 jsonlite_2.0.0    compiler_4.6.1   
#>  [5] rpart_4.1.27      jquerylib_0.1.4   splines_4.6.1     yaml_2.3.12      
#>  [9] fastmap_1.2.0     lattice_0.23-1    R6_2.6.1          Formula_1.2-6    
#> [13] knitr_1.52        MASS_7.3-66       inum_1.0-5        maketools_1.3.2  
#> [17] bslib_0.12.0      rlang_1.3.0       partykit_1.3-0    cachem_1.1.0     
#> [21] xfun_0.62         nnls_1.6          quadprog_1.5-8    sass_0.4.10      
#> [25] sys_3.4.3         cli_3.6.6         digest_0.6.39     grid_4.6.1       
#> [29] mvtnorm_1.4-2     sandwich_3.1-3    lifecycle_1.0.5   evaluate_1.0.5   
#> [33] libcoin_1.0-13    buildtools_1.0.0  zoo_1.9-1         survival_3.8-12  
#> [37] tools_4.6.1       htmltools_0.5.9