Skip to contents

Bootstrap covariance, correlation, communality, and ICC summaries

Source: R/bootstrap-sigma.R
bootstrap_Sigma.Rd

Use bootstrap_Sigma() when Hessian, Wald, or profile intervals are unavailable or unsafe but the fitted model still has useful point estimates. A pdHess = FALSE or skipped sdreport() is an inference warning, not automatic proof that the fitted mean or rotation-invariant covariance summaries are unusable. Inspect check_gllvmTMB() / gllvmTMB_diagnose() first, then use this helper for parametric simulate-refit uncertainty.

Usage

bootstrap_Sigma(
  fit,
  n_boot = 200,
  level = c("unit", "unit_obs", "phy", "B", "W"),
  what = c("Sigma", "R", "communality", "ICC", "cross_corr"),
  conf = 0.95,
  seed = NULL,
  n_cores = 1,
  progress = TRUE,
  keep_draws = FALSE,
  link_residual = c("auto", "none")
)

Arguments

fit

A fit returned by gllvmTMB().

n_boot

Integer; number of bootstrap replicates. Default 200.

level

Character vector; which tier(s) to bootstrap. Use the canonical levels c("unit", "unit_obs", "phy"); legacy aliases "B" and "W" are still accepted. Levels absent from the fit are silently dropped. Default: all available levels.

what

Character vector; which summaries to compute. Subset of c("Sigma", "R", "communality", "ICC", "cross_corr"). Default: all. "ICC" only makes sense at the site level and requires both B and W tiers in the fit. "cross_corr" bootstraps the aggregate multiple_r between a multinomial() trait and each partner (see extract_cross_correlations()); it is stored as a plain named numeric per tier (multiple_r_B, ...) and is silently skipped for fits without a nominal trait.

conf

Numeric in (0, 1); confidence level for percentile CIs. Default 0.95.

seed

Optional RNG seed for reproducibility.

n_cores

Integer; number of cores for parallel refits. Default 1 (sequential). >= 2 uses future::multisession.

progress

Logical; print a one-line status message at each replicate (sequential only). Default TRUE.

keep_draws

Logical; if TRUE, the full n_boot x ... matrices of bootstrap draws are returned as $draws. Default FALSE (CIs only — saves memory for large n_boot).

How to treat family-specific link-implicit residual variance when extracting Sigma, R, communality, and ICC. "auto" (default) matches extract_Sigma() and adds the family/link residual to non-Gaussian trait diagonals; "none" returns the fitted model covariance without link residuals. Use "none" when validating against a DGP target defined as \(\Lambda\Lambda^\top + \Psi\).

Value

A list with components:

point_est

Named list of point estimates for each requested summary at each requested level (e.g. Sigma_B, R_B, communality_B, ICC_site).

ci_lower, ci_upper

Named lists of percentile CI bounds, element-wise the same shape as the corresponding point_est.

n_effective

Named list (vector summaries only) giving the per-element count of finite bootstrap draws. For multiple_r_* entries, any element whose effective count is below 80% of n_boot has its CI bounds set to NA (minimum-effective-B floor).

boot_median

Named list (multiple_r_* entries only) of the per-element bootstrap median – a cheap sanity comparator against point_est; a large gap flags gross bootstrap corruption and is not a coverage certifier.

ci_method

Character; currently "percentile".

link_residual

Character; the link-residual convention used in point estimates and bootstrap refit summaries.

conf, n_boot, n_failed

Configuration metadata.

draws

NULL unless keep_draws = TRUE; otherwise a named list of bootstrap draw arrays.

Details

The function generates n_boot parametric bootstrap replicates of a fitted model and returns percentile confidence intervals for the canonical biological summaries: trait covariance matrices \(\hat\Sigma_\mathrm{unit}\), \(\hat\Sigma_\mathrm{unit\_obs}\); the corresponding correlation matrices; per-trait communalities \(c_t^2 = (\Lambda \Lambda^\top)_{tt} / \Sigma_{tt}\); and per-trait site-level ICCs \(R_t = (\Sigma_\mathrm{unit})_{tt} / [(\Sigma_\mathrm{unit})_{tt} + (\Sigma_\mathrm{unit\_obs})_{tt}]\).

Scope: Gaussian bootstrap summaries and mixed-family refit plumbing are covered by current tests. Non-Gaussian bootstrap calibration is experimental: its repeated-sampling behaviour has not been certified, and there is no production calibration evidence for mixed-family intervals. Treat those intervals as indicative.

Each bootstrap replicate (1) draws a new response vector from simulate(fit, nsim = 1), (2) refits the model with the same formula on the simulated data, (3) extracts the requested summaries via extract_Sigma(), extract_communality(), and extract_ICC_site(). Replicates whose refit fails to converge are recorded but excluded from CI calculation.

Multicore is dispatched via future + future.apply; pass n_cores >= 2 to enable parallel refits. When parallel, replicates use future.apply's L'Ecuyer-CMRG seed stream so the answers are reproducible given a fixed seed, but they are NOT bit-identical to an n_cores = 1 run with the same seed (different RNG streams).

Caveats

  • Uses the existing simulate.gllvmTMB_multi() method. By default that method REDRAWS the random effects from the fitted covariance, which is what makes these intervals span between-unit variability. CIs reflect parametric simulate-refit variability, not a Bayesian posterior distribution for variance components.

  • Intervals are too narrow when the simulator cannot redraw a tier. Redraw is not implemented for every random-effect tier — notably the SPDE spatial tier and the diagonal phylogenetic tier. For a fit using one of those, simulate.gllvmTMB_multi() falls back to reusing the fitted random-effect modes and emits a one-shot warning. That fallback understates between-unit variability, so the intervals returned here are not calibrated for such fits. Treat the warning as a signal that these intervals cannot be trusted. Unsupported families fall back through the simulator's own warning path.

  • Refits use the same formula reconstructed from fit$formula and fit$covstructs, and forward the fit's auxiliary structure (phylo_vcv, phylo_tree, mesh, lambda_constraint) so each refit matches the original model.

  • Convergence: replicates whose refit fails or whose optimiser does not return convergence == 0 are counted in n_failed and excluded from CIs.

Examples

if (FALSE) { # \dontrun{
set.seed(1)
s <- simulate_site_trait(n_sites = 30, n_species = 1, n_traits = 3,
                         mean_species_per_site = 1,
                         Lambda_B = matrix(c(1, .5, -.4), 3, 1),
                         psi_B = c(.2, .15, .1))
fit <- gllvmTMB(value ~ 0 + trait + latent(0 + trait | site, d = 1),
                data  = s$data,
                trait = "trait",
                unit  = "site")
boot <- bootstrap_Sigma(fit, n_boot = 50, level = "unit",
                        what = c("Sigma", "R"), seed = 42)
boot$point_est$Sigma_B
boot$ci_lower$Sigma_B
boot$ci_upper$Sigma_B
} # }