Confidence intervals for a fitted gllvmTMB model

Returns 95% (or other-level) confidence intervals for fixed effects, variance components, and trait covariance matrices, with three method choices:

Usage

# S3 method for class 'gllvmTMB_multi'
confint(
  object,
  parm,
  level = 0.95,
  method = c("profile", "wald", "bootstrap"),
  nsim = 500L,
  seed = NULL,
  ...
)

Arguments

object

A fit returned by gllvmTMB().

parm

One of:

"Sigma_unit" – between-unit trait covariance matrix.
"Sigma_unit_obs" – within-unit trait covariance matrix.
"Sigma_cluster" and "Sigma_cluster2" – diagonal covariance matrices for the first and second extra grouping tiers. "profile" and "wald" return diagonal variance intervals; "bootstrap" is intentionally blocked until a simulate-refit calibration gate is added for these tiers.
"sigma_phy" – per-trait phylogenetic standard deviations.
Legacy aliases "Sigma_B" and "Sigma_W", retained for existing scripts.
"Lambda" (all free entries of Lambda_unit), "Lambda:i,j" (single entry), or "Lambda:i,j;k,l" (multiple entries, semicolon-separated). Routes to loading_ci() / loading_profile(). For these tokens the method choices are c("wald", "wald_asym", "profile") with default "wald". Wald and profile rows target raw Lambda; "wald_asym" targets the standardised loading rho and labels rows accordingly.
"icc" (all traits), "icc:<trait_name>" (one trait by name), "icc:<t1>;<t2>" (multiple by name), or "icc:[1,3]" (1-based trait indices). Routes to extract_repeatability() with "wald" (default) or "bootstrap". The former profile token is withdrawn because it estimated a diagonal-only ratio rather than canonical full-covariance repeatability.
"phylo_signal" / "phylo_signal:<trait>" etc. – same grammar as "icc". Routes to profile_ci_phylo_signal() for "profile" and the companion Wald/bootstrap helpers for "wald" / "bootstrap". For multi-component fits, "profile" falls back to numerical Wald bounds with an explicit method label until the full fix-and-refit profile is implemented.
"communality:<tier>" (one tier, all traits) or "communality:<tier>:<trait>" (one tier, one trait). Tier is one of "unit" / "unit_obs" / "phy" (legacy "B" / "W"). Defaults to Wald; bootstrap is available explicitly. The nonlinear profile prototype is withheld.
"rho:<tier>:i,j" (one pair) or "rho:<tier>:i,j;k,l" (multiple pairs). Tier is one of "unit" / "unit_slope" / "unit_obs" / "cluster" / "cluster2" / "phy" / "spatial" (legacy "B" / "W" / "spde"). "unit_slope" indexes the augmented 2T coefficient vector by numeric position and is not currently an interval target. "cluster" and "cluster2" are recognised but profile intervals are withheld (see the error from extract_correlations()). Other tiers route to extract_correlations() ("fisher-z" / "wald" / "bootstrap"); Fisher-z is the default for this confint() token. The nonlinear profile prototype is withheld.
"proportion" (all components, all traits), "proportion:<component>" (one component, all traits), "proportion:<component>:<trait>" (one (component, trait)), "proportion:<component>:<t1>;<t2>" (one component, multiple traits), or "proportion:<c1>;<c2>" (multiple components). Components are by name ("shared_unit", "unique_unit", "shared_unit_obs", "unique_unit_obs", "unique_cluster", "unique_cluster2", "shared_phy", "unique_phy", "link_residual"). Defaults to Wald; bootstrap is available explicitly. The nonlinear profile prototype is withheld.
An integer index vector or character vector of fixed-effect term names (same as the standard confint() interface).
Missing (default) – all fixed-effect parameters.

level

Confidence level in (0, 1). Default 0.95.

method

One of "profile" (default), "wald", "bootstrap". For parm = "Lambda..." the accepted methods are c("wald", "wald_asym", "profile") with default "wald" (matching the base R confint() convention). "wald_asym" returns Fisher-z intervals for standardised loadings rather than raw Lambda.

nsim

Number of bootstrap replicates passed to bootstrap_Sigma() when method = "bootstrap". Default 500. Use a small value (e.g. 50) during development or testing.

seed

Optional integer RNG seed forwarded to bootstrap_Sigma() (only meaningful when method = "bootstrap").

...

Additional arguments currently unused.

Value

Sigma path – a data.frame with columns parameter (character, e.g. "Sigma_unit[t1,t1]"), estimate (point estimate), lower, upper, and method (the method used for that row, or "structural_zero" for off-diagonal entries that are fixed at zero by a diagonal-only tier). The parameter prefix follows the requested parm, so legacy calls still return legacy "Sigma_B[...]" or "Sigma_W[...]" labels.
Lambda path – a data.frame with columns parameter (e.g. "Lambda[trait_1,LV1]" for raw rows or "rho[trait_1,LV1]" for standardised rows), estimate, lower, upper, method, loading_scale, pd_hessian, and ci_status. When pdHess = FALSE on a Wald path, lower/upper are NA and ci_status flags the reason).
Derived-quantity path ("icc" / "phylo_signal" / "communality" / "rho" / "proportion") – a numeric matrix with two columns named after the requested level (e.g. "2.5 %" / "97.5 %") and rownames identifying the entry, e.g. "icc:trait_1", "communality:unit:trait_1", "rho:unit:1,2", or "proportion:shared_unit:trait_1".
Fixed-effects / variance-component path – a numeric matrix with rownames = parameter names and two columns named "2.5 %" / "97.5 %" (or the analogous quantiles for the requested level).

Details

method = "profile" (default): profile-likelihood bounds via TMB::tmbprofile() + stats::uniroot(). A successful numerical route can represent asymmetric likelihood shape, but its reference distribution and coverage remain target-specific.
method = "wald": Gaussian-approximation CIs from sd_report. Fastest; poor near boundaries.
method = "bootstrap": fitted-model parametric bootstrap via bootstrap_Sigma(). Slowest; inspect failed refits, Monte Carlo resolution, and whether simulation covers every active model component.

Fixed-effect and direct-parameter intervals use the requested method where supported, and Sigma-matrix intervals accept canonical parm = "Sigma_unit" / "Sigma_unit_obs" names, diagonal extra-grouping "Sigma_cluster" / "Sigma_cluster2" names, plus legacy "Sigma_B" / "Sigma_W" aliases. Profile intervals for full decomposed Sigma entries fall back to bootstrap because those entries are nonlinear functions of rotation-equivalent loadings and diagonal \(\Psi\). Non-Gaussian bootstrap and broader derived-target coverage are not universally calibrated; inspect the returned method and the target-specific article before reporting bounds. Fits using non-unit likelihood weights target a weighted estimating objective, not the ordinary likelihood used by these interval routes. confint() therefore refuses them; no sandwich interval is currently certified.

Main parm-class dispatch paths:

Sigma matrices – when parm is one of "Sigma_unit", "Sigma_unit_obs", "Sigma_cluster", "Sigma_cluster2", or "sigma_phy" (legacy aliases "Sigma_B" and "Sigma_W" still work), returns a tidy data.frame with columns parameter, estimate, lower, upper, method. Profile is computed element-wise via TMB::tmbprofile() for the diagonal entries; structural-zero off-diagonal entries in pure-diagonal tiers are labelled "structural_zero"; off-diagonals for reduced-rank tiers fall back to bootstrap (full Sigma sampling) since they mix two parameters in a non-linear way.
Fixed effects / variance components – when parm is missing, an integer index, or a character vector of fixed-effect term names, returns a numeric matrix with rows = parameters and columns = lower / upper bounds (same shape as stats::confint()). Method choice applies.
Lambda and derived summaries – "Lambda...", "icc...", "phylo_signal...", "communality...", "rho...", and "proportion..." tokens dispatch to their corresponding direct-profile, Wald, or bootstrap helper where that route is implemented. Nonlinear penalty-profile prototypes are withheld.

References

Pawitan, Y. (2001). In All Likelihood: Statistical Modelling and Inference Using Likelihood, Oxford University Press, ch. 9.

Venzon, D. J. & Moolgavkar, S. H. (1988). A method for computing profile-likelihood-based confidence intervals. Applied Statistics 37, 87-94. doi:10.2307/2347496

McCune, K. B., et al. (2024) coxme_icc_ci() – the Nakagawa-authored coxme-based profile-CI helper that inspired this work, in https://github.com/kelseybmccune/Time-to-Event_Repeatability/blob/main/R/rptRsurv.R.

Examples

set.seed(1)
# Four traits measured twice on each of 80 units (e.g. individuals).
# All four traits load on one shared latent axis, and each trait also
# has its own unit-level variance.
sim <- simulate_unit_trait(n_units = 80, n_obs_per_unit = 2,
                           n_traits = 4,
                           Lambda_B = matrix(c(0.9, 0.8, -0.7, 0.5)),
                           psi_B = 0.3, sigma2_eps = 0.3)
fit <- gllvmTMB(value ~ 0 + trait +
                  latent(0 + trait | unit, d = 1, unique = TRUE),
                data = sim$data, trait = "trait", unit = "unit")

# Wald intervals for the fixed effects (here the four trait means).
# Without `method`, confint() profiles the likelihood instead, which
# is slower: one profile per coefficient.
confint(fit, method = "wald")
#>                      2.5 %     97.5 %
#> traittrait_1 -0.8194989650 -0.3936216
#> traittrait_2 -0.0005766122  0.4573286
#> traittrait_3 -1.0772043198 -0.6825416
#> traittrait_4  1.4955066296  1.8201174

# Derived summaries are requested by name: the unit-level correlation
# between traits 1 and 3 (a Fisher-z interval), and each trait's
# communality, the share of its unit-level variance that comes from
# the shared axis (Wald intervals).
confint(fit, parm = "rho:unit:1,3")
#>                   2.5 %     97.5 %
#> rho:unit:1,3 -0.7891776 -0.5528873
confint(fit, parm = "communality:unit")
#>                               2.5 %    97.5 %
#> communality:unit:trait_1 0.45324851 0.9337573
#> communality:unit:trait_2 0.43949414 0.8873570
#> communality:unit:trait_3 0.35614548 0.8199124
#> communality:unit:trait_4 0.08083795 0.5489518

# Trait covariance matrices. With indep(), each diagonal entry of
# Sigma_unit is one variance parameter, so its profile-likelihood
# interval is computed directly; the off-diagonal entries are zero by
# construction.
sim2 <- simulate_unit_trait(n_units = 50, n_obs_per_unit = 2,
                            n_traits = 3, psi_B = c(1, 0.6, 0.4),
                            sigma2_eps = 0.3)
fit2 <- gllvmTMB(value ~ 0 + trait + indep(0 + trait | unit),
                 data = sim2$data, trait = "trait", unit = "unit")
confint(fit2, parm = "Sigma_unit")
#>                     parameter  estimate     lower     upper          method
#> 1 Sigma_unit[trait_1,trait_1] 1.1047205 0.7181747 1.7564188         profile
#> 2 Sigma_unit[trait_1,trait_2] 0.0000000 0.0000000 0.0000000 structural_zero
#> 3 Sigma_unit[trait_2,trait_2] 0.6075270 0.3727453 1.0006026         profile
#> 4 Sigma_unit[trait_1,trait_3] 0.0000000 0.0000000 0.0000000 structural_zero
#> 5 Sigma_unit[trait_2,trait_3] 0.0000000 0.0000000 0.0000000 structural_zero
#> 6 Sigma_unit[trait_3,trait_3] 0.4721189 0.2783097 0.7948474         profile

# With latent(), each Sigma_unit entry combines several parameters, so
# confint(fit, parm = "Sigma_unit") switches to a parametric bootstrap
# that refits the model `nsim` times. Use nsim = 999 or more (smaller
# values give a warning); that takes minutes, so it is not run here.