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 (reporting decision 2C):

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) puts the nominal block on the observation scale via the \((\pi^2/6)(I+J)\) softmax residual, making multiple_r commensurable with single-scale partners; "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 single-number summary; read the per-pair latent correlations from extract_Sigma() / extract_correlations() instead. An ordinal-probit partner is not admitted with "auto": its latent residual is already fixed at 1 by the threshold model, so adding another unit residual would double-count it. 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

if (FALSE) { # \dontrun{
set.seed(12)
n_unit <- 120L
n_rep <- 3L
unit_id <- rep(seq_len(n_unit), each = n_rep)
score <- stats::rnorm(n_unit)
score_obs <- score[unit_id]
gaussian_y <- 0.8 * score_obs + stats::rnorm(length(unit_id), sd = 0.4)
category_y <- vapply(seq_along(unit_id), function(i) {
  eta <- c(0, 0.9 * score_obs[i], -0.6 * score_obs[i])
  prob <- exp(eta - max(eta))
  sample.int(3L, 1L, prob = prob / sum(prob))
}, integer(1L))
dat <- data.frame(
  unit = factor(c(unit_id, unit_id)),
  trait = factor(rep(c("g", "cat"), each = length(unit_id)),
                 levels = c("g", "cat")),
  family = factor(rep(c("g", "m"), each = length(unit_id)),
                  levels = c("g", "m")),
  value = c(gaussian_y, category_y)
)
families <- list(g = gaussian(), m = multinomial())
attr(families, "family_var") <- "family"
fit <- gllvmTMB(
  value ~ 0 + trait + latent(0 + trait | unit, d = 1, unique = TRUE),
  data = dat, family = families, trait = "trait", unit = "unit"
)
extract_cross_correlations(fit, contrasts = TRUE)
} # }