Cross-family correlations between a nominal (multinomial) trait and partners

For a fit where a multinomial() trait shares a latent factor with other-family traits, report the association between the nominal trait and each partner trait. A nominal trait spans K-1 baseline-category contrasts, so its association with a single-scale partner is a vector, summarized two ways:

multiple_r: the reference-invariant multiple correlation \(R = \sqrt{\Sigma_{pc}\,\Sigma_{cc}^{-1}\,\Sigma_{cp}/\sigma_{pp}}\) between partner p and the whole K-1 contrast block c. Invariant to the baseline category; magnitude in [0, 1].
contrast_r (when contrasts = TRUE): the (K-1)-vector of individual contrast correlations, labelled by category-vs-baseline.

Usage

extract_cross_correlations(
  fit,
  level = "unit",
  contrasts = FALSE,
  link_residual = c("auto", "none"),
  method = c("point", "bootstrap", "profile", "wald"),
  conf = 0.95,
  nsim = 500L,
  seed = NULL
)

Arguments

fit: a fitted gllvmTMB_multi with one multinomial() trait and a shared ordinary unit tier (latent(0 + trait | unit, d)).
level: covariance tier. Only the ordinary "unit" tier is admitted.
contrasts: if TRUE, also return the per-contrast (K-1)-vector (as a list column).
link_residual: passed to extract_Sigma(); "auto" (default) adds the \((\pi^2/6)(I+J)\) softmax residual to form a residual-augmented latent-liability (model-scale) association, not a correlation of observed responses; it applies the documented residual convention to the nominal block and each admitted partner; "none" uses the latent (loadings) scale. Note: with "none" and a loadings-only fit whose shared factor the contrast block spans (e.g. the block dimension >= the partner's), multiple_r degenerates to 1 (the partner is then a deterministic function of the block). This is correct but uninformative, which is why "auto" (full-rank via the residual) is the sensible default for the model-scale single-number summary; read the per-pair latent correlations from extract_Sigma() / extract_correlations() instead. An ordinal (probit or logit) partner is not admitted with "auto" for this one-number summary. This is a scope limit, not a double count: "auto" adds the ordinal liability residual (1 for probit, pi^2/3 for logit) exactly once. Use a supported non-ordinal partner set for the one-number summary, or report latent-scale pairwise quantities with extract_Sigma(..., part = "shared", link_residual = "none").
method: one of "point" (default; point estimates only, no interval columns), "wald" (fast Fisher-z intervals on BOTH multiple_r and contrast_r; no refits or root-finding; heuristic bounds can be unavailable when the effective sample size is too small or an estimate is on a boundary), "bootstrap" (parametric-bootstrap interval on the aggregate multiple_r via bootstrap_Sigma()), or "profile". The profile token is retained for compatibility but is withdrawn: it always stops with a typed error rather than returning an interval. Use point estimates, "wald", or "bootstrap" instead. None of the available intervals is coverage-calibrated yet — they are recovery-oriented; "wald" is the most approximate (multiple_r is not a Pearson correlation). Use them to explore + report bugs, not for inference.
conf: confidence level for the interval columns. Default 0.95.
nsim: number of parametric-bootstrap replicates when method = "bootstrap". Default 500.
seed: optional RNG seed for the bootstrap.

Value

a data.frame, one row per (nominal, partner) pair. With method = "point" the columns are nominal, partner, multiple_r (and contrast_r when contrasts = TRUE). With method = "wald", "bootstrap", per-estimand interval columns are added: multiple_r_lower/multiple_r_upper/multiple_r_method/multiple_r_interval_status (scalar), and when contrasts = TRUE contrast_r_lower/contrast_r_upper (list columns) plus contrast_r_method/contrast_r_interval_status. The estimand NOT served by the requested method keeps point-only interval columns (NA bounds, method = "point"). Computed intervals carry an interval_status flag ("heuristic_unvalidated" for "wald", "target_specific_uncalibrated" for "bootstrap") – coverage is not yet certified for any route (recovery evidence for this route is not yet complete).

Examples

## Simulate 120 animals, each watched on three occasions. On each occasion
## we record a body-condition score (continuous) and the behaviour seen
## (1 = resting, 2 = foraging, 3 = travelling). One shared animal-level
## score raises condition and the odds of foraging and lowers the odds of
## travelling; condition also varies between animals for other reasons.
set.seed(12)
n_unit <- 120L
n_rep <- 3L
unit_id <- rep(seq_len(n_unit), each = n_rep)
score <- rnorm(n_unit)[unit_id]
own <- rnorm(n_unit, sd = 0.5)[unit_id]
condition <- 0.8 * score + own + rnorm(length(unit_id), sd = 0.4)
behaviour <- vapply(score, function(s) {
  ## log-odds of each behaviour relative to resting
  p <- exp(c(0, 0.9 * s, -0.6 * s))
  sample.int(3L, 1L, prob = p / sum(p))
}, integer(1L))
dat <- data.frame(
  unit   = factor(rep(unit_id, 2)),
  trait  = factor(rep(c("condition", "behaviour"), each = length(unit_id)),
                  levels = c("condition", "behaviour")),
  family = rep(c("gaussian", "multinomial"), each = length(unit_id)),
  value  = c(condition, behaviour)
)
## The `family` column picks the response family for each row.
families <- list(gaussian = gaussian(), multinomial = multinomial())
## One shared factor across animals. With unique = TRUE, condition also
## gets its own between-animal variance; the behaviour categories do not
## (the fit may print a note saying so), which matches the simulation.
fit <- gllvmTMB(
  value ~ 0 + trait + latent(0 + trait | unit, d = 1, unique = TRUE),
  data = dat, family = families, trait = "trait", unit = "unit"
)
#> ℹ Skipping the default between-unit Psi for 2 binary / categorical-contrast /
#>   single-observation ordinal traits under the family-specific identifiability
#>   gate.
#> ℹ Affected traits: behaviour:2, behaviour:3.
#> ℹ For multinomial contrast traits behaviour:2, behaviour:3, replication can
#>   identify a contrast-specific diagonal in principle, but the current engine
#>   conservatively maps it off and rejects explicit multinomial
#>   `unique()`/`indep()` terms. Use the shared `latent()` block for admitted
#>   cross-family covariance.
#> • Mapped `theta_diag_B[t]` and the corresponding `s_B` row off.
#> This message is displayed once per session.
## Both summaries are between-animal correlations on the model's latent
## scale, after adding the residual variance of the categorical
## (softmax) link, as the printed message reports; they are not
## correlations between the raw observations.
## multiple_r (0 to 1): how strongly condition goes with behaviour as a
## whole; it does not depend on which category is the baseline.
## contrast_r: the correlation of condition with the log-odds of each
## other category against the baseline (category 1, resting); here
## positive for foraging and negative for travelling, as simulated.
cc <- extract_cross_correlations(fit, contrasts = TRUE)
#> Added per-trait link-implicit residual variance to diag(Sigma): - condition
#> (gaussian): 0.000 - behaviour:2 (multinomial): 3.290 - behaviour:3
#> (multinomial): 3.290
#> Added the multinomial softmax off-diagonal link residual pi^2/6 ~= 1.645 within
#> each (K-1)-contrast block (with the pi^2/3 diagonal above, the full block
#> residual is (pi^2/6)(I + J); McFadden 1974).
cc
#>     nominal   partner multiple_r   contrast_r
#> 1 behaviour condition  0.4914992 0.281302....
cc$contrast_r[[1]]
#> behaviour:2 behaviour:3 
#>   0.2813021  -0.2941398