Bivariate nodes: joint estimation, residual correlation, and moderation (drm

In standard piecewise structural equation models (Shipley 2000, 2009; Lefcheck 2016), each response variable is modelled in its own univariate node. However, in behavioural ecology, evolutionary genetics, and physiology, researchers frequently encounter pairs of traits—such as activity and boldness in animal personality (Brooks et al. 2017) or growth and reproduction in life-history tradeoffs—that remain coupled even after conditioning on external covariates.

This coupling can reflect either within-observation residual covariance ( $\varepsilon_1 \leftrightarrow \varepsilon_2$ ) or between-individual random covariance ( $u_{1} \leftrightarrow u_{2}$ ).

drmSEM provides the bivariate node constructor drm_pair() to model these joint response distributions. This vignette explains:

The crucial distinction between directed paths, residual correlations ( $\rho_{12}$ ), and higher-level correlations (corpairs).
Declaring joint bivariate nodes with drm_pair().
Moderation paths on the $\tanh$ link ( $x \to \rho_{12}$ , declared as rho12 = ~ x).
Basis-set suppression: why estimating $\rho_{12}$ removes the independence claim from d-separation.
Post-estimation inspection: rho12(), corpairs(), paths(), and plotting.

The three ways traits couple

When two traits $Y_1$ and $Y_2$ covary in an observational dataset under an environmental predictor $X$ , there are three fundamentally different mechanisms:

Coupling class	Causal interpretation	Model representation
1. Directed path ( $Y_1 \to Y_2$ )	Trait 1 causes Trait 2 (e.g. Activity directly induces Boldness).	Univariate piecewise nodes: `boldness ~ activity + x`.
2. Residual correlation ( $\rho_{12}$ )	Unmodelled common shocks or physiology couple $Y_1$ and $Y_2$ at the observation level ( $\varepsilon_{1} \leftrightarrow \varepsilon_{2}$ ).	Joint bivariate node via `drm_pair(..., rho12 = ...)`.
3. Between-group correlation (`corpair`)	Individual syndromes or genetic linkages couple $Y_1$ and $Y_2$ at the group/individual level ( $u_{\text{id}, 1} \leftrightarrow u_{\text{id}, 2}$ ).	Shared random intercepts: `(1 \\| id)` in both formulas.

Conflating these three processes is a common pitfall. In drmSEM, they are modelled and reported strictly separately.

1. Declaring a bivariate node: `drm_pair()`

drm_pair() specifies a joint two-response node. It takes two response formulas, an optional model for the residual correlation parameter $\rho_{12}$ , and optional grouping levels:

# Declare a bivariate node for activity and boldness
pair_spec <- drm_pair(
  activity ~ temp + (1 | id),
  boldness ~ temp + (1 | id),
  rho12 = ~ temp    # environmental moderation of coupling strength
)
print(pair_spec)
#> 
#> ── <drm_pair> bivariate node "activity" & "boldness"
#> activity [gaussian]: `activity ~ temp + (1 | id)`
#> boldness [gaussian]: `boldness ~ temp + (1 | id)`
#> residual correlation: rho12 ~ temp [directed path into rho12]
#> higher-level correlation: corpair at 1 level ("id")
#> declaration only; pass this pair to drm_sem() for a joint bivariate fit

Accessing declarations before fitting

Before model fitting, rho12() and corpairs() report the declared structure. Estimates are honestly reported as NA (never fabricated):

# Residual within-observation correlation declaration:
rho12(pair_spec)
#> <residual correlation (rho12): 1 edge>
#>        y1       y2 term estimate std.error statistic p.value link predictors
#>  activity boldness <NA>       NA        NA        NA      NA tanh       temp
#>  constant
#>     FALSE
#> estimate NA: declaration only; fit the pair with drm_sem() or pass a bivariate
#> drmTMB fit to drm_psem().

# Higher-level between-individual random-effect correlation declaration:
corpairs(pair_spec)
#> <higher-level correlation (corpair): 1 edge>
#>  level       y1       y2 estimate std.error p.value
#>     id activity boldness       NA        NA      NA
#> estimate NA: declaration only; fit the pair with drm_sem() or pass a bivariate
#> drmTMB fit to drm_psem().

2. Moderation paths on the $\tanh$ link ( $x \to \rho_{12}$ )

In classical SEM and multivariate regression, residual correlation $\rho_{12}$ is assumed to be constant across all observations. In biological systems, however, environmental stress or temperature often modulates the coupling strength between traits—for example, coupling may tighten under high temperature and weaken under benign conditions.

drmSEM models $\rho_{12}$ as a first-class distributional component on the Fisher $z$ ( $\text{atanh}$ ) link:

$\text{atanh}(\rho_{12}) = \beta_0 + \beta_1 x$

$\rho_{12}(x) = \tanh(\beta_0 + \beta_1 x) = \frac{\exp(2(\beta_0 + \beta_1 x)) - 1}{\exp(2(\beta_0 + \beta_1 x)) + 1} \in (-1, 1)$

When declared via rho12 = ~ x, this edge represents a directed path into the residual correlation component. It appears in paths(sem) under component = "rho12", distinguishing it from both mean effects and $Y_1 \to Y_2$ regressions.

3. Joint fitting with `drm_sem()`

When passed to drm_sem(), a drm_pair() declaration is fitted as a single joint bivariate model in the drmTMB engine:

# Simulate bivariate location-scale-correlation data
n <- 400
temp <- rnorm(n)
# Temperature moderates residual correlation via tanh link:
# eta_rho = 0.20 + 0.50 * temp
eta_rho <- 0.20 + 0.50 * temp
rho_true <- 0.999 * tanh(eta_rho)

# Generate coupled residuals
e1 <- rnorm(n)
e2 <- rho_true * e1 + sqrt(pmax(1 - rho_true^2, 1e-6)) * rnorm(n)

dat <- data.frame(
  temp = temp,
  activity = 0.5 * temp + e1,
  boldness = 0.3 * temp + e2
)

# Downstream fitness response depending on both traits
dat$fitness <- 0.6 * dat$activity + 0.4 * dat$boldness + rnorm(n, sd = 0.5)

# Fit the SEM with joint bivariate node + downstream fitness node:
sem_biv <- drm_sem(
  pair = drm_pair(
    activity ~ temp,
    boldness ~ temp,
    rho12 = ~ temp
  ),
  fitness = drm_node(
    drmTMB::bf(fitness ~ activity + boldness),
    family = stats::gaussian()
  ),
  data = dat
)
#> ℹ Fitting bivariate pair "activity & boldness"
#> ✔ Fitting bivariate pair  [115ms]
#> 
#> ℹ Fitting node "fitness"
#> ✔ Fitting node "fitness" [72ms]
#> 

summary(sem_biv)
#> 
#> ── drmSEM summary ──────────────────────────────────────────────────────────────
#> 3 endogenous nodes; order "activity", "boldness", and "fitness"
#> 
#> ── Edges by component ──
#> 
#> mu: 2
#> mu1: 1
#> mu2: 1
#> rho12: 1
#> 
#> ── Paths ──
#> 
#> <drmSEM paths: 5 component-labelled coefficients>
#>      from       to component     link     term estimate std.error statistic
#>      temp activity       mu1 identity     temp   0.4699    0.0396    11.874
#>      temp activity     rho12     tanh     temp   0.5103    0.0452    11.291
#>      temp boldness       mu2 identity     temp   0.3390    0.0408     8.315
#>  activity  fitness        mu identity activity   0.5986    0.0212    28.201
#>  boldness  fitness        mu identity boldness   0.4311    0.0214    20.145
#>    p.value endogenous
#>   1.62e-32      FALSE
#>   1.46e-29      FALSE
#>   9.21e-17      FALSE
#>  5.72e-175       TRUE
#>   3.00e-90       TRUE

4. Post-estimation accessors

Once fitted, rho12() and paths() report the estimated correlation coefficients with standard errors, $z$ -statistics, and $p$ -values on the $\tanh$ link scale:

# Residual correlation parameter estimates on tanh link
rho12(sem_biv)
#> <residual correlation (rho12): 2 edges>
#>        y1       y2        term  estimate  std.error statistic      p.value link
#>  activity boldness (Intercept) 0.2125505 0.04556636  4.664636 3.091638e-06 tanh
#>  activity boldness        temp 0.5103055 0.04519613 11.290912 1.455185e-29 tanh
#>  predictors constant
#>        temp    FALSE
#>        temp    FALSE
#> rho12 coefficients on the tanh (atanh_guarded) link; not a y1 -> y2 path.

# Full structural paths table (includes directed paths to mu and to rho12)
paths(sem_biv)
#> <drmSEM paths: 5 component-labelled coefficients>
#>      from       to component     link     term estimate std.error statistic
#>      temp activity       mu1 identity     temp   0.4699    0.0396    11.874
#>      temp activity     rho12     tanh     temp   0.5103    0.0452    11.291
#>      temp boldness       mu2 identity     temp   0.3390    0.0408     8.315
#>  activity  fitness        mu identity activity   0.5986    0.0212    28.201
#>  boldness  fitness        mu identity boldness   0.4311    0.0214    20.145
#>    p.value endogenous
#>   1.62e-32      FALSE
#>   1.46e-29      FALSE
#>   9.21e-17      FALSE
#>  5.72e-175       TRUE
#>   3.00e-90       TRUE

5. Basis-set suppression in d-separation

In an ordinary piecewise SEM without a bivariate node, omitting an arrow between activity and boldness asserts that they are conditionally independent given temp ($Y_1 \perp\mkern-10mu\perp Y_2 \mid X$). This claim would enter the d-separation basis set and be tested via Fisher’s C (Shipley 2000).

When activity and boldness are declared as a drm_pair(), the residual covariance $\rho_{12}$ is explicitly estimated. Consequently:

The conditional independence claim activity _||_ boldness | temp is honestly suppressed from basis_set(sem).
Fisher’s C test tests only genuine conditional independencies outside the joint pair.

basis_set(sem_biv)
#>                                      claim    x       y              given
#> 1 temp _||_ fitness | {activity, boldness} temp fitness activity, boldness

Notice that activity versus boldness does not appear in the basis set.

Summary of best practices

Distinguish the three couplings: Use a directed arrow $Y_1 \to Y_2$ for causal causation; use drm_pair() for residual correlation ( $\rho_{12}$ ) and shared random effects (corpairs).
Model moderated coupling: Use rho12 = ~ x when environmental factors alter the tightness of correlation between traits.
Interpret link scales: Remember that $\rho_{12}$ regression coefficients operate on the $\text{atanh}$ link scale; transform back via $\tanh(\cdot)$ to obtain correlations on $[-1, 1]$ .
Inspect d-separation: Verify with basis_set(sem) that estimated covariance pairs are properly excluded from conditional independence tests.

References

Brooks, Mollie E., Kasper Kristensen, Koen J. van Benthem, et al. 2017. “glmmTMB Balances Speed and Flexibility Among Packages for Zero-Inflated Generalized Linear Mixed Models.” The R Journal 9 (2): 378–400. https://doi.org/10.32614/RJ-2017-066.

Lefcheck, Jonathan S. 2016. “piecewiseSEM: Piecewise Structural Equation Modelling in R for Ecology, Evolution, and Systematics.” Methods in Ecology and Evolution 7 (5): 573–79. https://doi.org/10.1111/2041-210X.12512.

Shipley, Bill. 2000. “A New Inferential Test for Path Models Based on Directed Acyclic Graphs.” Structural Equation Modeling 7 (2): 206–18. https://doi.org/10.1207/S15328007SEM0702_4.

Shipley, Bill. 2009. “Confirmatory Path Analysis in a Generalized Multilevel Context.” Ecology 90 (2): 363–68. https://doi.org/10.1890/08-1034.1.

Bivariate nodes: joint estimation, residual correlation, and moderation (drm_pair)