Confidence intervals on individual entries of the loading matrix
Source:R/loading-ci.R
loading_ci.RdReturn per-entry CIs on either the raw reduced-rank loading matrix
Lambda_<level> or the standardised loading
rho[t,k] = Lambda[t,k] / sqrt(Sigma_total[t,t]) of a confirmatory
gllvmTMB() fit. A numerical Jacobian of the complete requested loading
vector is combined with the TMB sdreport() covariance. On the
standardised scale this includes covariance among axes and uncertainty in
fitted denominator components such as Psi.
Usage
loading_ci(
fit,
level = "unit",
method = c("wald", "wald_asym", "profile"),
conf_level = 0.95,
sigma_d2 = 1,
loading_scale = NULL
)Arguments
- fit
A multivariate
gllvmTMB()fit object.- level
Which loading matrix to summarise:
"unit"(default) or"unit_obs". Legacy aliases"B"/"W"are accepted with a one-shot deprecation warning.- method
CI method:
"wald"(default)Symmetric Wald via the joint delta method on the requested
loading_scale."wald_asym"Asymmetric Wald via the Fisher-z transformation on \(\rho_{tk} = \Lambda_{tk}/\sqrt{\Sigma_{total,tt}}\). This method is available only with
loading_scale = "standardized"; its bounds stay in(-1, 1). It captures bounded-support asymmetry but not higher-order log-likelihood curvature."profile"Profile-likelihood inversion through
loading_profile(). This refits across a grid for each free loading entry and can be used when Wald inference is blocked by a non-positive-definite Hessian. The current profile target is rawLambda, soloading_scale = "standardized"is refused.
- conf_level
Confidence level. Defaults to 0.95.
- sigma_d2
Deprecated compatibility argument; ignored. Standardised inference now derives each denominator from model-implied total variance via
extract_Sigma()instead of accepting one scalar residual variance.- loading_scale
NULL,"raw", or"standardized".NULLuses"raw"for"wald"/"profile"and"standardized"for"wald_asym". Estimates, SEs, bounds, and downstream null regions are always on this recorded scale.
Value
A data frame (one row per Lambda entry) with columns
trait, axis, estimate, se, lower, upper, method,
loading_scale, and pinned. pinned = TRUE records that the raw
loading was fixed by lambda_constraint. Its raw SE is zero, but its
standardised SE can be nonzero because total variance can still vary.
Details
Per-entry CIs are only well-defined for confirmatory fits — i.e.
fits supplied with a lambda_constraint that fixes enough entries
to pin the rotation. Exploratory fits leave Lambda identified only
up to a d x d orthogonal rotation, so the SE on any single
Lambda[i, k] depends on the rotation convention and is not a
biological quantity. This function therefore errors on exploratory
fits and points the user at confirmatory_lambda() /
suggest_lambda_constraint(), or at extract_communality() /
extract_Sigma() for rotation-invariant summaries.
Interval calibration
Standardized symmetric joint-delta Wald intervals have an exact empirical
certificate only for the structurally free strict-lower targets in native
Laplace, pinned, unrotated, lower-triangular ordinary Gaussian unit-tier
three-trait cells (n_units=150,d=2), (n_units=400,d=1), and
(n_units=400,d=2), at conf_level = 0.95. The (n_units=150,d=1) cell
failed its lower-band gate. This is one frozen DGP: trait intercepts
(-0.20, 0.10, 0.25), unique standard deviations
(0.70, 0.80, 0.90), and loading vector (0.80, 0.45, -0.35) for d=1,
with the second column (0, 0.70, 0.40) added for d=2. Coverage is
conditional on eligible fits (optimizer convergence, converged fit health,
an available sdreport(), and a positive-definite Hessian); availability
was 98.82%, 93.38%, and 96.18% in the three certified cells, respectively.
The certificate does not extend to another loading, unique-variance, or
intercept regime. Pinned diagnostic rows, Fisher-z Wald
(method = "wald_asym"), arbitrary confirmatory constraints, raw-loading
intervals, profiles, bootstraps, other rotations, and every neighbouring
family, tier, rank, sample size, or confidence level remain exploratory.
The worked example below is binomial-probit and is outside the certified
cells. See NEWS.md for the current coverage status.
See also
flag_unreliable_loadings() for a decision-aid summary;
confirmatory_lambda() to build a confirmatory constraint matrix;
extract_communality() for rotation-invariant alternatives.
Examples
set.seed(1)
# Presence (1) or absence (0) of eight species at 200 sites. Odd-numbered
# species respond to one latent environmental gradient and even-numbered
# species to another; sp7 responds only weakly.
n_sites <- 200
Lambda <- cbind(c(1.0, 0, 0.8, 0, 0.7, 0, 0.2, 0),
c(0, 1.0, 0, 0.9, 0, 0.7, 0, 0.6))
scores <- matrix(rnorm(2 * n_sites), n_sites, 2)
dat <- expand.grid(site = factor(seq_len(n_sites)),
trait = factor(paste0("sp", 1:8)))
liability <- rowSums(scores[dat$site, ] * Lambda[dat$trait, ])
dat$present <- rbinom(nrow(dat), 1, pnorm(liability))
# Confirmatory hypothesis: group A (odd species) loads only on axis 1 and
# group B (even species) only on axis 2. These fixed zeros pin down the
# axes (up to sign); anchors = NA leaves all other loadings free.
M <- confirmatory_lambda(species = levels(dat$trait),
group = rep(c("A", "B"), times = 4), d = 2,
loads_on = list(A = 1, B = 2),
anchors = c(NA, NA))
M
#> LV1 LV2
#> sp1 NA 0
#> sp2 0 NA
#> sp3 NA 0
#> sp4 0 NA
#> sp5 NA 0
#> sp6 0 NA
#> sp7 NA 0
#> sp8 0 NA
# One 0/1 record per species and site cannot separate extra
# species-specific variation from coin-flip noise, so the latent term
# carries loadings only (unique = FALSE).
fit <- gllvmTMB(present ~ 0 + trait +
latent(0 + trait | site, d = 2, unique = FALSE),
data = dat, family = binomial(link = "probit"),
trait = "trait", unit = "site",
lambda_constraint = list(unit = M))
# 95% Wald intervals for every loading, on the probit (liability) scale.
# Fixed entries (pinned = TRUE) have zero width; the interval for sp7's
# weak loading is wide relative to its size.
ci <- loading_ci(fit)
ci[, c("trait", "axis", "estimate", "se", "lower", "upper", "pinned")]
#> trait axis estimate se lower upper pinned
#> 1 sp1 LV1 0.9633789 0.2794047 0.41575581 1.5110021 FALSE
#> 2 sp2 LV1 0.0000000 0.0000000 0.00000000 0.0000000 TRUE
#> 3 sp3 LV1 0.5854798 0.1780835 0.23644247 0.9345171 FALSE
#> 4 sp4 LV1 0.0000000 0.0000000 0.00000000 0.0000000 TRUE
#> 5 sp5 LV1 0.7235048 0.2075922 0.31663160 1.1303779 FALSE
#> 6 sp6 LV1 0.0000000 0.0000000 0.00000000 0.0000000 TRUE
#> 7 sp7 LV1 0.2247761 0.1348475 -0.03952026 0.4890724 FALSE
#> 8 sp8 LV1 0.0000000 0.0000000 0.00000000 0.0000000 TRUE
#> 9 sp1 LV2 0.0000000 0.0000000 0.00000000 0.0000000 TRUE
#> 10 sp2 LV2 0.8214755 0.2261517 0.37822636 1.2647247 FALSE
#> 11 sp3 LV2 0.0000000 0.0000000 0.00000000 0.0000000 TRUE
#> 12 sp4 LV2 0.8054023 0.2265470 0.36137831 1.2494263 FALSE
#> 13 sp5 LV2 0.0000000 0.0000000 0.00000000 0.0000000 TRUE
#> 14 sp6 LV2 0.7417159 0.2046233 0.34066165 1.1427702 FALSE
#> 15 sp7 LV2 0.0000000 0.0000000 0.00000000 0.0000000 TRUE
#> 16 sp8 LV2 0.4217512 0.1499908 0.12777456 0.7157278 FALSE
# Standardized loadings divide each loading by the species' total
# latent-scale standard deviation (shared part plus the probit residual
# of 1), so they read as the correlation between the species' liability
# and the axis.
ci_std <- loading_ci(fit, loading_scale = "standardized")
ci_std[!ci_std$pinned, c("trait", "axis", "estimate", "lower", "upper")]
#> trait axis estimate lower upper
#> 1 sp1 LV1 0.6937971 0.48925276 0.8983414
#> 3 sp3 LV1 0.5052524 0.28093599 0.7295688
#> 5 sp5 LV1 0.5861732 0.36979553 0.8025509
#> 7 sp7 LV1 0.2193042 -0.02615639 0.4647648
#> 10 sp2 LV2 0.6347611 0.43026047 0.8392617
#> 12 sp4 LV2 0.6272572 0.41750598 0.8370083
#> 14 sp6 LV2 0.5957332 0.38793326 0.8035331
#> 16 sp8 LV2 0.3886037 0.15863732 0.6185700
# Symmetric intervals can spill past +/-1 on this scale; "wald_asym"
# builds them on Fisher's z scale, so they stay inside (-1, 1).
ci_asym <- loading_ci(fit, method = "wald_asym")
ci_asym[!ci_asym$pinned, c("trait", "axis", "estimate", "lower", "upper")]
#> trait axis estimate lower upper
#> 1 sp1 LV1 0.6937971 0.43078353 0.8481776
#> 3 sp3 LV1 0.5052524 0.24973017 0.6949897
#> 5 sp5 LV1 0.5861732 0.32941642 0.7622061
#> 7 sp7 LV1 0.2193042 -0.03492307 0.4468740
#> 10 sp2 LV2 0.6347611 0.38579307 0.7975532
#> 12 sp4 LV2 0.6272572 0.37228241 0.7941953
#> 14 sp6 LV2 0.5957332 0.34907257 0.7651932
#> 16 sp8 LV2 0.3886037 0.13838913 0.5921858