Bootstrap covariance, correlation, communality, and ICC summaries

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 = 999,
  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 999. Values below 2 / (1 - conf) - 1 (39 at the default conf = 0.95) trigger a warning: the widest interval B draws can produce is [min, max], whose coverage is at most (B - 1) / (B + 1), so a smaller B cannot reach the requested level whatever the data are. The fit still runs. Values below 999 warn for endpoint noise.

Two different things bound n_boot, and they bite at different scales. Below the arithmetic floor the arithmetic ceiling binds: the interval cannot reach the requested level however much data you have. Above it, what binds is Monte Carlo error in the endpoints, which is a precision question and shrinks with B. Conflating the two is how a replicate count comes to be treated as a free budget lever when the low end of it is a correctness constraint.

The default is 999 rather than 1000 so that (1 - conf) / 2 * (B + 1) is a whole number at conf = 0.95 (25), letting the percentile bounds land on order statistics rather than being interpolated between them.

Raising n_boot reduces Monte Carlo error in the endpoints but does not make a percentile interval second-order accurate — its coverage asymptotes slightly below nominal regardless. Use profile_ci_total_variance() for the Sigma diagonals if that matters. For exploratory work a smaller B (200, say) is a reasonable time-for-precision trade; for a number you intend to publish, do not go below the default.

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).

link_residual

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(). Native Laplace/AGHQ refits retain the fitted ridge, integration, optimiser and REML settings; only the replicate standard-error calculation is disabled. Older fits without the full control record recover retained ridge and integration metadata, but their unrecorded optimiser options cannot be recovered. 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).

Automatic refits of structured trait-intercept rho fits are not supported. Use point covariance and source-strength extraction; no rho interval is provided.

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
} # }