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. The default is not calibrated: it flags most datasets simulated from the paired model itself at realistic species counts, so a flag marks disagreement only.
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().
lrt: One-row data frame with the same columns as in compare_dep_vs_two_psi(), but the direction is reversed: the phylo_indep + indep refit is the paired fit with its shared loadings set to zero, so the paired fit is the larger model. chi2 is twice the log-likelihood gain of the paired fit over the refit and df its extra free parameters. A small p_value means the shared loadings improve the fit. The chi-square reference is approximate, because the loadings are not identified under the refit.
boundary: 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. A flag means the two fits disagree beyond threshold; it does not show that the paired model is wrong. The baseline carries its own sampling noise and the 0.10 default is not calibrated (on the issue 1571 simulations from the paired model, this check also flagged, with rel_disagreement of 0.24 and 0.16 for the two components at seed 18).

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.14510863 0.6001295        0.2417955 TRUE
#> 2 Sigma_non_diag 0.05749739 0.3689156        0.1558551 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.7537365 0.6644219 0.7495611
#> per_trait 0.5599785 0.6467232 0.5904570
## Likelihood-ratio test: does the paired fit's shared factor improve on
## the per-trait refit?
check$lrt
#>   loglik_joint loglik_alt     chi2 df    p_value
#> 1    -641.5846  -646.0889 9.008496  3 0.02917814