Cheap diagonal cross-check for the paired phylogenetic decomposition (large T)

Refits the user's data with phylo_indep + indep – a stacked multivariate generalisation of T univariate phylogenetic GLMMs in the spirit of Williams et al. (2025) Eq. 3 – using the same engine, family, link, and grouping as the supplied two-psi fit. Compares the per-trait diagonals of the joint two-psi fit's implied $\boldsymbol\Sigma_{\mathrm{phy}}$ (and $\boldsymbol\Sigma_{\mathrm{non}}$) against the marginal-only baseline. Tests only the diagonals (cheap); does not test the off-diagonal cross-trait covariances.

Usage

compare_indep_vs_two_psi(
  fit_two_psi,
  threshold = 0.1,
  phylo_vcv = NULL,
  phylo_tree = NULL
)

Arguments

fit_two_psi: A joint two-psi fit returned by gllvmTMB(), e.g. one produced by gllvmTMB(value ~ 0 + trait + phylo_latent(species, d = K_phy, unique = TRUE) + indep(0 + trait | species), ...). The cross-check refits the same data and family with phylo_dep + dep and compares.
threshold: Numeric (default 0.10): relative-disagreement threshold (Frobenius RMSE / Frobenius magnitude of the unstructured estimate) above which a component is flagged.
phylo_vcv, phylo_tree: Optional. The phylogenetic correlation matrix or ape::phylo tree, only needed if fit_two_psi was produced by an older package version that did not store the phylogeny on the fit. Default is to recover them from the fit (fit_two_psi$phylo_vcv / fit_two_psi$phylo_tree).

Value

A list with components:

joint: Joint two-psi fit's implied diag(Sigma_phy) and diag(Sigma_non) (named numeric vectors of length T).
indep: The marginal-only baseline's diag(Sigma_phy) and diag(Sigma_non) (same shape; may contain NULL if alt fit failed).
agreement: Data frame with rows for Sigma_phy_diag and Sigma_non_diag, columns rmse, indep_mag, rel_disagreement, and flag.
flag, threshold, alt_fit: As in compare_dep_vs_two_psi().

Details

Use this diagnostic when T is large (T $\geq$ 30) and the unstructured fit needed by compare_dep_vs_two_psi() is intractable. Otherwise prefer compare_dep_vs_two_psi() (gold standard).

Cross-check intent: the per-trait diagonals are the simplest summary the joint two-psi fit and the marginal baseline both target. If the diagonals already disagree beyond the threshold, no inspection of the cross-trait off-diagonals is warranted.

References

Williams, M. J., McGillycuddy, M., Drobniak, S. M., Bolker, B. M., Warton, D. I., & Nakagawa, S. (2025). Fast phylogenetic generalised linear mixed-effects modelling using the glmmTMB R package. bioRxiv 2025.12.20.695312. doi:10.64898/2025.12.20.695312

Examples

set.seed(18)
## A balanced phylogeny of 64 species, built in base R: six rounds of
## evenly spaced splits. Two species' phylogenetic correlation is the
## share of the tree's depth they spent as one lineage (5/6 for sister
## species, 0 for species on opposite sides of the root).
n_sp <- 64
Cphy <- matrix(0, n_sp, n_sp)
for (level in 1:6) {
  clade <- (seq_len(n_sp) - 1) %/% 2^(6 - level)
  Cphy <- Cphy + outer(clade, clade, "==") / 6
}
sp <- paste0("sp", seq_len(n_sp))
dimnames(Cphy) <- list(sp, sp)

## Three traits. A species' value for each trait is the sum of a shared
## phylogenetic factor (with a trait-specific loading), a trait-specific
## phylogenetic effect and a non-phylogenetic species effect. Three
## individuals are measured per species, so these species effects can be
## told apart from measurement noise.
L <- t(chol(Cphy))  # turns independent draws into tree-correlated ones
phy <- outer(drop(L %*% rnorm(n_sp)), c(0.8, 0.6, -0.5)) +
  L %*% matrix(rnorm(n_sp * 3, sd = 0.6), n_sp, 3)
non <- matrix(rnorm(n_sp * 3, sd = 0.5), n_sp, 3)
df <- expand.grid(trait = c("t1", "t2", "t3"), individual = 1:3,
                  species = sp)
idx <- cbind(match(df$species, sp), as.integer(df$trait))
df$value <- phy[idx] + non[idx] + rnorm(nrow(df), sd = 0.5)

## The paired ("two-psi") fit: the phylogenetic covariance is one shared
## factor plus trait-specific phylogenetic variances; the
## non-phylogenetic species variances are independent across traits.
fit <- gllvmTMB(
  value ~ 0 + trait +
    phylo_latent(species, d = 1, unique = TRUE, vcv = Cphy) +
    indep(0 + trait | species),
  data = df, trait = "trait", unit = "species", cluster = "species"
)

## Refit with per-trait variances only (phylo_indep + indep) and compare
## the two fits' per-trait phylogenetic and non-phylogenetic variances.
## This cheap check is meant for fits with many traits, where the
## unrestricted refit of compare_dep_vs_two_psi() is too slow; three
## traits keep the example fast.
check <- compare_indep_vs_two_psi(fit)
check$alt_fit$opt$convergence  # 0 means the refit converged
#> [1] 0
## rel_disagreement is the gap between the two fits relative to the size
## of the per-trait estimates; flag marks gaps above `threshold`
## (0.10 by default).
check$agreement
#>        component       rmse indep_mag rel_disagreement flag
#> 1 Sigma_phy_diag 0.14510652 0.6001295        0.2417920 TRUE
#> 2 Sigma_non_diag 0.05749641 0.3689156        0.1558525 TRUE
check$flag
#> [1] TRUE
## Both components are flagged here, even though the data were simulated
## from the paired model itself: with 64 species each trait's
## phylogenetic and non-phylogenetic variances are only loosely
## determined, and the paired fit (which also uses the covariances among
## traits) and this per-trait refit land some way apart. A flag means the
## two fits disagree; on its own it does not show that the paired model
## is wrong. The phylogenetic variances of the two fits:
rbind(two_psi = check$joint$Sigma_phy_diag,
      per_trait = check$indep$Sigma_phy_diag)
#>                  t1        t2        t3
#> two_psi   0.7537366 0.6644161 0.7495559
#> per_trait 0.5599785 0.6467232 0.5904570