Changing residual coupling with rho12 • drmTMB

This tutorial introduces bivariate location-scale regression with predictor-dependent residual correlation. The scientific idea is that two responses can change in mean, change in residual variance, and change in residual coupling.

Current fixed-effect syntax uses one formula per distributional parameter:

drmTMB(
  drm_formula(
    mu1 = y1 ~ x1 + x2,
    mu2 = y2 ~ x1,
    sigma1 = ~ x1 + x2,
    sigma2 = ~ x1,
    rho12 = ~ x1 + x2
  ),
  family = c(gaussian(), gaussian()),
  data = dat
)

The central modelling question is whether predictors change residual coupling between two responses, not only their means or variances.

The key parameter is rho12. A location-scale model asks whether predictors change residual SDs. A location-coscale model asks whether predictors also change the residual covariance structure between two responses; in the current bivariate Gaussian model, that structure is the residual correlation rho12. In ecology and evolution, this can ask whether the residual association between traits such as body mass and litter size differs among terrestrial, aquatic, and aerial lifestyles after modelling response-specific means and dispersion.

Equations for the implemented model

For observation i, the implemented fixed-effect bivariate Gaussian model is:

$\begin{bmatrix} y_{1i} \\ y_{2i} \end{bmatrix} \mid \mu_{1i}, \mu_{2i}, \sigma_{1i}, \sigma_{2i}, \rho_{12i} \sim \operatorname{MVN} \left( \begin{bmatrix} \mu_{1i} \\ \mu_{2i} \end{bmatrix}, \Omega_i \right).$

The five distributional predictors are:

$\begin{aligned} \mu_{1i} &= X_{\mu 1,i}\beta_{\mu 1}, \\ \mu_{2i} &= X_{\mu 2,i}\beta_{\mu 2}, \\ \log(\sigma_{1i}) &= X_{\sigma 1,i}\beta_{\sigma 1}, \\ \log(\sigma_{2i}) &= X_{\sigma 2,i}\beta_{\sigma 2}, \\ \eta_{\rho12,i} &= X_{\rho12,i}\beta_{\rho12}, \\ \rho_{12i} &= \tanh(\eta_{\rho12,i}). \end{aligned}$

The residual covariance matrix is:

$\Omega_i = \begin{bmatrix} \sigma_{1i}^2 & \rho_{12i}\sigma_{1i}\sigma_{2i} \\ \rho_{12i}\sigma_{1i}\sigma_{2i} & \sigma_{2i}^2 \end{bmatrix}.$

Component meanings:

Component	Meaning
`mu1`, `mu2`	expected values for responses 1 and 2
`sigma1`, `sigma2`	residual SDs for responses 1 and 2
`eta_rho12`	unconstrained linear predictor for residual correlation
`rho12`	residual response-response correlation within observation `i`
`Omega_i`	residual covariance matrix after means and residual SDs are modelled

For implementation cross-checking, the same model can also be written in a compact code-like form:

mu1_i = X_mu1[i, ] beta_mu1
mu2_i = X_mu2[i, ] beta_mu2
log(sigma1_i) = X_sigma1[i, ] beta_sigma1
log(sigma2_i) = X_sigma2[i, ] beta_sigma2
eta_rho12_i = X_rho12[i, ] beta_rho12
rho12_i = tanh(eta_rho12_i)
Omega_i[1, 1] = sigma1_i^2
Omega_i[2, 2] = sigma2_i^2
Omega_i[1, 2] = Omega_i[2, 1] = rho12_i * sigma1_i * sigma2_i

Implementation note: drmTMB internally evaluates residual correlations as 0.99999999 * tanh(eta_rho12_i). The tiny multiplier keeps covariance matrices strictly positive definite near the boundaries. It is not a biological scaling factor, and it should not change how users interpret the model.

The matching drmTMB syntax is:

drmTMB(
  bf(
    mu1 = y1 ~ x1 + x2,
    mu2 = y2 ~ x1,
    sigma1 = ~ x1 + x2,
    sigma2 = ~ x1,
    rho12 = ~ x1 + x2
  ),
  family = c(gaussian(), gaussian()),
  data = dat
)

The mapping is:

Symbol	Meaning	R syntax source
`y1_i`, `y2_i`	observed responses	left-hand sides of `mu1` and `mu2` formulas
`X_mu1`, `X_mu2`	response-specific mean design matrices	right-hand sides of `mu1` and `mu2`
`mu1_i`, `mu2_i`	fitted response-specific means	columns of `fitted(fit)`
`X_sigma1`, `X_sigma2`	response-specific residual scale design matrices	right-hand sides of `sigma1` and `sigma2`
`X_rho12`	residual correlation design matrix	right-hand side of `rho12`
`sigma1_i`, `sigma2_i`	residual SDs for responses 1 and 2	`predict(fit, dpar = "sigma1")`, `predict(fit, dpar = "sigma2")`
`rho12_i`	residual response-response correlation	`rho12(fit)`

The fixed effects may be the same, overlapping, or different across parameters. For example, x1 can appear in all five formulae, while x2 might appear only in mu1, sigma1, and rho12. Each parameter still gets its own coefficient vector.

Read those coefficient vectors as separate biological statements:

Formula	Example coefficient	Interpretation
`mu1 = y1 ~ x`	`mu1:x`	predictor effect on the expected value of response 1
`mu2 = y2 ~ x`	`mu2:x`	predictor effect on the expected value of response 2
`sigma1 = ~ x`	`sigma1:x`	multiplicative change in residual SD for response 1
`sigma2 = ~ x`	`sigma2:x`	multiplicative change in residual SD for response 2
`rho12 = ~ x`	`rho12:x`	change in residual response-response coupling after both means and residual SDs are modelled

The rho12 slope is not a slope for either response mean. It is a slope on the correlation-link scale. Use rho12(fit, newdata = grid) or a rho12 curve when the reader needs the fitted residual correlation on the response scale.

If both responses have the same location predictors, mvbind() is a shorthand for the two location formulas:

drmTMB(
  bf(
    mvbind(y1, y2) ~ x1 + x2,
    sigma1 = ~ x1 + x2,
    sigma2 = ~ x1,
    rho12 = ~ x1 + x2
  ),
  family = c(gaussian(), gaussian()),
  data = dat
)

This expands internally to mu1 = y1 ~ x1 + x2 and mu2 = y2 ~ x1 + x2. Use explicit mu1 and mu2 formulas whenever the responses need different location predictors.

Worked example: behaviour coupling and disturbance

Suppose an ecologist measures two behaviours, activity and boldness, and wants to know whether an environmental gradient changes their residual coupling after accounting for response-specific means and residual SDs. The code below uses simulated data, but the model structure matches that scientific question.

For this example, the fitted biological model is:

$\begin{aligned} \begin{bmatrix} \text{activity}_i\\ \text{boldness}_i \end{bmatrix} &\sim \operatorname{MVN} \left( \begin{bmatrix} \mu_{\text{activity},i}\\ \mu_{\text{boldness},i} \end{bmatrix}, \Omega_i \right),\\ \mu_{\text{activity},i} &= \beta_{a0} + \beta_{a1}\text{food}_i + \beta_{a2}\text{temperature}_i,\\ \mu_{\text{boldness},i} &= \beta_{b0} + \beta_{b1}\text{food}_i,\\ \log(\sigma_{\text{activity},i}) &= \gamma_{a0} + \gamma_{a1}\text{food}_i + \gamma_{a2}\text{temperature}_i,\\ \log(\sigma_{\text{boldness},i}) &= \gamma_{b0} + \gamma_{b1}\text{food}_i,\\ \eta_{\rho12,i} &= \delta_0 + \delta_1\text{disturbance}_i,\\ \rho_{12,i} &= \tanh(\eta_{\rho12,i}). \end{aligned}$

Here rho12 is not the ordinary correlation between raw activity and raw boldness. It is the residual correlation after the model has already adjusted for the mean and residual SD predictors in each response.

set.seed(1)
n <- 180
dat <- data.frame(
  food = rnorm(n),
  temperature = rnorm(n),
  disturbance = rnorm(n)
)
mu1 <- 0.2 + 0.5 * dat$food + 0.2 * dat$temperature
mu2 <- -0.1 + 0.4 * dat$food
sigma1 <- exp(-0.2 + 0.2 * dat$food - 0.1 * dat$temperature)
sigma2 <- exp(0.1 - 0.2 * dat$food)
eta_rho12 <- -0.1 + 0.4 * dat$disturbance
rho12 <- tanh(eta_rho12)
e1 <- rnorm(n)
e2 <- rho12 * e1 + sqrt(1 - rho12^2) * rnorm(n)
dat$activity <- mu1 + sigma1 * e1
dat$boldness <- mu2 + sigma2 * e2

The model lets food and temperature explain mean activity, lets food explain mean boldness, lets food and temperature explain residual SD in activity, lets food explain residual SD in boldness, and lets disturbance explain residual correlation:

fit_biv <- drmTMB(
  drm_formula(
    mu1 = activity ~ food + temperature,
    mu2 = boldness ~ food,
    sigma1 = ~ food + temperature,
    sigma2 = ~ food,
    rho12 = ~ disturbance
  ),
  family = c(gaussian(), gaussian()),
  data = dat
)

Check the fit and inspect the coefficient table:

check_drm(fit_biv)
#> <drm_check: 11 checks>
#> ok: 11; notes: 0; warnings: 0; errors: 0
#>                      check status
#>      optimizer_convergence     ok
#>           optimizer_budget     ok
#>           finite_objective     ok
#>             fixed_gradient     ok
#>            sdreport_status     ok
#>  hessian_positive_definite     ok
#>     standard_errors_finite     ok
#>               dropped_rows     ok
#>             positive_scale     ok
#>             rho12_boundary     ok
#>   fixed_effect_design_size     ok
#>                                                                                   value
#>                                                                                       0
#>                                                 iterations=19; function=32; gradient=19
#>                                                                                   495.7
#>                                                    max=0.0001227; component=beta_mu1[2]
#>                                                                                      ok
#>                                                                                    TRUE
#>                                                                 range=[0.05053,0.08091]
#>                                                                     nobs=180; dropped=0
#>                                                                              min=0.5408
#>                                                                                  0.9633
#>  total_mb=0.07536; max_cols=3; largest=mu1; largest_class=matrix; largest_density=1.000
#>                                                                              message
#>                                                        nlminb convergence code is 0.
#>  Optimizer evaluation counts recorded; no eval.max or iter.max control was supplied.
#>                                             Objective and log-likelihood are finite.
#>       Maximum absolute fixed gradient is <= 0.001; largest component is beta_mu1[2].
#>                                              TMB::sdreport() completed successfully.
#>                                        sdreport reports a positive-definite Hessian.
#>                                         All fixed-effect standard errors are finite.
#>                   No rows were dropped by model-frame or known-covariance filtering.
#>                                     All fitted scale values are finite and positive.
#>                        All fitted residual correlations have absolute value <= 0.98.
#>                          Dense fixed-effect design matrices are modest for this fit.
summary(fit_biv)
#> <summary.drmTMB>
#>                        estimate  std_error
#> mu1:(Intercept)     0.130704430 0.05416607
#> mu1:food            0.565332984 0.06286362
#> mu1:temperature     0.217923200 0.05437147
#> mu2:(Intercept)    -0.081405830 0.07948982
#> mu2:food            0.330438602 0.08090952
#> sigma1:(Intercept) -0.155133496 0.05052620
#> sigma1:food         0.208914220 0.05454897
#> sigma1:temperature  0.003849881 0.05566019
#> sigma2:(Intercept)  0.209287360 0.05234177
#> sigma2:food        -0.192523348 0.05744740
#> rho12:(Intercept)  -0.004443165 0.07014707
#> rho12:disturbance   0.523480249 0.06355269
#> Distributional, random-effect, scale, and correlation parameters:
#>                          component   dpar         term    estimate    minimum
#> fitted:sigma1 distributional-scale sigma1 fitted range  0.88343668  0.5408469
#> fitted:sigma2 distributional-scale sigma2 fitted range  1.23862903  0.7764057
#> fitted:rho12  residual-correlation  rho12 fitted range -0.02545558 -0.9184608
#>                 maximum    scale
#> fitted:sigma1 1.4181053 response
#> fitted:sigma2 1.8882846 response
#> fitted:rho12  0.9633258 response
#> logLik: -495.7
#> convergence: 0

How to read this output:

Rows beginning with mu1: describe expected activity, and rows beginning with mu2: describe expected boldness.
Rows beginning with sigma1: and sigma2: describe log residual SDs for activity and boldness.
Rows beginning with rho12: describe the residual activity-boldness correlation on the linear-predictor scale. Use rho12(fit_biv) or rho12(fit_biv, newdata = ...) to read that effect as a correlation.

The rho12 coefficients are on the linear-predictor scale:

coef(fit_biv, "rho12")
#>  (Intercept)  disturbance 
#> -0.004443165  0.523480249
head(rho12(fit_biv))
#> [1] -0.8767348  0.5937324 -0.3249194 -0.2255810 -0.0928101  0.3059228

corpairs() gives the same residual correlation in the long correlation-pair format. The same table shape now also holds fitted ordinary group-level and phylogenetic rows; richer spatial and study-level correlation rows remain planned:

corpairs(fit_biv)[
  c("level", "from_response", "to_response", "class", "parameter",
    "estimate", "min", "max", "modelled")
]
#>      level from_response to_response    class parameter    estimate        min
#> 1 residual      activity    boldness residual     rho12 -0.02545558 -0.9184608
#>         max modelled
#> 1 0.9633258     TRUE

The rho12 coefficients are fitted on the linear-predictor scale eta_rho12; the response transform is rho12 = tanh(eta_rho12) with the tiny internal guard described above. Use rho12(fit) to read them as residual correlations, and use corpairs(fit) to place residual rho12 in a table that is ready to hold other correlation levels later. In this example, a positive coefficient for disturbance means that observations with larger disturbance have stronger positive residual coupling between activity and boldness, after the model has already accounted for the predictors in mu1, mu2, sigma1, and sigma2.

newdat <- data.frame(
  food = 0,
  temperature = 0,
  disturbance = c(-1, 0, 1)
)
rho_table <- data.frame(
  disturbance = newdat$disturbance,
  sigma_activity = predict(fit_biv, newdata = newdat, dpar = "sigma1"),
  sigma_boldness = predict(fit_biv, newdata = newdat, dpar = "sigma2"),
  rho12 = rho12(fit_biv, newdata = newdat)
)
rho_table$residual_covariance <- with(
  rho_table,
  rho12 * sigma_activity * sigma_boldness
)
rho_table$residual_variance_activity <- rho_table$sigma_activity^2
rho_table$residual_variance_boldness <- rho_table$sigma_boldness^2
round(rho_table, 3)
#>   disturbance sigma_activity sigma_boldness  rho12 residual_covariance
#> 1          -1          0.856          1.233 -0.484              -0.511
#> 2           0          0.856          1.233 -0.004              -0.005
#> 3           1          0.856          1.233  0.477               0.503
#>   residual_variance_activity residual_variance_boldness
#> 1                      0.733                       1.52
#> 2                      0.733                       1.52
#> 3                      0.733                       1.52

This table gives three quantities a reader can report directly. sigma_activity and sigma_boldness are the fitted residual SDs for each response at the chosen covariate values. rho12 is the fitted residual correlation between the two responses. residual_covariance is the off-diagonal element of \Omega_i, computed as rho12 * sigma_activity * sigma_boldness. The two residual variance columns are the diagonal elements of \Omega_i.

What should I report?

For a location-coscale model, a useful result sentence usually needs all three distributional pieces: means, residual SDs, and residual coupling. A minimal table for the coscale part is:

report_grid <- data.frame(
  food = 0,
  temperature = 0,
  disturbance = c(-1.5, 0, 1.5)
)
report_table <- data.frame(
  disturbance = report_grid$disturbance,
  eta_rho12 = predict(
    fit_biv,
    newdata = report_grid,
    dpar = "rho12",
    type = "link"
  ),
  rho12 = rho12(fit_biv, newdata = report_grid),
  sigma_activity = predict(fit_biv, newdata = report_grid, dpar = "sigma1"),
  sigma_boldness = predict(fit_biv, newdata = report_grid, dpar = "sigma2")
)
report_table$residual_covariance <- with(
  report_table,
  rho12 * sigma_activity * sigma_boldness
)
report_table$residual_variance_activity <- report_table$sigma_activity^2
report_table$residual_variance_boldness <- report_table$sigma_boldness^2
round(report_table, 3)
#>   disturbance eta_rho12  rho12 sigma_activity sigma_boldness
#> 1        -1.5    -0.790 -0.658          0.856          1.233
#> 2         0.0    -0.004 -0.004          0.856          1.233
#> 3         1.5     0.781  0.653          0.856          1.233
#>   residual_covariance residual_variance_activity residual_variance_boldness
#> 1              -0.695                      0.733                       1.52
#> 2              -0.005                      0.733                       1.52
#> 3               0.689                      0.733                       1.52

The raw correlation between activity and boldness is not the same estimand as rho12. The raw correlation mixes mean structure, residual SDs, and residual coupling. The fitted rho12 is the residual correlation after the model has accounted for the formulae in mu1, mu2, sigma1, and sigma2:

round(data.frame(
  raw_activity_boldness_correlation = stats::cor(dat$activity, dat$boldness),
  mean_fitted_residual_rho12 = mean(rho12(fit_biv))
), 3)
#>   raw_activity_boldness_correlation mean_fitted_residual_rho12
#> 1                             0.021                     -0.025

When reporting this model, avoid saying only “activity and boldness were correlated”. A better sentence is: “After modelling response-specific means and residual SDs, the fitted residual activity-boldness correlation increased with disturbance.” That sentence tells the reader which correlation has been modelled.

The same response-scale result can be visualised as a rho12 curve:

rho_grid <- data.frame(
  food = 0,
  temperature = 0,
  disturbance = seq(-2, 2, length.out = 80)
)
rho_grid$rho12 <- rho12(fit_biv, newdata = rho_grid)

plot(
  rho12 ~ disturbance,
  data = rho_grid,
  type = "l",
  lwd = 2,
  ylim = c(-1, 1),
  xlab = "Disturbance",
  ylab = "Fitted residual correlation (rho12)"
)
abline(h = 0, lty = 2, col = "grey60")

In this example, a positive disturbance coefficient for rho12 means that activity and boldness become more tightly positively coupled at higher disturbance, after their means and residual SDs have been modelled. It does not say that disturbance changes among-individual personality or plasticity correlations; those require group-level covariance blocks and should appear as separate rows in corpairs(fit).

A concise reporting sentence would be: after accounting for food, temperature, and response-specific residual SDs, the fitted residual correlation between activity and boldness increased along the disturbance gradient. In biological terms, individuals that were more active than expected were also more bold than expected under higher disturbance.

Where residual rho12 stops

The fitted model above has no repeated-individual structure. It estimates residual coupling: whether two measurements from the same observation are closer together or farther apart than expected after modelling their means and residual SDs. That is useful, but it is not the same question as whether some individuals have consistently high activity and consistently high boldness.

For ecological and evolutionary examples, residual rho12 can represent changing residual coupling between traits such as activity and boldness, body size and fecundity, or leaf area and seed number. A group-level covariance block asks a different question: do groups, individuals, species, or sites share latent deviations across fitted components?

Worked example: individual differences

Suppose the same ecologist repeatedly measures each individual and now asks an individual-difference question. After accounting for food and disturbance, do individuals with higher average activity also tend to have higher average boldness?

The first fitted group-level version adds matching labelled random intercepts in the mean structure. The shared p label creates one covariance block for the mu1 and mu2 random intercepts:

drmTMB(
  formula = bf(
    mu1 = activity ~ food + disturbance + (1 | p | ID),
    mu2 = boldness ~ food              + (1 | p | ID),
    sigma1 = ~1,
    sigma2 = ~1,
    rho12 = ~1
  ),
  family = biv_gaussian(),
  data = dat_group
)

Here sigma1 and sigma2 are residual scale parameters. The fitted random-intercept SDs are group-level standard deviations in sdpars$mu, not residual sigma.

Symbolically, the implemented random-intercept mean part is:

$\begin{aligned} \mu_{1ij} &= X_{\mu 1}[ij, ]\beta_{\mu 1} + b_{0,1j},\\ \mu_{2ij} &= X_{\mu 2}[ij, ]\beta_{\mu 2} + b_{0,2j},\\ \begin{bmatrix} b_{0,1j}\\ b_{0,2j} \end{bmatrix} &\sim \operatorname{MVN}(\mathbf{0}, \Sigma_{\mu,ID}). \end{aligned}$

The covariance matrix Sigma_mu_ID contains two group-level standard deviations and one group-level correlation. That correlation answers whether individuals with higher average activity also tend to have higher average boldness. It is not rho12_i. The parameter rho12_i remains the observation-level residual correlation in Omega_i.

The small simulation below creates repeated observations from that model:

set.seed(98)
n_ID <- 55
n_each <- 7
ID <- factor(rep(seq_len(n_ID), each = n_each))
id_index <- as.integer(ID)
n <- length(ID)

food <- rnorm(n)
disturbance <- rnorm(n)

sigma_activity <- 0.35
sigma_boldness <- 0.45
sd_activity_ID <- 0.65
sd_boldness_ID <- 0.65
rho_ID <- 0.55
rho_residual <- 0.15

u_activity <- rnorm(n_ID)
u_boldness <- rho_ID * u_activity + sqrt(1 - rho_ID^2) * rnorm(n_ID)

e_activity <- rnorm(n)
e_boldness <- rho_residual * e_activity +
  sqrt(1 - rho_residual^2) * rnorm(n)

dat_group <- data.frame(ID = ID, food = food, disturbance = disturbance)
dat_group$activity <- 0.25 + 0.45 * food + 0.15 * disturbance +
  sd_activity_ID * u_activity[id_index] + sigma_activity * e_activity
dat_group$boldness <- -0.15 + 0.35 * food +
  sd_boldness_ID * u_boldness[id_index] + sigma_boldness * e_boldness

Fit the grouped model with matching (1 | p | ID) terms in the two location formulas:

fit_group <- drmTMB(
  bf(
    mu1 = activity ~ food + disturbance + (1 | p | ID),
    mu2 = boldness ~ food + (1 | p | ID),
    sigma1 = ~1,
    sigma2 = ~1,
    rho12 = ~1
  ),
  family = biv_gaussian(),
  data = dat_group,
  control = list(eval.max = 500, iter.max = 500)
)

In a real analysis, inspect every row of check_drm(fit_group). The filtered rows below show the convergence and covariance checks most directly tied to this individual-difference question:

group_checks <- check_drm(fit_group)
group_checks[
  group_checks$check %in% c(
    "optimizer_convergence",
    "random_effect_sd_boundary",
    "rho12_boundary",
    "biv_mu_random_effect_covariance"
  ),
  c("check", "status", "value", "message")
]
#> <drm_check: 4 checks>
#> ok: 4; notes: 0; warnings: 0; errors: 0
#>                            check status
#>            optimizer_convergence     ok
#>        random_effect_sd_boundary     ok
#>                   rho12_boundary     ok
#>  biv_mu_random_effect_covariance     ok
#>                                                               value
#>                                                                   0
#>            min=0.6263; boundary=0.0001000; term=mu.mu1:(1 | p | ID)
#>                                                              0.1743
#>  n_groups=55; min_group_n=7; singleton_groups=0; min_sd_ratio=1.414
#>                                                                                                                       message
#>                                                                                                 nlminb convergence code is 0.
#>  All fitted random-effect standard deviations are finite, positive, and above the requested lower-boundary warning threshold.
#>                                                                 All fitted residual correlations have absolute value <= 0.98.
#>             Bivariate group-level covariance has replicated groups and non-negligible fitted SDs relative to residual scales.

Reading group-level covariance

corpairs() keeps residual and group-level correlations in the same long table without giving them the same meaning:

pair_table <- corpairs(fit_group)
pair_table[
  ,
  c(
    "level", "group", "block", "from_dpar", "to_dpar", "class",
    "parameter", "estimate", "min", "max", "modelled"
  )
]
#>      level group block from_dpar  to_dpar     class
#> 1 residual  <NA>  <NA>  residual residual  residual
#> 2    group    ID     p       mu1      mu2 mean-mean
#>                                       parameter  estimate       min       max
#> 1                                         rho12 0.1743289 0.1743289 0.1743289
#> 2 cor(mu1:(Intercept),mu2:(Intercept) | p | ID) 0.4752257 0.4752257 0.4752257
#>   modelled
#> 1    FALSE
#> 2    FALSE

The optional plotting helper consumes that explicit table. Faceting by level keeps the residual rho12 row visually separate from the group-level random-intercept correlation:

if (requireNamespace("ggplot2", quietly = TRUE)) {
  plot_corpairs(pair_table, facet = "level")
}

The residual row is rho12: within-observation activity-boldness coupling. The group row is the mu1/mu2 random-intercept correlation: among-individual covariation in average activity and average boldness.

summary(fit_group)$covariance gives the same group-level row with the component SDs and covariance on the interpretation scale:

group_covariance <- summary(fit_group)$covariance
report_group_covariance <- group_covariance[
  ,
  c(
    "level", "group", "block", "from_response", "to_response", "class",
    "correlation", "from_sd", "to_sd", "covariance", "from_scale",
    "to_scale", "covariance_conf.status"
  )
]
numeric_columns <- vapply(report_group_covariance, is.numeric, logical(1))
report_group_covariance[numeric_columns] <- lapply(
  report_group_covariance[numeric_columns],
  round,
  3
)
report_group_covariance
#>   level group block from_response to_response     class correlation from_sd
#> 1 group    ID     p      activity    boldness mean-mean       0.475   0.626
#>   to_sd covariance from_scale to_scale covariance_conf.status
#> 1 0.643      0.191   identity identity          not_requested

For this mean-mean block, from_sd is the among-individual SD in average activity and to_sd is the among-individual SD in average boldness. The correlation column is the individual-level activity-boldness correlation, and covariance is correlation * from_sd * to_sd. The identity scale columns tell you these quantities are on the response-location scale, not on a log residual-SD scale.

For this fitted block, check_drm() reports whether group levels are replicated and whether either group-level SD is tiny relative to the matching residual scale. A note there is not a proof that the model is wrong; it tells you to inspect whether the group-level SDs and correlation are supported by the data before treating them as biological conclusions.

Use profile_targets(fit_group) before requesting profile intervals for this block:

group_targets <- profile_targets(fit_group)
group_target_names <- c(
  "sd:mu:mu1:(1 | p | ID)",
  "sd:mu:mu2:(1 | p | ID)",
  "cor:mu:cor(mu1:(Intercept),mu2:(Intercept) | p | ID)"
)
group_targets[
  match(group_target_names, group_targets$parm),
  c("parm", "profile_ready", "profile_note")
]
#>                                                    parm profile_ready
#> 12                               sd:mu:mu1:(1 | p | ID)          TRUE
#> 13                               sd:mu:mu2:(1 | p | ID)          TRUE
#> 14 cor:mu:cor(mu1:(Intercept),mu2:(Intercept) | p | ID)          TRUE
#>    profile_note
#> 12        ready
#> 13        ready
#> 14        ready

The group-level SD targets are named like sd:mu:mu1:(1 | p | id) and sd:mu:mu2:(1 | p | id), while the group-level correlation target is named like cor:mu:cor(mu1:(Intercept),mu2:(Intercept) | p | id). These are separate from the residual-correlation target rho12. Returned profile rows include profile.boundary and profile.message; a non-"ok" message is a prompt to inspect check_drm(fit) and the profile path before treating the interval as a stable biological estimate.

Boundaries for this covariance layer

The middle p in (1 | p | ID) follows grouped covariance-block syntax for correlated group-level effects. This syntax is distinct from residual rho12. Use rho12(fit) for within-observation residual coupling and corpairs(fit, level = "group") for fitted ordinary group-level covariance rows.

This example teaches only the mu1/mu2 random-intercept covariance row. Other implemented ordinary bivariate random-intercept rows include sigma1/sigma2, one same-response mu/sigma pair, and the all-four mu1/mu2/sigma1/sigma2 block. Read the all-four block as six rows: one mu1-mu2 row, four mean-scale rows (mu1-sigma1, mu1-sigma2, mu2-sigma1, and mu2-sigma2), and one sigma1-sigma2 row.

Bivariate random slopes, slope1-slope2 plasticity-syndrome correlations, random effects in rho12, bivariate meta_known_V() plus random effects, mixed-response families, and ordinary spatial group-level covariance remain future work unless a later status map says otherwise. For phylogenetic covariance, use the structural-dependence tutorial and keep those rows separate from both ordinary group-level covariance and residual rho12.