Cross-family correlations and environmental predictors

Can one environmental gradient explain a shared biological axis while the responses use different distributions? Yes: one ordinary latent() block can couple continuous, binary, count, ordered, and unordered traits, while lv = ~ x associates the shared axes with a unit-level predictor.

The long-data call is the canonical worked path. The unevaluated wide-data illustration below is a structural translation through the same gllvmTMB() entry point; this article does not demonstrate long/wide fit parity.

fit_long <- gllvmTMB(
  value ~ 0 + trait + latent(0 + trait | unit, d = 3, unique = FALSE, lv = ~ x),
  data = dat, family = family, trait = "trait", unit = "unit"
)
fit_wide <- gllvmTMB(
  traits(g, b, p, o, cat) ~ 1 + latent(1 | unit, d = 3, unique = FALSE, lv = ~ x),
  data = dat_wide,
  family = list(g = gaussian(), b = binomial(), p = poisson(),
                o = ordinal_probit(), cat = multinomial()),
  unit = "unit"
)

Traits do not arrive on one scale. A single organism carries continuous measurements, presence/absence marks, counts, ordered scores, and unordered categories at once; a single site records abundance, occupancy, and a dominant class. The question this article answers is narrow and concrete: can a Gaussian, a binary, a count, an ordinal, and a nominal trait sit on one shared latent model, can an observed gradient be associated with its axes, and can the model report the correlations between those families?

This page covers only the ordinary-unit cross-family route: one multinomial() trait and its non-nominal partners share latent(0 + trait | unit, d, lv = ~ x). One seeded realization provides a known-truth teaching check for the shared covariance and predictor effect; it is not repeated-simulation recovery evidence. Separate route-health canaries exercise every registered native family, but they are not recovery or interval studies. This page does not admit multiple nominal traits, structured tiers, augmented slopes, fixed X + X_lv, missing responses or LV predictors, REML, or calibrated arbitrary-composition intervals. The five-family latent covariance is valid, but the one-number residual-augmented latent-liability (model-scale) nominal summary is shown on a non-ordinal subset because the current automatic residual convention is not coherent for an ordinal-probit partner.

Why put different families on one latent scale

The motivation is the same one that drives joint species distribution models and phenotypic-integration analyses: the interesting quantity is not each response on its own but how they covary. If a species that is more abundant is also more likely to be present, and its dominant habitat class shifts with both, that joint structure is the signal. A shared latent factor $\mathbf z_i$ is the device that couples the families: every response loads on the same axes, so a correlation between families is well-defined even though the families live on different observation scales.

The alternative — fit each response separately and correlate the point predictions — discards the measurement model and gives correlations that are biased by each family’s link and dispersion. Putting the families on one latent scale keeps the coupling in the model where it belongs.

The model

For unit $i$ , predictors $\mathbf M_i$ shift the mean of a $d$ -dimensional latent score, while $\mathbf e_i$ records remaining unit-to-unit innovation:

$\mathbf z_i = \mathbf M_i\boldsymbol\alpha + \mathbf e_i, \qquad \mathbf e_i \sim \mathcal N(\mathbf 0, \mathbf I_d).$

Each pseudo-trait $t$ uses its own response family and link:

$\eta_{it} = \mu_t + \boldsymbol\lambda_t^{\!\top}\mathbf z_i, \qquad y_{it} \mid \mathbf z_i \sim \text{family}_t(\,g_t^{-1}(\eta_{it})\,).$

Gaussian, binomial, and Poisson traits each occupy a single row of the loading matrix $\boldsymbol\Lambda$ . An ordinal_probit() trait occupies one row and crosses a set of ordered thresholds. A multinomial() trait is different: an unordered response with $K$ categories is not one latent dimension but $K-1$ of them — the baseline-category contrasts — so it expands into $K-1$ rows of $\boldsymbol\Lambda$ , named "<trait>:<cat>". The shared low-rank covariance among all pseudo-traits is $\boldsymbol\Lambda\boldsymbol\Lambda^{\!\top}$ , and the cross-family correlations are its standardised off-diagonals — rotation-invariant, and the quantity to report.

The predictor effect on each trait’s linear-predictor scale is also rotation-invariant:

$\mathbf B_{lv} = \boldsymbol\Lambda\boldsymbol\alpha^\top.$

Raw alpha and Lambda entries depend on the chosen latent-axis orientation; B_lv and Lambda %*% t(Lambda) are the cross-fit scientific targets.

latent() can carry a between-unit trait-specific variance $\boldsymbol\Psi$ (unique = TRUE), so the full tier covariance is $\boldsymbol\Sigma = \boldsymbol\Lambda\boldsymbol\Lambda^\top + \boldsymbol\Psi$ — $\boldsymbol\Psi$ is the variance each trait has among units beyond the shared factor. For a nominal contrast with one categorical draw per unit, the corresponding between-unit $\boldsymbol\Psi$ is not identified. Replication can identify that variance in principle, but the current engine conservatively maps it off even for the replicated design used here. Admitted partner traits retain fitted diagonal companions only when their physical loading, latent-predictor, and engine-free diagonal parameters pass a necessary joint mean-covariance dimension screen. Passing that screen does not itself establish stable recovery or calibrated uncertainty. This rank-3 teaching model requests the loadings-only form with unique = FALSE so its focus stays on the shared covariance; the cross-family correlations assessed below live in $\boldsymbol\Lambda\boldsymbol\Lambda^\top$ .

A five-family worked example

We simulate 500 units, each observed 6 times, with five traits sharing one 3-dimensional factor. A standardised environmental gradient shifts the mean of all three axes; unit innovations create the remaining correlation:

trait	family	value column holds
`g`	`gaussian()`	a real number
`b`	`binomial()`	0 / 1
`p`	`poisson()`	a count
`o`	`ordinal_probit()`	an ordered code `1..Ko`
`cat`	`multinomial()`	a category code `1..K`

The truth is a single loading matrix over the six pseudo-traits g, b, p, o, cat:2, cat:3 (the nominal trait cat contributes the last two rows). Everything the model must recover is contained in $\boldsymbol\Lambda\boldsymbol\Lambda^{\!\top}$ .

library(gllvmTMB)

Simulate

set.seed(20260718)

N     <- 500L   # units
reps  <- 6L     # observations per unit (replication sharpens recovery)
d     <- 3L     # shared latent dimensions
units <- seq_len(N)

# True loadings for the six pseudo-traits on the three shared axes.
# The nominal trait `cat` (K = 3) contributes the rows cat:2 and cat:3.
Lambda <- rbind(
  g       = c( 1.2,  0.2,  0.0),
  b       = c( 0.9,  0.6,  0.1),
  p       = c( 0.3,  1.0,  0.4),
  o       = c( 0.7,  0.4,  0.9),
  `cat:2` = c( 1.0,  0.5,  0.2),
  `cat:3` = c(-0.5,  0.8,  0.6)
)

# Trait means / baseline intercepts (ordinal handled by thresholds instead).
mu  <- c(g = 0, b = 0, p = log(2), o = 0, `cat:2` = 0, `cat:3` = 0)
tau <- c(-1.6, -0.2, 1.2)   # 3 thresholds -> Ko = 4 ordered categories

# One predictor-informed latent position per unit, shared across replicates.
x <- as.numeric(scale(seq(-1, 1, length.out = N)))
alpha <- c(0.55, -0.30, 0.20)
E <- matrix(rnorm(N * d), N, d)
U_mean <- outer(x, alpha)
U_total <- U_mean + E
Eta <- U_total %*% t(Lambda)           # N x 6 unit-level predictors
colnames(Eta) <- rownames(Lambda)
B_lv_true <- drop(Lambda %*% alpha)

sim_one_rep <- function(rep_id) {
  g_val <- mu["g"] + Eta[, "g"] + rnorm(N)                      # Gaussian
  b_val <- rbinom(N, 1L, plogis(mu["b"] + Eta[, "b"]))         # binary
  p_val <- rpois(N, exp(mu["p"] + Eta[, "p"]))                  # count
  o_val <- findInterval(mu["o"] + Eta[, "o"] + rnorm(N), tau) + 1L  # ordinal
  e2 <- mu["cat:2"] + Eta[, "cat:2"]                            # nominal
  e3 <- mu["cat:3"] + Eta[, "cat:3"]                            #   contrasts
  P  <- cbind(1, exp(e2), exp(e3)); P <- P / rowSums(P)         # softmax
  cat_val <- apply(P, 1L, function(pr) sample.int(3L, 1L, prob = pr))
  rbind(
    data.frame(unit = units, replicate = rep_id, trait = "g",   family = "g", x = x, value = g_val),
    data.frame(unit = units, replicate = rep_id, trait = "b",   family = "b", x = x, value = b_val),
    data.frame(unit = units, replicate = rep_id, trait = "p",   family = "p", x = x, value = p_val),
    data.frame(unit = units, replicate = rep_id, trait = "o",   family = "o", x = x, value = o_val),
    data.frame(unit = units, replicate = rep_id, trait = "cat", family = "m", x = x, value = cat_val)
  )
}

dat <- do.call(rbind, lapply(seq_len(reps), sim_one_rep))
dat$unit   <- factor(dat$unit)
dat$trait  <- factor(dat$trait,  levels = c("g", "b", "p", "o", "cat"))
dat$family <- factor(dat$family, levels = c("g", "b", "p", "o", "m"))

table(dat$trait)
#> 
#>    g    b    p    o  cat 
#> 3000 3000 3000 3000 3000

dat_wide <- reshape(
  dat[c("unit", "replicate", "x", "trait", "value")],
  idvar = c("unit", "replicate", "x"), timevar = "trait", direction = "wide"
)
names(dat_wide) <- sub("^value\\.", "", names(dat_wide))

The family column maps each row to its response family; the nominal trait is the single column cat, whose values are category codes 1..K. It expands into the pseudo-traits cat:2 and cat:3 inside the fit, not in the data.

Fit

One gllvmTMB() call carries all five families on the shared factor. The family argument is a named list keyed by the family column’s levels, and attr(family, "family_var") names that column.

family <- list(
  g = gaussian(), b = binomial(), p = poisson(),
  o = ordinal_probit(), m = multinomial()
)
attr(family, "family_var") <- "family"

fit <- gllvmTMB(
  value ~ 0 + trait + latent(0 + trait | unit, d = 3, unique = FALSE, lv = ~ x),
  data   = dat,
  family = family,
  trait  = "trait",
  unit   = "unit",
  silent = TRUE
)

c(convergence = fit$opt$convergence, pdHess = isTRUE(fit$sd_report$pdHess))
#> convergence      pdHess 
#>           0           1

The explicit unique = FALSE requests the loadings-only covariance. The same shape with automatic $\Psi$ passes the parser’s necessary joint mean-covariance dimension screen, but that screen is not recovery or interval evidence. We keep the teaching fit loadings-only because the scientific target here is the shared cross-family covariance.

Recover the predictor effect

extract_lv_effects(type = "trait_effect") returns $B_{lv}$ rather than raw axis coefficients. Each row is an association with x on that response’s linear-predictor scale.

B_hat <- extract_lv_effects(fit, type = "trait_effect")
B_compare <- data.frame(
  trait = B_hat$trait,
  true = unname(B_lv_true[B_hat$trait]),
  recovered = B_hat$estimate
)
B_compare
#>   trait   true   recovered
#> 1     g  0.600  0.59548297
#> 2     b  0.335  0.35899423
#> 3     p -0.055 -0.07988727
#> 4     o  0.445  0.42960278
#> 5 cat:2  0.440  0.50470563
#> 6 cat:3 -0.395 -0.40649681

For Gaussian traits this is an associated mean change. For Poisson traits it is a log expected-count change; exponentiate it for a multiplicative contrast. For ordinal and nominal rows it remains a latent linear-predictor association. The model is associational: changing x is not interpreted as an intervention.

Recover the cross-family correlations

extract_Sigma() with part = "shared" returns the latent low-rank covariance $\boldsymbol\Lambda\boldsymbol\Lambda^{\!\top}$ over all six pseudo-traits. We standardise it and compare against the truth built from Lambda.

Sigma_true <- Lambda %*% t(Lambda)
R_true <- cov2cor(Sigma_true)

Sig_hat <- extract_Sigma(fit, level = "unit", part = "shared", link_residual = "none")
R_hat   <- cov2cor(Sig_hat$Sigma)

# align the recovered matrix to the truth's ordering by name
ord   <- rownames(R_true)
R_hat <- R_hat[ord, ord]

pairs <- rbind(
  c("g", "b"), c("g", "p"), c("g", "o"),
  c("g", "cat:2"), c("g", "cat:3"),
  c("b", "cat:2"), c("p", "cat:3"), c("o", "cat:2")
)

data.frame(
  pair      = paste(pairs[, 1], pairs[, 2], sep = " ~ "),
  true      = round(mapply(function(i, j) R_true[i, j], pairs[, 1], pairs[, 2]), 3),
  recovered = round(mapply(function(i, j) R_hat[i, j], pairs[, 1], pairs[, 2]), 3),
  row.names = NULL
)
#>        pair   true recovered
#> 1     g ~ b  0.908     0.887
#> 2     g ~ p  0.412     0.403
#> 3     g ~ o  0.626     0.627
#> 4 g ~ cat:2  0.941     0.954
#> 5 g ~ cat:3 -0.323    -0.309
#> 6 b ~ cat:2  0.989     0.976
#> 7 p ~ cat:3  0.712     0.713
#> 8 o ~ cat:2  0.787     0.797

In this seeded realization, the estimated correlations track the known truth across all five families, including the two nominal contrasts cat:2 and cat:3 and their sign — a strong positive g ~ cat:2, a negative g ~ cat:3. Differences here can reflect sampling, estimation, and optimization. They do not establish repeated-simulation recovery or guarantee that increasing the number of units or replicates will reduce a discrepancy. Every off-diagonal pair, not just the eight highlighted here, is available in R_hat.

offd <- upper.tri(R_true)
plot(R_true[offd], R_hat[offd],
     xlab = "true correlation", ylab = "recovered correlation",
     xlim = c(-1, 1), ylim = c(-1, 1), pch = 19, col = "#3366aa")
abline(0, 1, lty = 2, col = "grey40")

Scatterplot of estimated versus known off-diagonal latent correlations for one seeded realization. A dashed diagonal identity line runs from minus one to one; points close to it have similar estimated and known correlations in this realization.

Known versus estimated latent cross-family correlations in one seeded realization, for all off-diagonal pairs. The line is the identity.

Summarising the nominal trait: the reference-invariant multiple correlation

The pairwise table above lists a separate correlation for each nominal contrast (cat:2, cat:3). Often a reader wants one number for “how strongly is the categorical trait associated with each partner?” that does not depend on which category was chosen as the reference. That is what extract_cross_correlations() reports: a reference-invariant multiple correlation between the categorical trait and each single-scale partner, optionally with the underlying $(K-1)$ -vector.

The five-family fit contains an ordinal-probit partner. Its threshold model already fixes the latent residual at one, so an "auto" summary would add that unit residual a second time. extract_cross_correlations() therefore refuses that combination rather than attenuating the reported association. We retain the valid five-family latent-covariance result above, then refit the same data without the ordinal trait for the residual-augmented latent-liability (model-scale) one-number summary.

ordinal_refusal <- tryCatch(
  extract_cross_correlations(fit, level = "unit", link_residual = "auto"),
  error = function(e) cli::ansi_strip(conditionMessage(e))
)
ordinal_refusal
#> [1] "The residual-augmented model-scale nominal summary is unavailable with\nan ordinal (probit/logit) partner.\nℹ Ordinal partner traits: o.\nℹ The threshold model already fixes each ordinal latent residual by\n  construction (1 for ordinal_probit, pi^2/3 for ordinal_logit); `link_residual\n  = \"auto\"` would add it again.\n→ Use a supported non-ordinal partner set for `extract_cross_correlations()`,\n  or report latent-scale pairwise quantities with `extract_Sigma(..., part =\n  \"shared\", link_residual = \"none\")`."

dat_summary <- droplevels(dat[dat$trait != "o", ])
family_summary <- family[c("g", "b", "p", "m")]
attr(family_summary, "family_var") <- "family"

fit_summary <- gllvmTMB(
  value ~ 0 + trait + latent(0 + trait | unit, d = 3, unique = FALSE, lv = ~ x),
  data   = dat_summary,
  family = family_summary,
  trait  = "trait",
  unit   = "unit",
  silent = TRUE
)

xc <- extract_cross_correlations(
  fit_summary, level = "unit", link_residual = "auto"
)
xc
#>   nominal partner multiple_r
#> 1     cat       g  0.6200802
#> 2     cat       b  0.2737271
#> 3     cat       p  0.3925946

For wide data, the non-ordinal summary subset has the following unevaluated structural translation. It is an illustration only: this article fits the summary subset in long form and does not demonstrate long/wide parity.

dat_summary_wide <- dat_wide[c("unit", "replicate", "x", "g", "b", "p", "cat")]
fit_summary_wide <- gllvmTMB(
  traits(g, b, p, cat) ~
    1 + latent(1 | unit, d = 3, unique = FALSE, lv = ~ x),
  data = dat_summary_wide,
  family = list(g = gaussian(), b = binomial(), p = poisson(),
                cat = multinomial()),
  unit = "unit",
  silent = TRUE
)

The single-number summary is a residual-augmented latent-liability (model-scale) association, not a correlation of observed responses. The "auto" route adds the applicable family-specific residual for each admitted partner and the full softmax residual matrix for the nominal block. This is the deliberate default, and the reason is in the next section.

Honest caveats

These are limits of the capability as it stands at 0.6, not incidental notes.

Fitted $\boldsymbol\Psi$ and the fixed link residual are different objects. The between-unit $\boldsymbol\Psi$ inside latent() is an estimated among-unit diagonal companion. The extractor’s $R_{\mathrm{link}}$ is a separate family-specific residual addition on the latent-liability/model scale; for the nominal block it is $(\pi^2/6)(\mathbf I+\mathbf J)$ . With one categorical draw per unit the nominal contrast $\boldsymbol\Psi$ is not identified. Replication can identify it in principle, but the current engine conservatively maps it off even in the replicated design here; admitted partner diagonals remain fitted where the model identifies them. A replication-aware multinomial refinement is future work. Either way the pairwise latent correlations (the matrix) are the quantities compared with known truth in this teaching check. See the vocabulary of $\boldsymbol\Psi$ and covariance and correlation.
The single-number summary is residual-augmented latent-liability (model) scale, not observed-response scale. For an admitted non-ordinal partner set, the reference-invariant multiple correlation uses each applicable family-specific residual and the softmax residual $(\pi^2/6)(\mathbf I + \mathbf J)$ — the covariance of the differenced Gumbel errors in the random-utility representation of the softmax (McFadden 1974), which reduces to the binomial’s $\pi^2/3$ when $K = 2$ . On the latent scale (link_residual = "none") the multiple correlation degenerates toward 1 whenever the $K-1$ contrast block spans the partner, because a $(K-1)$ -dimensional block can align with almost any single direction. That is why the residual-augmented model-scale association is the sensible default for the one-number summary, while the per-pair latent correlations remain the right object for the matrix. Ordinal-probit is not admitted on this automatic summary route because its residual is already fixed by the threshold model; the typed refusal above prevents double counting.
Nominal-trait covariance can be weakly informed. A single categorical draw carries little information about a $(K-1)$ -dimensional liability, and a reduced-rank estimate can rail at $\pm 1$ . This article’s 500 units $\times$ 6 replicates provide one seeded teaching realization, not repeated-simulation evidence about the effects of more units or replication. Do not extrapolate its agreement with known truth to one observation per unit or to a new data-generating regime.
Intervals are cell-specific, not universal. Existing evidence calibrates selected rank-1 Gaussian-anchor B_lv cells. It does not calibrate the rank-3 arbitrary composition shown here or cross-family correlation intervals. Treat these outputs as point/recovery-oriented.
Gaussian and lognormal residual scales are separate when they coexist. Gaussian residual SD uses raw response units; lognormal residual SD uses log(y). A joint fit therefore reports two sigma_eps slots, shared within family rather than across these incomparable scales. Pure fits retain the historical one-slot contract.
One multinomial trait per fit. The current implementation supports a single multinomial() trait in a fit. Two nominal responses on one shared factor is not yet supported.