Per-trait repeatability with confidence intervals

Returns the per-trait repeatability $$R_t = v_{B,t} / (v_{B,t} + v_{W,t}),$$ where $$v_{B,t} = [\Lambda_B\Lambda_B^\top]_{tt} + \sigma^2_{B,t}$$ and $$v_{W,t} = [\Lambda_W\Lambda_W^\top]_{tt} + \sigma^2_{W,t} + \sigma^2_{d,t}.$$ The first two terms include the shared latent and diagonal companion variance at the unit and observation tiers; $\sigma^2_{d,t}$ is the family-specific link residual used for non-Gaussian traits. The function is intended for a fit returned by gllvmTMB() with the relevant ordinary unit and observation-level components. Also known as the intraclass correlation coefficient (ICC) at the unit level.

Usage

extract_repeatability(
  fit,
  level = 0.95,
  method = c("wald", "profile", "bootstrap"),
  nsim = 500L,
  seed = NULL
)

Arguments

fit: A fit returned by gllvmTMB. A bootstrap_Sigma object is also accepted when it contains an ICC_site summary; in that case the function reuses the stored point estimates and percentile bounds rather than refitting.
level: Confidence level. Default 0.95.
method: One of "wald" (default), "profile", or "bootstrap". "profile" is accepted only for backwards compatibility and aborts because a defensible profile interval for canonical full-covariance repeatability is not available.
nsim: Number of bootstrap replicates when method = "bootstrap". Default 500.
seed: Optional RNG seed for the bootstrap.

Value

A data frame with columns trait, R (point estimate), lower, upper, method.

Interval calibration

The point estimates are the supported claim. The Wald and bootstrap interval methods here are provided for exploration: their empirical coverage is not certified for this estimand, so treat intervals as exploratory rather than coverage-calibrated. See NEWS.md for the current coverage status.

Method choice

"wald" (default): Gaussian-approximation CI via the delta method on log(v_B) - log(v_W), transformed with plogis().
"profile": withdrawn. It aborts rather than silently substituting a different interval method or estimand.
"bootstrap": parametric bootstrap via bootstrap_Sigma().

References

Nakagawa, S. & Schielzeth, H. (2010) Repeatability for Gaussian and non-Gaussian data: a practical guide for biologists. Biological Reviews 85, 935-956. doi:10.1111/j.1469-185X.2010.00141.x

Examples

## Simulate 4 traits recorded on a few species at each of 50 sites.
## Every trait varies between sites (the "unit" tier) and between the
## species records within a site (the "unit_obs" tier). A zero beta
## matrix turns off the simulator's environmental effects, which this
## model does not include.
set.seed(3)
sim <- simulate_site_trait(
  n_sites = 50, n_species = 10, n_traits = 4,
  mean_species_per_site = 3,
  beta = matrix(0, 4, 2),
  Lambda_B = matrix(c(0.8, 0.6, -0.5, 0.4), 4, 1),
  psi_B = rep(0.3, 4),
  Lambda_W = matrix(c(0.5, -0.4, 0.4, 0.3), 4, 1),
  sigma2_eps = 0.3
)
## One latent factor plus trait-specific variances at each tier. Each
## species record has one value per trait, so the within-site
## trait-specific variances also absorb the measurement noise
## (gllvmTMB prints a note saying so).
fit <- gllvmTMB(
  value ~ 0 + trait +
    latent(0 + trait | site, d = 1, unique = TRUE) +
    latent(0 + trait | site_species, d = 1, unique = TRUE),
  data     = sim$data,
  trait    = "trait",
  unit     = "site",
  unit_obs = "site_species"
)
#> ℹ Auto-suppressing `sigma_eps`: `indep(0 + trait | site_species)` is at the
#>   per-row level, so it already absorbs the observation residual.
#> • Fixed at 0.00122 (~1/1000 of sd(y)) to keep the Gaussian density
#>   well-defined; the row-level residual variance is fully captured by the
#>   per-row diagonal term.
## Repeatability of each trait: the share of its variance that lies
## between sites, R = between-site / (between-site + within-site),
## where each variance adds up the shared (latent) and trait-specific
## parts of its tier. The Wald interval is built on the log ratio of
## the two variances, so it stays between 0 and 1; it is not
## guaranteed to cover the true value 95% of the time.
extract_repeatability(fit)
#>     trait         R     lower     upper method
#> 1 trait_1 0.7270923 0.6099918 0.8194405   wald
#> 2 trait_2 0.6563691 0.5307142 0.7633797   wald
#> 3 trait_3 0.6242875 0.4895909 0.7421582   wald
#> 4 trait_4 0.5815677 0.4371436 0.7132431   wald
## The repeatabilities used to simulate the data, for comparison.
## With only 50 sites the estimates above are noisy.
tr <- sim$truth
v_between <- tr$Lambda_B[, 1]^2 + tr$psi_B
v_within <- tr$Lambda_W[, 1]^2 + tr$sigma2_eps
round(v_between / (v_between + v_within), 2)
#> [1] 0.63 0.59 0.54 0.54
# \donttest{
## Bootstrap intervals: bootstrap_Sigma() refits the model to 50 data
## sets simulated from the fit, so this part takes tens of seconds.
## 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 = c("unit", "unit_obs"),
                        what = "ICC", 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.
## Read the stored bootstrap intervals (no refitting). To do both steps
## in one call, use method = "bootstrap" with nsim and seed instead:
## extract_repeatability(fit, method = "bootstrap", nsim = 50, seed = 1)
extract_repeatability(boot)
#>     trait         R     lower     upper    method
#> 1 trait_1 0.7270923 0.6132359 0.7798397 bootstrap
#> 2 trait_2 0.6563691 0.5275736 0.7706151 bootstrap
#> 3 trait_3 0.6242875 0.4972395 0.7362115 bootstrap
#> 4 trait_4 0.5815677 0.4229509 0.7210811 bootstrap
# }