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). A flag means the two fits disagree; it does not show that the paired two-psi model is wrong or not identified.

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. 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: 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.
lrt: One-row data frame: the likelihood-ratio test of the paired fit against the unstructured refit, which contains it. Columns loglik_joint, loglik_alt, chi2 (twice the log-likelihood gain of the refit), df (the refit's extra free parameters) and p_value (upper tail of a chi-square with df degrees of freedom). All NA when the refit failed, when fit_two_psi has no non-phylogenetic term (the refit always adds one, so the fits are not nested) or when a log-likelihood is unavailable. The chi-square reference is approximate because the paired model sits on a reduced-rank and variance boundary of the unstructured one; use p_value as a guide, after checking alt_fit$opt$convergence. rel_disagreement and flag stay as descriptive summaries; the likelihood ratio is the comparison that accounts for the refit's extra flexibility.
boundary: List with character vectors joint and alt: the boundary flags of the paired fit and of the refit (from gllvmTMB_diagnose(), for example "near_zero_sd_B"). A variance estimated at zero can still report pdHess = TRUE. When either vector is non-empty the chi-square reference of the lrt is unreliable, because the null value lies on the edge of the parameter space; read p_value as a rough guide only. The function also prints a one-line note in that case.
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. Agreement between the two fits is reassuring. Disagreement is a prompt to look further, not a verdict: the more flexible baseline estimates many more parameters (for example non-phylogenetic off-diagonals whose true value may be 0) and carries its own sampling noise at realistic species counts. In simulations from the paired model itself (T = 3, d = 1; see issue 1571), the default threshold flagged nearly every dataset: 24/24 seeds at 64 species, 4/4 at 128 and 3/4 at 256, with the paired fit usually closer to the true non-phylogenetic Sigma than the baseline. The 0.10 default of threshold is therefore not calibrated; treat a flag as a comparison to inspect, for example by checking the boundary and convergence diagnostics of both fits.

The lrt component is the likelihood-based comparison: because the paired fit is nested in the unstructured refit, it reports their likelihood-ratio test, which accounts for the refit's many extra parameters. rel_disagreement and flag are kept as descriptive summaries only.

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 what is compared (see the caveat above).
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)
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.1810554 0.3928673        0.4608563 TRUE
#> 2 Sigma_non 0.0933282 0.2285131        0.4084152 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 
## The likelihood-ratio test of the two fits (the unrestricted refit
## contains the paired fit); this accounts for its extra parameters.
check$lrt
#>   loglik_joint loglik_alt     chi2 df  p_value
#> 1    -641.5846  -639.7612 3.646863  3 0.302208
check$joint$Sigma_phy  # paired fit
#>            t1         t2         t3
#> t1  0.7537365  0.2841576 -0.4519319
#> t2  0.2841576  0.6644219 -0.2607492
#> t3 -0.4519319 -0.2607492  0.7495611
check$dep$Sigma_phy    # unrestricted refit
#>             t1          t2         t3
#> t1  0.66630576  0.03233858 -0.2236047
#> t2  0.03233858  0.69424559 -0.1763252
#> t3 -0.22360468 -0.17632525  0.5467096