Extract the implied trait covariance / correlation at one tier

Implements the decomposition $$\boldsymbol\Sigma_\text{tier} \;=\; \underbrace{\boldsymbol\Lambda_\text{tier}\boldsymbol\Lambda_\text{tier}^\top}_{\text{shared (latent)}} \;+\; \underbrace{\boldsymbol\Psi_\text{tier}}_{\text{unique}},$$ where ordinary latent() now carries both $\boldsymbol\Lambda$ and the diagonal $\boldsymbol\Psi$ companion by default. This is the same decomposition the behavioural-syndromes / phenotypic-integration literature uses (Bartholomew et al. 2011).

Usage

extract_Sigma(
  fit,
  level = c("unit", "unit_slope", "unit_obs", "phy", "phy_slope", "spatial",
    "spde_slope", "cluster", "cluster2", "B", "B_slope", "W", "spde"),
  part = c("total", "shared", "unique", "psi"),
  link_residual = c("auto", "none"),
  .skip_warn = FALSE
)

Arguments

fit: A fit returned by gllvmTMB(). Admitted engine = "julia" bridge fits expose the ordinary unit tier only: link_residual = "none" reconstructs $\Lambda\Lambda^\top$ from retained loadings, while link_residual = "auto" uses the retained GLLVM.jl residual-augmented payload where available and keeps Gaussian / lognormal rows on the native no-op convention. unit_obs, structured tiers, and augmented-slope tiers remain gated for Julia bridge extractors.
level: One of "unit" (between-unit), "unit_obs" (within-unit), "phy" (phylogenetic), "spatial", or "cluster". Legacy aliases "B", "W", and "spde" are accepted with a soft-deprecation message.
part: One of "total" (default), "shared", "unique". "psi" is an alias for "unique" (the per-trait residual Psi diagonal, now folded into latent() by default).
link_residual: For non-Gaussian fits. "auto" (default) adds a per-trait link-specific implicit residual variance to the diagonal of Sigma, giving the marginal latent-scale interpretation; in mixed- family fits each unambiguous trait gets the residual implied by its family/link (see "Family-aware link residuals" below for the full table). "none" returns the fitted model covariance without link-residual additions ($\boldsymbol\Lambda \boldsymbol\Lambda^\top + \boldsymbol\Psi$ where a latent decomposition is present). For Gaussian or lognormal-only fits this argument is effectively a no-op (their implied $\sigma^2_d = 0$).
.skip_warn: Internal flag (default FALSE). When TRUE, suppresses the once-per-session deprecation message that the internal .normalise_level() helper emits for legacy level aliases ("B", "W", "spde", "total"). Used by internal callers (e.g. extract_Omega()) that have already issued the deprecation message themselves; not part of the public API.

Value

For part = "total" or "shared": a list with components Sigma (T x T matrix), R (T x T correlation matrix; only for "total"), level, part, and note (character vector of advisory messages, e.g. about a no-Psi unique = FALSE fit).

For part = "unique": a list with s (length-T named numeric vector of unique variances), level, part, note.

For a phylo_dep(1 + x1 + ... + xs | species) fit with one or more slopes, call with level = "phy": the result is the single full unstructured (1+s)T x (1+s)T covariance over the trait-stacked (intercept, slope_1, ..., slope_s) random-effect columns – a list with Sigma and R carrying INTERLEAVED dimnames (per trait: intercept.<t>, then slope.<t> for s == 1 or slope.<x_j>.<t> for s >= 2), level = "phy_dep", part = "dep", and a note. The part and link_residual arguments do not apply to this single unstructured block and are ignored. (The unit / unit_obs tiers return NULL for a dep-only fit, as it carries no between/within-unit covariance term.)

Current phylo_indep(1 + x | species) uses the same interleaved 2T report channel with all cross-trait blocks fixed to zero. Its return has level = "phy_indep_slope" and part = "indep"; each trait retains its own intercept-slope block. This is distinct from the legacy shared 2 x 2 augmented phylo_unique channel.

For augmented spatial slopes, call with level = "spatial". spatial_unique(1 + x | coords) retains the shared 2 x 2 field-scale channel and returns level = "spde_base_slope". Current spatial_indep(1 + x | coords) returns the interleaved, cross-trait block-diagonal 2T x 2T field covariance with level = "spde_indep_slope" and part = "indep"; spatial_dep returns the full unstructured block with level = "spde_dep" and part = "dep". These SPDE matrices are on the fitted field-covariance scale; the returned note records the marginal-scale conversion.

For an ordinary individual-level random-regression fit with latent(1 + x | unit, d = K), indep(1 + x | unit), or their long-form equivalents, call with level = "unit_slope": the result is the augmented 2T x 2T covariance over trait-specific intercept and slope coefficients, with row names intercept.<trait> and slope.<x>.<trait>. As for level = "unit", part = "shared" returns Lambda_aug Lambda_aug^T, part = "unique" returns the augmented diagonal, and part = "total" returns their sum.

Details

When a fit has no Psi component

If the formula deliberately uses latent(0 + trait | unit, d = K, unique = FALSE), the engine fits only the $\boldsymbol\Lambda \boldsymbol\Lambda^\top$ component. Calling extract_Sigma(fit, level, part = "total") then returns just the shared component. This is useful for an intentionally loadings-only target. It should not be described as a full latent-liability covariance unless the relevant family/link residual is also included. Whether the diagonal is incomplete depends on the estimand and response family.

For Gaussian / lognormal / Gamma fits this function emits an advisory note when a reduced-rank tier has no Psi component. Use the ordinary latent(..., unique = TRUE) default for $\boldsymbol\Lambda\boldsymbol\Lambda^\top + \boldsymbol\Psi$.

For non-Gaussian families (binomial, Poisson, Gamma) the latent-scale residual variance has a closed-form approximation that should be added to the diagonal of $\boldsymbol\Sigma$ – see the link_residual argument below.

The `part` argument

"total" (default): $\boldsymbol\Sigma_\text{tier} = \boldsymbol\Lambda \boldsymbol\Lambda^\top + \boldsymbol\Psi$ – the matrix users almost always want for reporting correlations.
"shared": $\boldsymbol\Lambda \boldsymbol\Lambda^\top$ only – the reduced-rank, rotation-invariant component. Diagonals are $\sum_\ell \Lambda_{t\ell}^2$; these are not the trait variances, they are the shared part of the trait variances.
"unique": $\boldsymbol\Psi_\text{tier}$ only – the trait-specific unique variances, returned as a length-T named numeric vector (the diagonal of $\boldsymbol\Psi$).

Caveat: `"shared"` vs `"unique"` partition is only weakly identified

The total $\boldsymbol\Sigma_\text{tier} = \boldsymbol\Lambda \boldsymbol\Lambda^\top + \boldsymbol\Psi$ is rotation-invariant, so part = "total" avoids arbitrary loading rotation. Rotation invariance does not guarantee that the data identify every covariance entry; fit health, sample design, and uncertainty still matter. The split between $\boldsymbol\Lambda \boldsymbol\Lambda^\top$ and $\boldsymbol\Psi$ is only weakly identified – different optimiser starts can flow trait $t$'s variance more into the shared ("shared") or unique ("unique") component, with the same total likelihood. In a synthetic Poisson + OLRE recovery run with target $\sigma^2_S = (0.5, 0.4, 0.3, 0.6)$, gllvmTMB returned $(0, 0.67, 0, 0.57)$ – total trait variances correct, but the partition arbitrary.

This is the standard rotation-and-shift indeterminacy of factor models (gllvm and Hmsc have it too); communality $c_t^2$ computed from a single fit therefore inherits the same indeterminacy. For interpretation, lean on part = "total" for correlations and use extract_communality() only when paired with bootstrap_Sigma() or the variance-decomposition consistency checks in gllvmTMB_diagnose().

Family-aware link residuals

For non-Gaussian responses each row carries an implicit observation-level residual on the latent (link) scale. Adding it to the per-trait diagonal of $\boldsymbol\Sigma$ defines one family-specific latent-scale convention. It can make the denominator explicit, but it does not make different response families directly commensurate without additional scientific assumptions.

Per-family formulas (Nakagawa & Schielzeth 2010; Nakagawa, Johnson & Schielzeth 2017):

`gaussian` (identity)	$\sigma^2_d = 0$
`binomial(link = "logit")`	$\sigma^2_d = \pi^2/3 \approx 3.290$
`binomial(link = "probit")`	$\sigma^2_d = 1$
`binomial(link = "cloglog")`	$\sigma^2_d = \pi^2/6 \approx 1.645$
`poisson(link = "log")`	$\sigma^2_d = \log(1 + 1/\hat\mu_t)$ (lognormal-Poisson approx.)
`lognormal(link = "log")`	$\sigma^2_d = 0$ (sigma_eps already models the log-scale residual)
`Gamma(link = "log")`	$\sigma^2_d = \psi'(\hat\nu)$ where $\hat\nu = \hat\phi_{\gamma,t}$ is the per-trait shape
`nbinom2(link = "log")`	$\sigma^2_d = \psi'(\hat\phi)$ where $\hat\phi$ is the per-trait NB2 dispersion
`tweedie(link = "log")`	$\sigma^2_d = \log(1 + \hat\phi \hat\mu_t^{\hat p - 2})$ (delta method)
`Beta(link = "logit")`	$\sigma^2_d = \psi'(\hat\mu_t \hat\phi) + \psi'((1 - \hat\mu_t)\hat\phi)$ (Smithson & Verkuilen 2006)
`betabinomial(link = "logit")`	$\sigma^2_d = \pi^2/3 + \psi'(\hat\mu_t \hat\phi) + \psi'((1 - \hat\mu_t)\hat\phi)$

multinomial() is the exception to the scalar-per-trait table: its distribution-specific residual is not a scalar but the $(K-1)\times(K-1)$ matrix $(\pi^2/6)(\mathbf{I} + \mathbf{J})$ over the category-contrast pseudo-traits – $\pi^2/3$ on the diagonal (each baseline contrast is a logit, exactly as for binomial(link = "logit")) and $\pi^2/6$ off-diagonal (the shared baseline category couples the contrasts). This is the residual covariance implied by the softmax's random-utility (Gumbel) representation (McFadden 1974), and it reduces to the binomial $\pi^2/3$ at $K = 2$. (It is distinct from MCMCglmm's arbitrary fixed identification residual $(1/K)(\mathbf{I}+\mathbf{J})$ and that package's MCMC-specific $c^2$ overdispersion correction, neither of which a direct softmax fit needs.) link_residual = "auto" (the default) applies this full-matrix residual: the $(\pi^2/6)(\mathbf{I}+\mathbf{J})$ block is added to the multinomial trait's $(K-1)$-contrast block of $\boldsymbol\Sigma$ (the diagonal $\pi^2/3$ entry-by-entry with the other families, the $\pi^2/6$ off-diagonal as a block); "none" returns the latent-scale V unchanged. The addition is tier-agnostic – it applies to whichever tier $\boldsymbol\Sigma$ is requested (level = "phy", a spatial or kernel tier, or an ordinary latent term) – and it makes the categorical block commensurable with any single-scale (e.g. Gaussian) trait sharing the fit. Relatedly, a diagonal among-category V is not independence: the independence null for the contrasts is itself $(\mathbf{I}+\mathbf{J})$-structured (equal diagonal, equal off-diagonal).

For mixed-family fits the residual is computed per trait when that trait has one unambiguous family/link, then added to the diagonal of $\boldsymbol\Sigma$ entry-by-entry. The default link_residual = "auto" applies this; "none" returns the fitted model covariance without link-residual additions ($\boldsymbol\Lambda \boldsymbol\Lambda^\top + \boldsymbol\Psi$ where a latent decomposition is present). For continuous-only Gaussian or lognormal fits "auto" is a no-op (their per-trait $\sigma^2_d$ is zero).

Citations: Nakagawa & Schielzeth (2010); Nakagawa, Johnson & Schielzeth (2017); McFadden (1974) for the multinomial-logit random-utility residual — see References below.

Future: 3+ latent tiers

The engine currently supports two latent() tiers ("unit" and "unit_obs"); legacy level = "B" and "W" aliases are still accepted with a soft-deprecation message. level = "phy" extracts the phylogenetic implied $\boldsymbol\Sigma_\text{phy}$ from phylo_latent(). If a future release adds 3+ latent tiers, level = "<colname>" will dispatch to the corresponding tier without API change. For now, custom strings error with a clear unsupported-level message.

References

Nakagawa, S. & Schielzeth, H. (2010). Repeatability for Gaussian and non-Gaussian data: a practical guide for biologists. Biological Reviews 85, 935-956. doi:10.1111/j.1469-185X.2010.00141.x

Nakagawa, S., Johnson, P. C. D., & Schielzeth, H. (2017). The coefficient of determination $R^2$ and intra-class correlation coefficient from generalized linear mixed-effects models revisited and expanded. Journal of the Royal Society Interface 14(134), 20170213. doi:10.1098/rsif.2017.0213

McFadden, D. (1974). Conditional logit analysis of qualitative choice behavior. In P. Zarembka (Ed.), Frontiers in Econometrics (pp. 105-142). Academic Press, New York. (Random-utility origin of the multinomial-logit $(\pi^2/6)(\mathbf{I}+\mathbf{J})$ contrast residual.)

Examples

if (FALSE) { # \dontrun{
fit <- gllvmTMB(
  value ~ 0 + trait +
          latent(0 + trait | unit, d = 2),
  data  = df,
  trait = "trait",
  unit  = "unit"
)
extract_Sigma(fit, level = "unit", part = "total")$Sigma   # full T x T cov
extract_Sigma(fit, level = "unit", part = "shared")$Sigma  # rr-only
extract_Sigma(fit, level = "unit", part = "unique")$s      # diag(s_unit)
} # }

`gaussian` (identity)	\(\sigma^2_d = 0\)
`binomial(link = "logit")`	\(\sigma^2_d = \pi^2/3 \approx 3.290\)
`binomial(link = "probit")`	\(\sigma^2_d = 1\)
`binomial(link = "cloglog")`	\(\sigma^2_d = \pi^2/6 \approx 1.645\)
`poisson(link = "log")`	\(\sigma^2_d = \log(1 + 1/\hat\mu_t)\) (lognormal-Poisson approx.)
`lognormal(link = "log")`	\(\sigma^2_d = 0\) (sigma_eps already models the log-scale residual)
`Gamma(link = "log")`	\(\sigma^2_d = \psi'(\hat\nu)\) where \(\hat\nu = \hat\phi_{\gamma,t}\) is the per-trait shape
`nbinom2(link = "log")`	\(\sigma^2_d = \psi'(\hat\phi)\) where \(\hat\phi\) is the per-trait NB2 dispersion
`tweedie(link = "log")`	\(\sigma^2_d = \log(1 + \hat\phi \hat\mu_t^{\hat p - 2})\) (delta method)
`Beta(link = "logit")`	\(\sigma^2_d = \psi'(\hat\mu_t \hat\phi) + \psi'((1 - \hat\mu_t)\hat\phi)\) (Smithson & Verkuilen 2006)
`betabinomial(link = "logit")`	\(\sigma^2_d = \pi^2/3 + \psi'(\hat\mu_t \hat\phi) + \psi'((1 - \hat\mu_t)\hat\phi)\)