Extract cross-trait correlations

Returns the implied cross-trait correlations from a fit returned by gllvmTMB() at one or more covariance levels. Use the canonical input names "unit", "unit_obs", "phy", and "spatial"; legacy aliases "B", "W", and "spde" still work. The helper returns point estimates only by default, with no confidence level. Fisher-z and bootstrap intervals can be requested explicitly with the method argument; the former profile token is retained only to return a clear withdrawal error:

Usage

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

Arguments

fit

A fit returned by gllvmTMB. Julia bridge (engine = "julia") fits expose the ordinary unit tier only and carry no correlation-interval payload. Their off-diagonal rows have lower and upper set to NA, method = "none", and interval_status = "none". With entries = "all", self-pairs instead have method = "fixed". Use engine = "tmb" when you need correlation confidence intervals.

tier

Character vector. Use "all" (the default) to request every level present in the fit. Canonical inputs are "unit", "unit_obs", "phy", and "spatial"; legacy aliases "B", "W", and "spde" are accepted.

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. Numeric indices must be finite, distinct whole numbers within the trait range; fractional indices are rejected. Default NULL (all pairs).

level

Confidence level in (0, 1). Default 0.95.

method

One of "none" (default), "fisher-z", "wald" (alias of "fisher-z"), or "bootstrap". The accepted "profile" token is withdrawn and stops with an explanation. 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. This is a user-supplied sensitivity parameter, not an estimated effective sample size. It is consulted for method %in% c("fisher-z", "wald") and for explicit Wald fallback rows when a requested bootstrap tier is not yet routed. 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. These counts are only transparent heuristics: grouped, phylogenetic, spatial, latent, and non-Gaussian structure can make the iid Fisher variance inappropriate. If the automatic tier count is missing or below 4, Fisher-z rows return point correlations with unavailable interval bounds rather than substituting an arbitrary sample size. Set n_eff only when the analysis supplies a defensible target-specific rationale, and report that rationale.

nsim

Number of bootstrap replicates when method = "bootstrap". Default 500.

seed

Optional RNG seed for the bootstrap.

link_residual

How to treat the family-specific link-residual variance on the diagonal of the implied $\boldsymbol\Sigma$ before computing correlations:

"auto" (default): For non-Gaussian fits, add the link-specific implicit residual (e.g. $\pi^2/3$ for binomial-logit; $1$ for probit; trigamma() terms for Gamma / NB2 / Beta / etc.) to the diagonal before computing correlations. Returned correlations are on the latent-liability scale; this is the convention most readers expect. Gaussian fits are unaffected (link residual is $0$).
"none": Use the fitted model-implied $\Sigma$ directly with no link-residual addition. For ordinary latent() fits this includes the default diagonal Psi companion; correlations come out on the model-implied scale without the family adjustment.

Behaviour change in this release: the previous version hardcoded link_residual = "none". Non-Gaussian callers who relied on that behaviour will see different correlation values under the new default. A one-shot warning fires the first time per session that a non-Gaussian fit is processed without an explicit link_residual argument. Pass link_residual = "auto" to suppress the warning and lock the new behaviour, or link_residual = "none" to restore the previous behaviour.

entries

Table layout. "unique" (default) returns each unordered pair once, without the diagonal. "offdiag" includes both directions, without the diagonal. "all" includes both directions and the diagonal. A single tier with p traits has p(p-1)/2, p(p-1), or p squared rows, respectively. Non-unique layouts cannot be combined with pair.

Value

A data frame (tibble-like) with columns:

tier: Character level label. The current output stores internal labels "B", "W", "phy", and "spde"; use "unit" and "unit_obs" as input names in new calls.
trait_i, trait_j: Trait names. The default unique layout has i < j in fitted trait order.
correlation: Point estimate.
lower, upper: Confidence-interval bounds.
method: Method used to compute the CI.
interval_status: Claim-boundary marker: "none" for point-only output, "heuristic_unvalidated" for Fisher-z/Wald bounds, and "target_specific_uncalibrated" for bootstrap bounds. The bootstrap route computes intervals, but this function does not certify their frequentist coverage for the fitted target.

Details

"none" (default): point correlations only. No universal interval calibration has been established across covariance tiers, families, and targets.
"fisher-z": Fisher's z-transform heuristic interval. 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), but the classical iid-correlation variance is not a calibrated mixed-model standard error. Treat these bounds as a sensitivity display and see n_eff.
"profile": withdrawn. The nonlinear penalty-profile prototype did not yet provide a sufficiently strict constraint and constrained-optimizer diagnostic contract. This token stops with an explanation rather than returning bounds.
"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 a fitted model gives useful point estimates but Hessian- or profile-based intervals are not the right uncertainty summary. Structured tiers not yet resampled by bootstrap_Sigma, currently the SPDE spatial tier, return an explicit Wald/Fisher-z fallback with a message rather than fake bootstrap support.

For T traits at one tier, there are T(T-1)/2 unique correlations. A fit with T = 6 and four covariance levels present has up to 60 cross-trait correlations to report.

Pair table versus correlation matrix

By default, each unordered pair appears once, in fitted trait order. With $p$ traits, a single tier returns $p(p-1)/2$ rows; reversed pairs and diagonal entries are omitted. The first trait therefore appears only in trait_i, and the last only in trait_j. No trait is missing from the combined pair table when at least two traits are present. This also applies to ordinal responses.

Plotting these rows directly with geom_tile() leaves the reversed cells blank; alphabetical axis ordering can scatter those blanks. Use extract_correlations(fit, tier = "unit", entries = "all") for a complete table, or entries = "offdiag" to omit self-pairs. Mirrored rows retain the same estimates and interval columns; they are duplicate representations of one correlation, rather than independent estimates. Diagonal rows are fixed at 1, with method = "fixed", interval_status = "none", and NA interval bounds. Use entries = "unique" for interval/eye plots and pairwise analyses that should count each correlation once.

Use plot_correlations(cors, style = "heatmap") for a symmetric display, extract_Sigma(fit, level = "unit")$R for a square matrix, or extract_Sigma_table(fit, level = "unit", measure = "correlation", entries = "all") for a full table including the diagonal. Use the same link_residual convention when comparing these outputs.

For off-diagonal rows from an engine = "julia" bridge fit, lower/upper are NA, method = "none", and interval_status = "none".

Caveats

The former nonlinear penalty-profile route is withheld pending an exact constraint solver and explicit constrained-fit diagnostics.
Bootstrap uses bootstrap_Sigma refits and is the practical fallback when point estimates are useful but Hessian- or profile-based intervals are unavailable. Inspect bootstrap warnings, failed replicates, and interval width before treating the intervals as final.

Examples

## Simulate 80 units (for example, individual animals), each measured
## three times on four traits. Between units, the traits share one latent
## factor plus trait-specific variation, so they are correlated.
## Lambda holds each trait's loading on the shared factor; psi holds each
## trait's own between-unit variance.
set.seed(1)
Lambda <- matrix(c(0.9, 0.4, -0.3, 0.5), ncol = 1)
psi <- c(0.4, 0.3, 0.5, 0.2)
sim <- simulate_unit_trait(
  n_units = 80, n_obs_per_unit = 3, n_traits = 4,
  Lambda_B = Lambda, psi_B = psi
)
## One latent factor across units plus a per-trait variance (Psi).
fit <- gllvmTMB(
  value ~ 0 + trait + latent(0 + trait | unit, d = 1, unique = TRUE),
  data = sim$data, trait = "trait", unit = "unit"
)

## Default: point estimates of the between-unit correlations, one row per
## pair of traits ("B" in the tier column means the unit level).
cors <- extract_correlations(fit, tier = "unit")
cors
#>   tier trait_i trait_j correlation lower upper method interval_status
#> 1 unit trait_1 trait_2   0.4300549    NA    NA   none            none
#> 2 unit trait_1 trait_3  -0.1917988    NA    NA   none            none
#> 3 unit trait_2 trait_3  -0.1320519    NA    NA   none            none
#> 4 unit trait_1 trait_4   0.5097888    NA    NA   none            none
#> 5 unit trait_2 trait_4   0.3509855    NA    NA   none            none
#> 6 unit trait_3 trait_4  -0.1565349    NA    NA   none            none
## The correlations used to simulate the data; the estimates above differ
## from them by sampling error.
round(cov2cor(tcrossprod(Lambda) + diag(psi)), 2)
#>       [,1]  [,2]  [,3]  [,4]
#> [1,]  1.00  0.48 -0.32  0.61
#> [2,]  0.48  1.00 -0.23  0.44
#> [3,] -0.32 -0.23  1.00 -0.29
#> [4,]  0.61  0.44 -0.29  1.00

## A symmetric heatmap drawn from the pair table (needs ggplot2).
if (requireNamespace("ggplot2", quietly = TRUE)) {
  plot_correlations(cors, style = "heatmap")
}


## The same estimates as a full table (both directions plus the
## diagonal: 4 x 4 = 16 rows), as a square matrix, or as a table of
## matrix entries.
full <- extract_correlations(fit, tier = "unit", entries = "all")
nrow(full)
#> [1] 16
R <- extract_Sigma(fit, level = "unit")$R
round(R, 2)
#>         trait_1 trait_2 trait_3 trait_4
#> trait_1    1.00    0.43   -0.19    0.51
#> trait_2    0.43    1.00   -0.13    0.35
#> trait_3   -0.19   -0.13    1.00   -0.16
#> trait_4    0.51    0.35   -0.16    1.00
full_sigma <- extract_Sigma_table(fit, level = "unit",
                                  measure = "correlation", entries = "all")

## Fisher-z bounds are opt-in and depend on an effective sample size, n_eff.
## By default it is the number of units (80 here). We pass a more
## cautious 60 as an illustration, because each unit is seen only through
## three noisy measurements; in your own analysis, choose n_eff from your
## design and report why. These are rough bounds, not intervals
## calibrated to give 95% coverage in a mixed model.
cors_z <- extract_correlations(fit, tier = "unit", method = "fisher-z",
                               n_eff = 60L)
cors_z
#>   tier trait_i trait_j correlation      lower     upper   method
#> 1 unit trait_1 trait_2   0.4300549  0.1977218 0.6166414 fisher-z
#> 2 unit trait_1 trait_3  -0.1917988 -0.4250238 0.0653065 fisher-z
#> 3 unit trait_2 trait_3  -0.1320519 -0.3734542 0.1261011 fisher-z
#> 4 unit trait_1 trait_4   0.5097888  0.2939102 0.6761830 fisher-z
#> 5 unit trait_2 trait_4   0.3509855  0.1065576 0.5554099 fisher-z
#> 6 unit trait_3 trait_4  -0.1565349 -0.3947683 0.1014211 fisher-z
#>         interval_status
#> 1 heuristic_unvalidated
#> 2 heuristic_unvalidated
#> 3 heuristic_unvalidated
#> 4 heuristic_unvalidated
#> 5 heuristic_unvalidated
#> 6 heuristic_unvalidated

## A parametric bootstrap simulates new data from the fitted model and
## refits it each time. nsim = 40 keeps this example short, so expect a
## warning (it calls nsim `n_boot`) that so few resamples give noisy
## bounds; it stays quiet only from 999 resamples, and real analyses need
## at least several hundred. Bootstrap bounds are not certified to reach
## 95% coverage either.
# \donttest{
cors_b <- extract_correlations(fit, tier = "unit", method = "bootstrap",
                               nsim = 40, seed = 42)
#> Warning: `n_boot` = 40 gives noisy percentile bounds.
#> ℹ Percentile CIs are first-order accurate, so more draws buy precision, not
#>   correctness -- but below ~999 the Monte Carlo error in the endpoints is
#>   itself material.
#> ℹ `n_cores` >= 2 parallelises the refits; the cost here is wall-clock, not
#>   accuracy.
#> ℹ For Sigma diagonals, `profile_ci_total_variance()` needs no resampling at
#>   all.
cors_b
#>   tier trait_i trait_j correlation       lower       upper    method
#> 1 unit trait_1 trait_2   0.4300549  0.09281779  0.61561624 bootstrap
#> 2 unit trait_1 trait_3  -0.1917988 -0.38813278 -0.02809850 bootstrap
#> 3 unit trait_2 trait_3  -0.1320519 -0.50974747 -0.01567524 bootstrap
#> 4 unit trait_1 trait_4   0.5097888  0.10748158  0.70646699 bootstrap
#> 5 unit trait_2 trait_4   0.3509855  0.10887754  0.60104241 bootstrap
#> 6 unit trait_3 trait_4  -0.1565349 -0.42099639 -0.02646512 bootstrap
#>                interval_status
#> 1 target_specific_uncalibrated
#> 2 target_specific_uncalibrated
#> 3 target_specific_uncalibrated
#> 4 target_specific_uncalibrated
#> 5 target_specific_uncalibrated
#> 6 target_specific_uncalibrated
# }