Extract cross-lineage coevolutionary modules — extract_coevolution

extract_coevolution_modules() standardizes a component-specific cross-lineage covariance block and decomposes it into coupled trait axes. For a named coevolution kernel tier, the helper computes R = Sigma_row^{-1/2} Gamma Sigma_col^{-1/2} from the covariance returned by extract_Sigma(), then applies a singular-value decomposition. Gamma is the cross-lineage block of the shared covariance. Sigma_row and Sigma_col are the within-lineage blocks of the TOTAL covariance (part = "total": shared Lambda Lambda^T plus the per-trait Psi diagonal when the tier carries one). Standardising by the shared blocks alone would give singular values of exactly 1 whenever d is at most the number of traits per lineage, so they would carry no information.

This is a point-estimate derived-output helper, covered by the package's validation tests, for fitted dense-kernel coevolution models whose component-specific Gamma blocks are already extractable with extract_Gamma(). It does not report uncertainty, choose the biological rank, estimate rho, or calibrate null thresholds. Use simulation, bootstrap, and kernel-separability checks before treating the returned axes as scientific evidence.

Usage

extract_coevolution_modules(
  fit,
  level,
  row_traits,
  col_traits,
  scale = c("shape", "effect"),
  n_modules = NULL,
  tol = sqrt(.Machine$double.eps)
)

Arguments

fit: A fitted gllvmTMB_multi object.
level: Character scalar naming the fitted covariance tier.
row_traits, col_traits: Character vectors naming the row-lineage and column-lineage traits.
scale: Character scalar. "shape" (default) uses Gamma_shape = Lambda_row Lambda_col^T; "effect" uses the fixed-rho Gamma_effect available for kernels built by make_cross_kernel().
n_modules: Optional positive integer limiting the number of returned singular axes. By default all axes are returned.
tol: Numerical tolerance for the generalized inverse square roots of the within-lineage total covariance blocks.

Value

A list with R, modules, row_axes, and col_axes. R is the standardized cross-lineage correlation-like block. modules has one row per singular axis with its singular value and squared-value share. The singular values lie in [0, 1] and are informative only when the within-lineage total covariance is larger than the shared part, that is, when the tier has a Psi component estimated above zero (for example kernel_latent(..., unique = TRUE)). For a latent-only tier with d at most the traits per lineage they are 1 by construction, and the trait loadings are the informative output. row_axes and col_axes contain trait loadings for each coupled axis.

Examples

set.seed(1)
## A small host-partner system: 12 host and 12 partner species. Each
## lineage has its own similarity matrix. Here similarity decays with
## distance along a species ordering, a stand-in for a phylogenetic
## correlation matrix such as ape::vcv(tree, corr = TRUE).
host <- paste0("H", 1:12)
partner <- paste0("P", 1:12)
A_H <- 0.6^abs(outer(1:12, 1:12, "-"))
A_P <- A_H
dimnames(A_H) <- list(host, host)
dimnames(A_P) <- list(partner, partner)
## W records which hosts interact with which partners: host i interacts
## with partners i - 1, i and i + 1.
W <- 1 * (abs(outer(1:12, 1:12, "-")) <= 1)
dimnames(W) <- list(host, partner)
## One kernel over all 24 species; the fixed rho sets how strongly
## interacting hosts and partners are linked.
K <- make_cross_kernel(A_H, A_P, W, rho = 0.9)

## Simulate one latent axis shared by all four traits. Host traits are
## measured on hosts only and partner traits on partners only: four
## replicate measurements per species, each with noise (SD 0.3).
lambda <- c(h_size = 0.9, h_defence = 0.6, p_size = 0.8, p_attack = -0.5)
z <- drop(t(chol(K)) %*% rnorm(24))
names(z) <- rownames(K)
dat <- expand.grid(rep = 1:4, trait = names(lambda),
                   species = rownames(K), stringsAsFactors = FALSE)
dat <- dat[startsWith(dat$species, "H") == startsWith(dat$trait, "h_"), ]
dat$value <- lambda[dat$trait] * z[dat$species] +
  rnorm(nrow(dat), sd = 0.3)
dat$trait <- factor(dat$trait, levels = names(lambda))
dat$species <- factor(dat$species, levels = rownames(K))
## Each record (one species, one replicate) is a sampling unit
dat$obs <- factor(paste(dat$species, dat$rep))

## One latent axis whose species scores are correlated as in K. The
## term is named "cross"; the functions below refer to it by name.
fit <- gllvmTMB(
  value ~ 0 + trait + kernel_latent(species, K = K, d = 1, name = "cross"),
  data = dat, trait = "trait", unit = "obs", cluster = "species"
)
## Coupled host and partner trait axes (point estimates only). With
## one latent axis (d = 1) there is a single module (n_modules = 1).
## This fit has no Psi, so the total covariance equals the shared part
## and the singular value is 1 by construction: read the trait
## loadings. Hosts with larger size and defence are linked with
## partners of larger size and lower attack. Flipping the signs of
## both axes gives the same module. A tier fitted with
## kernel_latent(..., unique = TRUE) has a Psi component; when Psi is
## estimated above zero the singular values fall below 1.
mods <- extract_coevolution_modules(
  fit,
  level = "cross",
  row_traits = c("h_size", "h_defence"),
  col_traits = c("p_size", "p_attack"),
  n_modules = 1
)
mods$modules
#>   component   module singular_value squared_share
#> 1     cross module_1              1             1
mods$row_axes
#>   component side   module     trait    loading
#> 1     cross  row module_1    h_size -0.8315610
#> 2     cross  row module_1 h_defence -0.5554335
mods$col_axes
#>   component   side   module    trait    loading
#> 1     cross column module_1   p_size -0.8799046
#> 2     cross column module_1 p_attack  0.4751505