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" and scale = "sd".
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 calculation is tested; its repeated-sampling coverage is not.

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 simulation). 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_*).

The point estimate (estimate, theta) is always returned. A non-positive-definite Hessian disqualifies standard errors, but not point estimates.

Examples

set.seed(1)
# Two traits (t1, t2) measured on five occasions for each of 60
# individuals. An environmental covariate x is recorded once per
# occasion and applies to both traits.
n_ind <- 60
dat <- expand.grid(occasion = 1:5, individual = factor(1:n_ind),
                   trait = factor(c("t1", "t2")))
dat$x <- rep(rnorm(5 * n_ind), times = 2)

# Each individual has its own intercept and its own slope on x for each
# trait, in the order (t1 intercept, t1 slope, t2 intercept, t2 slope).
# Two latent axes generate them: the first shifts both intercepts, the
# second shifts both slopes. The true slope SDs are 0.5 and 0.4.
Lambda <- cbind(c(0.8, 0, 0.6, 0), c(0, 0.5, 0, 0.4))
coefs <- matrix(rnorm(2 * n_ind), n_ind) %*% t(Lambda)
ind <- as.integer(dat$individual)
tr <- as.integer(dat$trait)
dat$value <- c(0.2, -0.1)[tr] + c(0.3, -0.2)[tr] * dat$x +
  coefs[cbind(ind, 2 * tr - 1)] + coefs[cbind(ind, 2 * tr)] * dat$x +
  rnorm(nrow(dat), sd = 0.3)

# Random intercepts and slopes for both traits. Their covariance is
# built from two latent axes; unique = FALSE adds no extra per-trait
# variance on top of the axes.
fit <- gllvmTMB(
  value ~ 0 + trait + (0 + trait):x +
    latent(0 + trait + (0 + trait):x | individual, d = 2,
           unique = FALSE),
  data = dat, trait = "trait", unit = "individual"
)

# 95% Wald intervals for how much individuals differ in their slopes,
# one row per trait, on the standard-deviation scale.
ci <- slope_sd_ci(fit)
ci[, c("trait", "term", "estimate", "lower", "upper")]
#> Wald (log-SD-scale) intervals on augmented random-slope standard
#> deviations. These are recovery-only, UNCALIBRATED intervals --
#> repeated-sampling coverage has not been measured for this estimand,
#> for any family. Consult the package's current limitations for
#> the applicable validation boundary.
#> 
#>   trait term  estimate     lower     upper
#> 1    t1    x 0.5029887 0.4106722 0.6160573
#> 2    t2    x 0.3871676 0.3129357 0.4790081

# The same quantity as 90% intervals on the variance scale.
ci_var <- slope_sd_ci(fit, level = 0.90, scale = "variance")
ci_var[, c("trait", "term", "estimate", "lower", "upper")]
#> Wald (log-SD-scale) intervals on augmented random-slope standard
#> deviations. These are recovery-only, UNCALIBRATED intervals --
#> repeated-sampling coverage has not been measured for this estimand,
#> for any family. Consult the package's current limitations for
#> the applicable validation boundary.
#> 
#>   trait term  estimate     lower     upper
#> 1    t1    x 0.2529976 0.1800143 0.3555706
#> 2    t2    x 0.1498987 0.1048662 0.2142696