Predict pair-specific cross-lineage covariance
Source:R/extract-sigma.R
predict_cross_covariance.Rdpredict_cross_covariance() combines a component-specific
Gamma_shape block from extract_Gamma() with entries of the fitted
dense kernel matrix. For a fixed cross-lineage kernel built with
make_cross_kernel(), the off-diagonal K entries already include the
supplied fixed rho, so the pair-specific covariance is
Gamma_shape * K[row_level, col_level]. This helper therefore uses the
shape-scale Gamma and does not multiply by extract_Gamma(scale = "effect").
Fixed dense kernel_latent() tiers, covered by the package's validation
tests, store their
aligned K matrices on the fit, and this helper returns point estimates for
named species/lineage pairs. It does not estimate rho, produce
intervals, calibrate null thresholds, or combine components into a universal
total Gamma.
Usage
predict_cross_covariance(
fit,
level,
row_levels = NULL,
col_levels = NULL,
row_traits,
col_traits
)Arguments
- fit
A fitted
gllvmTMB_multiobject.- level
Character scalar naming the fitted
kernel_*()tier.- row_levels, col_levels
Character vectors naming row and column levels in the fitted kernel matrix. If both are omitted for a kernel built by
make_cross_kernel(), host and partner levels from the kernel metadata are used.- row_traits, col_traits
Character vectors of trait names defining the rows and columns of the component-specific
Gamma_shapeblock.
Value
A data frame with one row per level-pair and trait-pair. Columns
include component, row_level, col_level, row_trait, col_trait,
kernel_value, gamma_shape, covariance, rho, and
kernel_includes_rho.
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"
)
## Predicted covariance between host traits in hosts H1 and H2 and
## partner traits in partners P1 and P12 (point estimates only). Each
## covariance is the kernel entry for that host-partner pair times
## the Gamma entry for that trait pair (kernel_value * gamma_shape).
## Both hosts interact with P1 but not with the distant P12, so their
## predicted covariances with P12 are close to zero.
predict_cross_covariance(
fit,
level = "cross",
row_levels = c("H1", "H2"),
col_levels = c("P1", "P12"),
row_traits = c("h_size", "h_defence"),
col_traits = c("p_size", "p_attack")
)
#> component row_level col_level row_trait col_trait kernel_value gamma_shape
#> 1 cross H1 P1 h_size p_size 0.510411409 0.7668029
#> 2 cross H2 P1 h_size p_size 0.548041718 0.7668029
#> 3 cross H1 P12 h_size p_size 0.003408946 0.7668029
#> 4 cross H2 P12 h_size p_size 0.005821662 0.7668029
#> 5 cross H1 P1 h_defence p_size 0.510411409 0.5121792
#> 6 cross H2 P1 h_defence p_size 0.548041718 0.5121792
#> 7 cross H1 P12 h_defence p_size 0.003408946 0.5121792
#> 8 cross H2 P12 h_defence p_size 0.005821662 0.5121792
#> 9 cross H1 P1 h_size p_attack 0.510411409 -0.4140754
#> 10 cross H2 P1 h_size p_attack 0.548041718 -0.4140754
#> 11 cross H1 P12 h_size p_attack 0.003408946 -0.4140754
#> 12 cross H2 P12 h_size p_attack 0.005821662 -0.4140754
#> 13 cross H1 P1 h_defence p_attack 0.510411409 -0.2765780
#> 14 cross H2 P1 h_defence p_attack 0.548041718 -0.2765780
#> 15 cross H1 P12 h_defence p_attack 0.003408946 -0.2765780
#> 16 cross H2 P12 h_defence p_attack 0.005821662 -0.2765780
#> covariance rho kernel_includes_rho
#> 1 0.3913849703 0.9 TRUE
#> 2 0.4202400017 0.9 TRUE
#> 3 0.0026139900 0.9 TRUE
#> 4 0.0044640673 0.9 TRUE
#> 5 0.2614220907 0.9 TRUE
#> 6 0.2806955508 0.9 TRUE
#> 7 0.0017459912 0.9 TRUE
#> 8 0.0029817338 0.9 TRUE
#> 9 -0.2113488109 0.9 TRUE
#> 10 -0.2269305962 0.9 TRUE
#> 11 -0.0014115608 0.9 TRUE
#> 12 -0.0024106069 0.9 TRUE
#> 13 -0.1411685481 0.9 TRUE
#> 14 -0.1515762622 0.9 TRUE
#> 15 -0.0009428394 0.9 TRUE
#> 16 -0.0016101433 0.9 TRUE