Canonical likelihood-based cross-check for the paired phylogenetic decomposition

Refits the user's data with full unstructured T x T phylogenetic and non-phylogenetic trait covariances (phylo_dep + dep) using the same engine, family, link, and unit / cluster grouping as the supplied two-psi fit. Compares the joint two-psi fit's implied $\boldsymbol\Sigma_{\mathrm{phy}} = \boldsymbol\Lambda_{\mathrm{phy}}\boldsymbol\Lambda_{\mathrm{phy}}^\top + \boldsymbol\Psi_{\mathrm{phy}}$ (and the analogous $\boldsymbol\Sigma_{\mathrm{non}}$) against the unstructured baseline. Per-component RMSE plus a flag field (TRUE when any component disagrees beyond threshold) identifies two-psi identifiability failures.

Usage

compare_dep_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: The two-psi fit's implied Sigma_phy and Sigma_non (T x T matrices).
dep: The unstructured phylo_dep + dep baseline's implied Sigma_phy and Sigma_non (T x T matrices). May contain NULL entries if the alt fit failed.
agreement: Data frame with rows for Sigma_phy and Sigma_non, columns rmse (Frobenius RMSE between joint and dep), dep_mag (Frobenius magnitude of the dep estimate), rel_disagreement (rmse / dep_mag), and flag (rel_disagreement > threshold, or NA when the alt refit failed so the component could not be tested).
flag: Logical: TRUE if any component is flagged, FALSE if all tested components agreed, NA if nothing could be tested (distinguishes agreement from a failed refit).
threshold: The threshold used.
alt_fit: The refit gllvmTMB_multi object (or NULL on failure), retained so users can inspect convergence and pull other extractor outputs.

Details

Conceptual basis: Williams et al. (2025) bioRxiv 2025.12.20.695312 Eq. 3, generalised across traits. The unstructured fit has $T(T+1)/2$ free parameters per tier and is the strongest available benchmark for the implied total covariance. When the joint two-psi fit is well-identified, its implied total Sigma matches the unstructured baseline; when it isn't, the diagnostic flags the disagreement and the user knows to relax the rank or revisit identifiability.

The cross-check separates three questions:

Level 1: can we fit?: Both fits converged.
Level 2: total Sigma_phy and Sigma_non agree?: This diagnostic answers Level 2 directly. The two estimators target the SAME total covariance with different parameterisations, so agreement is the identifiability check.
Level 3: split into Lambda Lambda^T + Psi?: Diagnostic does not test Level 3 directly; Level 2 agreement is the pre-requisite for Level 3 to be meaningful.

Computational scope: tractable for T <= ~30. For larger T, prefer compare_indep_vs_two_psi().

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

Hadfield, J. D. & Nakagawa, S. (2010). General quantitative genetic methods for comparative biology: phylogenies, taxonomies and multi- trait models for continuous and categorical characters. Journal of Evolutionary Biology 23, 494-508. doi:10.1111/j.1420-9101.2009.01915.x

Meyer, K. & Kirkpatrick, M. (2008). Perils of parsimony: properties of reduced-rank estimates of genetic covariance matrices. Genetics 180, 1153-1166. doi:10.1534/genetics.108.090159

Felsenstein, J. (2005). Using the quantitative genetic threshold model for inferences between and within species. Philosophical Transactions of the Royal Society B 360, 1427-1434. doi:10.1098/rstb.2005.1669

Felsenstein, J. (2012). A comparative method for both discrete and continuous characters using the threshold model. American Naturalist 179, 145-156. doi:10.1086/663681

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 unrestricted phylogenetic and non-phylogenetic covariances
## (phylo_dep + dep) and compare the two fits' covariance matrices.
check <- compare_dep_vs_two_psi(fit)
#> ! `phylo_vcv = ...` as a global argument to `gllvmTMB()` is deprecated.
#> ℹ Pass `vcv = ...` inside the phylo_*() keyword instead, e.g.
#>   `phylo_scalar(species, vcv = Cphy)`.
#> → Prefer `tree = ...` where you have a phylo object: it triggers the
#>   sparse-A^-1 path (~24x faster at n_species = 1000).
#> This message is displayed once per session.
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 unrestricted estimate; flag marks gaps above `threshold`
## (0.10 by default).
check$agreement
#>   component      rmse   dep_mag rel_disagreement flag
#> 1 Sigma_phy 0.1810531 0.3928663        0.4608517 TRUE
#> 2 Sigma_non 0.0933278 0.2285130        0.4084136 TRUE
check$flag
#> [1] TRUE
## Both components are flagged here, even though the data were simulated
## from the paired model itself. The unrestricted refit gains little
## log-likelihood for its three extra parameters, yet it splits the trait
## covariance differently between the phylogenetic and non-phylogenetic
## parts: with 64 species the data only loosely determine that split.
## So a flag means the two fits disagree; on its own it does not show
## that the paired model is wrong.
c(two_psi = logLik(fit), unrestricted = logLik(check$alt_fit))
#>      two_psi unrestricted 
#>    -641.5846    -639.7612 
check$joint$Sigma_phy  # paired fit
#>            t1         t2         t3
#> t1  0.7537366  0.2841587 -0.4519279
#> t2  0.2841587  0.6644161 -0.2607528
#> t3 -0.4519279 -0.2607528  0.7495559
check$dep$Sigma_phy    # unrestricted refit
#>             t1          t2         t3
#> t1  0.66629780  0.03234647 -0.2236005
#> t2  0.03234647  0.69424786 -0.1763350
#> t3 -0.22360047 -0.17633500  0.5467058