Wald confidence intervals on augmented random-slope standard deviations

Per-trait confidence intervals on the standard deviation of an augmented random-slope term. All intervals are delta-method Wald intervals built on the log-SD scale so the bounds stay positive by construction. Three underlying routes are covered, dispatched automatically from the fit's structure:

The ordinary augmented diagonal (Psi) route (theta_diag_B_slope, from latent(0 + trait + (0 + trait):x | unit, d = K) or unique(...)). theta_diag_B_slope is a genuine univariate log-SD (sd = exp(theta)), so the interval is read directly off sdreport()'s cov.fixed – no Jacobian.
The phylogenetic Cholesky augmented-slope route (theta_dep_chol, from phylo_dep(1 + x1 + ... + xs | species) / phylo_indep(1 + x | species)). The marginal slope SD is a nonlinear function of several packed Cholesky entries ( Var(slope_t) = L21_t^2 + exp(diag_slope_t)^2 for a single slope), so its standard error comes from an ADREPORT()ed quantity in the C++ template (sd_b), letting sdreport() run the multivariate delta method against the authoritative packing.
The loadings-only augmented random-slope route (theta_rr_B_slope with no diagonal companion, e.g. latent(..., unique = FALSE)). The marginal slope variance is a quadratic form in several loading entries (Lambda_B_slope \%*\% t(Lambda_B_slope)), again handled via an ADREPORT()ed quantity (sd_rr_B_slope).

Usage

slope_sd_ci(fit, level = 0.95, scale = c("sd", "variance"))

Arguments

fit: A multi-trait gllvmTMB() fit with at least one augmented random-slope structure covered above.
level: Confidence level. Defaults to 0.95.
scale: "sd" (default) returns the standard deviation scale; "variance" returns the variance scale (estimate^2, bounds exp(2 * (theta +/- z * se(theta)))). Both scales are computed from the identical theta / se_theta, so they agree by construction. total_sd is always reported on the SD scale regardless of scale; total_lower / total_upper follow scale like lower / upper.

Value

A data frame of class gllvmTMB_slope_ci, one row per slope SD (trait x slope term), with columns:

trait: Trait factor level.
term: The slope covariate name.
estimate, lower, upper: On the requested scale. lower / upper are NA whenever status != "ok".
component: "unique_psi" when the diagonal route's fit also carries a shared loadings block, so estimate is only part of the marginal slope SD; "total" otherwise, meaning estimate already is the total marginal slope SD.
total_sd: Point estimate of the TOTAL marginal slope SD, on the SD scale regardless of scale. Equal to estimate when component == "total".
total_lower, total_upper: Interval on total_sd, on the requested scale. NA whenever total_status != "ok" or the supporting ADREPORT is unavailable.
total_status: Kill-switch status for the total_sd interval; see below. "unavailable" when no supporting ADREPORT exists (only possible for the diagonal route on a stale cached fit).
theta, se_theta: The log-SD estimate and its standard error, for auditability. For the diagonal route these are the raw fixed-effect parameter and its cov.fixed SE (unchanged from slice 1); for the phylogenetic and loadings-only routes these are derived from an ADREPORTed natural-scale quantity via the delta method se(log(X)) = se(X) / X.
method: "wald_log_scale" for the diagonal route (direct parameter read), "wald_log_scale_adreport" for the phylogenetic and loadings-only routes (ADREPORT delta method).
interval_status: Always "wald_uncalibrated".
status: "ok", "no_pd_hessian", "se_nonfinite", "se_blowup", or "near_zero_relative" – see Kill-switch guard below.
scale: The requested scale, "sd" or "variance".

The returned object also carries a hard-coded attribute calibrated = FALSE (never an argument; only a future coverage certificate can flip this in the source), and the requested level.

Scope boundary – what this interval IS and IS NOT

Every interval returned by this function is a transformed Wald interval read off sdreport(), the same inferential instrument as loading_ci() and confint_inspect(): it uses the asymptotic normal approximation and propagates it through a monotone log/exp transform. It is not a calibrated coverage statement – no repeated-sampling coverage campaign has been run for any of these estimands, and interval_status = "wald_uncalibrated" marks every row as such. The implementation is tested, but repeated-sampling coverage is not certified.

When the fit's random-slope term carries BOTH a diagonal Psi companion and a shared loadings component (theta_rr_B_slope, the default alongside theta_diag_B_slope for ordinary latent()), estimate for that route is the per-trait UNIQUE (Psi) component of slope variance only – it excludes the shared loadings contribution to the marginal slope variance, and can understate the total marginal slope SD substantially (measured up to ~45% in a recovery study). This restriction travels with the data, not only the print method: the component column reads "unique_psi" in that case ("total" otherwise – including for the phylogenetic and loadings-only routes, where the reported quantity already IS the total marginal slope SD). total_sd always reports the total marginal slope SD as a point estimate. total_lower / total_upper report a genuine interval on that total wherever the C++ template ADREPORT()s the combined expression directly (both routes as of this version); total_status names the same kill-switch statuses as status when that support is present, or "unavailable" when it is not (e.g. a fit refit under an older package version whose ADREPORT does not yet exist). A cli::cli_warn() fires on the call itself whenever the shared loadings block is present alongside the diagonal route, not only when the result is printed.

Kill-switch guard – point estimates always, intervals never from a fit that cannot support one

Applied identically to every route. lower/upper (and, where supported, total_lower/total_upper) are set to NA (with status explaining why and a cli::cli_warn()) whenever any of:

fit$sd_report$pdHess is FALSE (status = "no_pd_hessian", applies to every row);
the log-SD estimate or its standard error is non-finite (status = "se_nonfinite") – for the ADREPORT-based routes this also covers a non-positive natural-scale estimate, which cannot be logged;
se_theta > 10 on the log scale – essentially information-free (an SE of 10 on the log scale means the interval would span $e^{40}$); status = "se_blowup";
sd_hat is at most 1% of the largest sd_hat among this fit's OTHER slope coordinates on the SAME route (only evaluated when there is more than one to compare against); status = "near_zero_relative", mirroring the package's relative-to-siblings convention for a collapsed variance component (psi_rel_thresh / near_zero_psi_*, R/diagnose.R).

The point estimate (estimate, theta) is always returned, following the house line that a non-PD Hessian disqualifies standard errors but not point estimates (R/bootstrap-sigma.R, R/cv-internal.R).

Examples

if (FALSE) { # \dontrun{
set.seed(1)
n_ind <- 30L
n_traits <- 2L
trait_levels <- paste0("t", seq_len(n_traits))
individuals <- paste0("id", seq_len(n_ind))
df <- expand.grid(
  individual = factor(individuals, levels = individuals),
  rep = 1:5,
  trait = factor(trait_levels, levels = trait_levels)
)
df$x <- stats::rnorm(nrow(df))
alpha <- c(0.2, -0.1)
beta <- c(0.3, -0.2)
psi_noise <- matrix(stats::rnorm(4L * n_ind, sd = 0.25), nrow = 4L)
eta <- numeric(nrow(df))
for (o in seq_len(nrow(df))) {
  tt <- as.integer(df$trait[o])
  ii <- as.integer(df$individual[o])
  base <- 2L * (tt - 1L)
  eta[o] <- alpha[tt] + beta[tt] * df$x[o] +
    psi_noise[base + 1L, ii] + psi_noise[base + 2L, ii] * df$x[o]
}
df$value <- eta + stats::rnorm(nrow(df), sd = 0.3)

options(gllvmTMB.quiet_grammar_notes = TRUE)
fit <- gllvmTMB(
  value ~ 0 + trait + (0 + trait):x +
    unique(0 + trait + (0 + trait):x | individual),
  data = df,
  trait = "trait",
  unit  = "individual",
  control = gllvmTMBcontrol(se = TRUE, optimizer = "optim",
                             optArgs = list(method = "BFGS"))
)
slope_sd_ci(fit)
slope_sd_ci(fit, level = 0.90, scale = "variance")
} # }