Phylogenetic-signal proportions per trait

For a phylogenetic species-level GLLVM, decomposes each trait's between-species latent variance into three additive components that sum to one: $$H_t^2 \;=\; \frac{[\boldsymbol\Sigma_\text{phy}]_{tt}}{V_{\eta,t}}, \qquad C^2_{\text{non},t} \;=\; \frac{[\boldsymbol\Sigma_\text{non,shared}]_{tt}}{V_{\eta,t}}, \qquad \psi_t \;=\; \frac{[\boldsymbol\Psi_\text{non}]_{tt}}{V_{\eta,t}},$$ where $V_{\eta,t} = [\boldsymbol\Sigma_\text{phy}]_{tt} + [\boldsymbol\Sigma_\text{non,shared}]_{tt} + [\boldsymbol\Psi_\text{non}]_{tt}$ is the total between-species latent variance for trait $t$ .

Usage

extract_phylo_signal(
  fit,
  ci = FALSE,
  conf_level = 0.95,
  method = c("profile", "wald", "bootstrap"),
  nsim = 500L,
  seed = NULL,
  link_residual = c("none", "auto")
)

Arguments

fit: A fit returned by gllvmTMB() with a phylo_latent() term.
ci: Logical. When TRUE, adds confidence-interval columns to the output for the H^2 column. Default FALSE for backward compatibility.
conf_level: Confidence level when ci = TRUE. Default 0.95.
method: One of "profile" (default), "wald", "bootstrap". Only used when ci = TRUE. For 2-component decompositions (a phylo diagonal vs species-level diagonal component only) profile uses a linear contrast; for 3-component decompositions (PGLLVM with phylo_latent plus a species-level latent decomposition with Psi) the full profile path is not yet implemented and falls back to numerical delta-method Wald bounds labelled "wald(numeric)".
nsim: Number of bootstrap replicates when method = "bootstrap". Default 500.
seed: Optional RNG seed for the bootstrap.
link_residual: "none" (default) keeps the historical species-level-latent denominator. "auto" adds each trait's distribution-specific link-scale residual $\sigma^2_d$ to the denominator and reports it as a fourth link_residual proportion column — the liability-scale phylogenetic heritability for threshold families (see the dedicated section below). Not yet compatible with ci = TRUE.

Value

A data frame with columns trait, H2, C2_non, Psi, V_eta (the denominator), one row per trait. The three proportions sum to 1.0 by construction. With link_residual = "auto" an extra link_residual proportion column is inserted before V_eta, the four proportions sum to 1.0, and V_eta includes $\sigma^2_d$. When ci = TRUE, three additional columns are added: H2_lower, H2_upper, H2_method.

Details

Interpretation:

$H_t^2$: phylogenetic signal — proportion of between- species latent variance attributable to phylogenetically structured variation ("evolutionary conservatism"). When the model uses the folded phylo_latent(..., unique = TRUE) decomposition, $\boldsymbol\Sigma_\text{phy}$ is the sum $\boldsymbol\Lambda_\text{phy} \boldsymbol\Lambda_\text{phy}^{\!\top} + \boldsymbol\Psi_\text{phy}$ and $H_t^2$ reflects the total phylogenetic variance.
$C^2_{\text{non},t}$: non-phylogenetic communality — proportion of variance attributable to shared non-phylogenetic axes ("coordinated tip-level lability" across traits).
$\psi_t$: uniqueness — proportion of variance not captured by any shared axis ("relative modularity").

Requires phylo_latent() (optionally with unique = TRUE for a phylogenetic Psi) plus a species-level latent() term, which carries its diagonal Psi companion by default. If the species-tier Psi is absent (a latent(..., unique = FALSE) subset), $\psi_t = 0$ for all traits and a cli::cli_inform() advisory fires.

Interval calibration

The point estimates are the supported claim. When ci = TRUE, the interval methods are provided for exploration: their empirical coverage is not certified for this estimand, so treat the intervals as exploratory rather than coverage-calibrated. See NEWS.md for the current coverage status.

Which variances enter the denominator

V_eta is the species-level latent variance: only components whose grouping is the cluster column contribute. Concretely $V_\eta = \sigma^2_{phy} + \sigma^2_{non}$, where $\sigma^2_{non}$ collects the species-grouped non-phylogenetic variance.

When unit == cluster (the usual unit = "species" setup) the ordinary latent() term and its Psi companion are species-level, so both enter.
When unit != cluster (a crossed site x species design) Lambda_B and its Psi are site-level and do not enter; the species-level non-phylogenetic variance is the cluster-tier q_it term (indep(0 + trait | <cluster>), reported as sd_q^2).

The denominator is defined at the species level. Other defensible definitions exist — for instance a total-variance denominator that also absorbs site-level variation, which would give a smaller H2 for the same fit. Before 2026-07-08 this function used the unit tier unconditionally, which silently reported H2 = 1 for every trait in a crossed design because the q_it variance was never read. Compare extract_proportions() if you want every component reported separately rather than folded into a species-level ratio.

Liability-scale heritability for categorical families

With the default link_residual = "none" the denominator contains only the estimated species-level components, so a phylogenetic-only fit reports H2 = 1 for every trait. For threshold/liability families the conventional phylogenetic heritability (Mizuno et al. 2025 J. Evol. Biol., eq 4/18/19; de Villemereuil & Nakagawa 2014) instead includes the distribution-specific latent residual $\sigma^2_d$ in the denominator. Request it with link_residual = "auto", which adds each trait's fixed link-scale residual (see extract_Sigma()'s "Family-aware link residuals" table) as a fourth component:

ordinal_probit() (and binomial probit): $\sigma^2_d = 1$ exactly, so a phylogenetic-only fit reports the paper's $H^2 = V_a / (V_a + 1)$ per trait.
multinomial(): each of the K-1 baseline-category contrasts is a logit, so $\sigma^2_d = \pi^2/3$ per contrast and a phylogenetic-only fit reports $H^2_{(k)} = V_{a(k)} / (V_{a(k)} + \pi^2/3)$ for contrast $k$. Contrast heritabilities are baseline-referenced: the softmax link residual is the full matrix $(\pi^2/6)(\mathbf{I} + \mathbf{J})$ — the shared baseline category couples the contrasts with $\pi^2/6$ off-diagonals — so each $H^2_{(k)}$ is the heritability of the "category k vs baseline" liability under that baseline choice, not of "the trait". Do not average or otherwise collapse the per-contrast values into a single scalar heritability; changing the baseline changes the contrasts.

When a categorical family (fid 14/16) is present and link_residual = "none", an advisory fires pointing at link_residual = "auto"; the default is unchanged for backward compatibility. ci = TRUE is not yet implemented for the liability-scale denominator and refuses with a typed error.

Structured trait-intercept rho fits are excluded: source strength is not phylogenetic variance share. Use extract_Sigma() and its source_strength metadata.

Examples

if (FALSE) { # \dontrun{
fit <- gllvmTMB(
  value ~ 0 + trait + phylo_latent(species, d = 2) +
                      latent(0 + trait | species, d = 2) +
                      indep(0 + trait | species),
  data       = df,
  trait      = "trait",
  unit       = "species",
  cluster    = "species",
  phylo_tree = tree
)
extract_phylo_signal(fit)
} # }