$$c_t^2 \;=\; \frac{(\boldsymbol\Lambda \boldsymbol\Lambda^{\!\top})_{tt}}{(\boldsymbol\Lambda \boldsymbol\Lambda^{\!\top})_{tt} + \psi_{tt}}.$$
Arguments
- fit
A fit returned by
gllvmTMB(). Abootstrap_Sigma()result is also accepted when it containscommunalitysummaries; in that case the function reuses the stored point estimates and percentile bounds rather than refitting.- level
"unit"(between-unit),"unit_obs"(within-unit), or"phy"(phylogenetic tier). Legacy aliases"B"and"W"are accepted with a deprecation warning.- link_residual
For non-Gaussian fits:
"auto"(default) adds each trait's family/link implicit residual to the denominator;"none"returns communalities on the fitted model covariance scale without link-residual additions.- ci
Logical. When
TRUE, returns a tidy data frame with confidence-interval columns; whenFALSE(the default), returns a plain named numeric vector for backward compatibility.- conf_level
Confidence level when
ci = TRUE. Default 0.95.- method
One of
"wald"(default),"bootstrap", or"profile". Only used whenci = TRUE. The nonlinear penalty-profile prototype is currently withheld from the public release; requesting"profile"stops with an explanation. Bootstrap is a slower simulation-and-refit option and is not automatically coverage-calibrated.- nsim
Number of bootstrap replicates when
method = "bootstrap". Default 500.- seed
Optional RNG seed for the bootstrap.
Value
When ci = FALSE: a numeric vector indexed by trait.
When ci = TRUE: a data frame with columns trait, tier, c2,
lower, upper, method. For a bootstrap_Sigma input, the interval
columns are copied from the bootstrap object.
Details
The proportion of trait \(t\)'s variance that is shared with the
other traits via the latent factors. Bounded between 0 and 1. Calls
extract_Sigma() internally for the chosen level, so the diagonal
uses the full \(\boldsymbol\Sigma = \boldsymbol\Lambda \boldsymbol\Lambda^{\!\top} + \boldsymbol\Psi\)
decomposition. Ordinary latent() includes \(\boldsymbol\Psi\)
by default.
Caveat: communality with no-Psi fits
If the fit uses latent(..., unique = FALSE) at the requested
level, then \(\boldsymbol\Psi = \mathbf 0\) and c_t^2 = 1 for
every trait. This is mathematically correct for the no-residual
subset but tells you nothing about trait integration. The
extract_Sigma() advisory message will fire to flag this. To get
meaningful communalities, use ordinary latent() with the default
unique = TRUE.
For non-Gaussian fits the per-trait link-implicit residual (e.g.
\(\pi^2/3\) for binomial logit, \(\log(1 + 1/\hat\mu)\) for Poisson
log) is added to the denominator by default when link_residual = "auto";
pass link_residual = "none" to suppress (see extract_Sigma()).
Examples
## Simulate 6 traits recorded on a few species at each of 80 sites.
## Between sites, the traits share one latent axis (the loadings in
## Lambda_B) and each trait also has its own unique variance (psi_B).
## sigma2_eps is the variance of the remaining noise in each
## observation. beta = 0 turns off the simulator's environmental
## effects, which this model does not include.
set.seed(2)
Lambda_B <- matrix(c(0.9, 0.7, 0.5, -0.6, 0.4, 0.8), ncol = 1)
sim <- simulate_site_trait(
n_sites = 80, n_species = 10, n_traits = 6,
mean_species_per_site = 3,
beta = matrix(0, 6, 2),
Lambda_B = Lambda_B,
psi_B = rep(0.3, 6),
sigma2_eps = 0.3
)
## One latent axis between sites. unique = TRUE also estimates each
## trait's unique variance; without it every communality would be 1.
fit <- gllvmTMB(
value ~ 0 + trait +
latent(0 + trait | site, d = 1, unique = TRUE),
data = sim$data,
trait = "trait",
unit = "site"
)
## Communality of each trait between sites: the share of its
## between-site variance that comes from the shared latent axis
extract_communality(fit, level = "unit")
#> trait_1 trait_2 trait_3 trait_4 trait_5 trait_6
#> 0.7643490 0.5091306 0.4373497 0.5543386 0.3083301 0.6185751
## The true values used to simulate the data, for comparison
shared <- rowSums(Lambda_B^2)
round(shared / (shared + 0.3), 2)
#> [1] 0.73 0.62 0.45 0.55 0.35 0.68
## Wald intervals, built on the logit scale so they stay between 0
## and 1. They are approximate and are not guaranteed to reach their
## nominal coverage, especially with few units.
extract_communality(fit, level = "unit", ci = TRUE, method = "wald")
#> trait tier c2 lower upper method
#> 1 trait_1 unit 0.7643490 0.5509782 0.8955493 wald
#> 2 trait_2 unit 0.5091306 0.2968136 0.7182040 wald
#> 3 trait_3 unit 0.4373497 0.2393765 0.6575174 wald
#> 4 trait_4 unit 0.5543386 0.3437614 0.7470629 wald
#> 5 trait_5 unit 0.3083301 0.1273185 0.5766419 wald
#> 6 trait_6 unit 0.6185751 0.4117122 0.7898300 wald
# \donttest{
## Bootstrap intervals: bootstrap_Sigma() refits the model to data
## simulated from the fit. A small n_boot keeps the run time down and
## makes bootstrap_Sigma() warn that so few resamples give noisy
## bounds; real analyses need several hundred (the default is 999).
## These percentile intervals are not guaranteed to reach their
## nominal coverage either.
boot <- bootstrap_Sigma(fit, n_boot = 50, level = "unit",
what = "communality", progress = FALSE)
#> Warning: `n_boot` = 50 gives noisy percentile bounds.
#> ℹ Percentile CIs are first-order accurate, so more draws buy precision, not
#> correctness -- but below ~999 the Monte Carlo error in the endpoints is
#> itself material.
#> ℹ `n_cores` >= 2 parallelises the refits; the cost here is wall-clock, not
#> accuracy.
#> ℹ For Sigma diagonals, `profile_ci_total_variance()` needs no resampling at
#> all.
extract_communality(boot, level = "unit", ci = TRUE)
#> trait tier c2 lower upper method
#> 1 trait_1 unit 0.7643490 0.5343264 0.9081314 bootstrap
#> 2 trait_2 unit 0.5091306 0.2088676 0.7680402 bootstrap
#> 3 trait_3 unit 0.4373497 0.2814722 0.6924666 bootstrap
#> 4 trait_4 unit 0.5543386 0.3643615 0.7706691 bootstrap
#> 5 trait_5 unit 0.3083301 0.1553823 0.4751865 bootstrap
#> 6 trait_6 unit 0.6185751 0.4365552 0.7717693 bootstrap
# }