Cross-trait correlations with confidence intervals

Returns the implied cross-trait correlations from a fitted gllvmTMB_multi model at one or more tiers (between-unit B, within-unit W, phylogenetic phy, spatial spde), with 95% (or other-level) confidence intervals. Four method names are supported via the method argument:

Usage

extract_correlations(
  fit,
  tier = "all",
  pair = NULL,
  level = 0.95,
  method = c("fisher-z", "profile", "wald", "bootstrap"),
  n_eff = NULL,
  nsim = 500L,
  seed = NULL
)

Arguments

fit: A gllvmTMB_multi fit returned by gllvmTMB.
tier: Character vector. Subset of c("B", "W", "phy", "spde"). Use "all" (the default) to request every tier present in the fit.
pair: Optional length-2 character or integer vector specifying one trait pair (c("trait_1", "trait_2") or c(1, 2)). When supplied, only that pair is returned for each requested tier. Default NULL (all pairs).
level: Confidence level in (0, 1). Default 0.95.
method: One of "fisher-z" (default), "profile", "wald" (alias of "fisher-z"), or "bootstrap". See Details.
n_eff: Optional positive integer (>= 4): override the effective sample size used in Fisher's $\widehat{\mathrm{SE}}(\hat z) = 1/\sqrt{n_{\text{eff}} - 3}$ formula. Only consulted for method %in% c("fisher-z", "wald"). The default heuristic uses fit$n_sites for tier "B" ("unit"), fit$n_site_species for "W" ("unit_obs"), fit$n_species for "phy", and fit$n_sites elsewhere — adequate for well-identified Gaussian fits, but can under-cover when latent structure or non-Gaussian likelihood reduces the effective N below the unit count. Set n_eff explicitly for non-trivial structure (e.g. n_eff = fit$n_species for a binomial joint-SDM where species mediates the cross-trait correlation).
nsim: Number of bootstrap replicates when method = "bootstrap". Default 500.
seed: Optional RNG seed for the bootstrap.

Value

A data frame (tibble-like) with columns:

tier: Character: "B", "W", "phy", or "spde".
trait_i, trait_j: Trait names with i < j.
correlation: Point estimate.
lower, upper: Confidence-interval bounds.
method: Method used to compute the CI.

Details

"fisher-z" (default): Fisher's z-transform Wald CI. Computes $\hat z = \mathrm{atanh}(\hat\rho)$, $\widehat{\mathrm{SE}}(\hat z) = 1/\sqrt{n_{\text{eff}} - 3}$, constructs the CI on z, then back-transforms via $\tanh(\cdot)$ (so bounds are guaranteed inside $[-1, 1]$). Fast (seconds for any T) and reasonable for most fits. See n_eff for the effective-sample-size override.
"profile": profile-likelihood CI via fix-and-refit. Most accurate for skewed sampling distributions, but slow ($O(T(T-1)/2)$ constrained refits per call). Use for final, publication-grade CIs on small T.
"wald": backward-compat alias of "fisher-z" with the same numerics. Emits a one-shot inform pointing at the canonical name. Kept for scripts that filter on method == "wald".
"bootstrap": parametric bootstrap via bootstrap_Sigma. Slowest (full sampling distribution); use when you need full uncertainty propagation to multiple downstream quantities.

For T traits at one tier, there are T(T-1)/2 unique correlations. A multi-trait fit at 5 tiers (B, W, phy, non, spde) with T = 6 has up to 75 cross-trait correlations to report.

Caveats

For tiers with latent() only and small ranks, the profile path can be unstable due to rotation indeterminacy of the factor model (the implied Σ is identifiable but the split into Lambda Lambda^T and S is not). Fall back to method = "bootstrap" when the profile fails.
Bootstrap uses bootstrap_Sigma which conditions on the fitted random effects (parametric simulation), so the bootstrap CIs reflect residual-level uncertainty rather than full posterior uncertainty in the variance components. For full posterior CIs use a Bayesian fit (e.g. rstanarm / brms).

Examples

if (FALSE) { # \dontrun{
set.seed(1)
s <- simulate_site_trait(
  n_sites = 80, n_species = 6, n_traits = 4,
  mean_species_per_site = 4,
  Lambda_B = matrix(c(0.9, 0.4, -0.3, 0.5), 4, 1),
  S_B = c(0.4, 0.3, 0.5, 0.2)
)
fit <- gllvmTMB(
  value ~ 0 + trait + latent(0 + trait | site, d = 1) +
          unique(0 + trait | site),
  data = s$data
)
## Default Fisher-z Wald CIs (fast, bounded inside \eqn{[-1, 1]}).
cors <- extract_correlations(fit, tier = "unit")
## With n_eff override (e.g. for binomial joint-SDM where species
## count is the relevant effective N).
cors2 <- extract_correlations(fit, tier = "unit", n_eff = 6L)
## Profile-likelihood (slow but accurate for skewed sampling).
cors_p <- extract_correlations(fit, tier = "unit", method = "profile")
## Bootstrap (B = 200).
cors_b <- extract_correlations(fit, tier = "unit", method = "bootstrap",
                               nsim = 200, seed = 42)
} # }