Approximate Leave-One-Out Cross-Validation for INLAvaan Models

Computes leave-one-out (LOO) cross-validation for a fitted INLAvaan model from a single fit, with no refitting and no sampling, via a Taylor approximation of the case-deletion posterior around the Laplace summary. Single-level models are scored per subject (leave-one-subject-out, LOSO); two-level models are scored per cluster (leave-one-cluster-out, LOCO).

Usage

loo(x, ...)

# S3 method for class 'INLAvaan'
loo(
  x,
  type = c("auto", "loso", "loco"),
  units = NULL,
  second_order = TRUE,
  theta = NULL,
  Omega = NULL,
  cores = NULL,
  verbose = FALSE,
  ...
)

# S3 method for class 'inlavaan_internal'
loo(
  x,
  type = c("auto", "loso", "loco"),
  units = NULL,
  second_order = TRUE,
  theta = NULL,
  Omega = NULL,
  cores = NULL,
  verbose = FALSE,
  ...
)

add_loo(object, cores = NULL, verbose = FALSE)

# S3 method for class 'inlavaan_loo'
print(x, ...)

# S3 method for class 'inlavaan_loo'
summary(object, ...)

Arguments

x

A fitted INLAvaan object (or its inlavaan_internal list).

...

Not used, beyond catching the deprecated argument name Sigma (see Omega).

type

Unit type: "auto" (default) resolves to "loso" (per-subject) for single-level models and "loco" (per-cluster, marginal predictive) for two-level models. "loco" cannot be forced on a model without clusters; "loso" on a two-level model scores the conditional (leave-one-unit-out) predictive instead (with a warning; see Details).

units

Optional integer vector of unit indices to score; defaults to all units. For LOSO these are case numbers (row numbers of the analysed dataset); for LOCO, cluster positions.

second_order

Logical; compute the second-order correction (default TRUE). FALSE skips the Hessian stage and reports first-order estimates, which cannot be compared across models of different dimension (see Details).

theta, Omega

Optional posterior mean vector and covariance matrix (in the unconstrained parameter space, as stored in theta_star and Sigma_theta) at which to evaluate the LOO instead of the fit's own Laplace summary. See Details.

cores

Number of cores for the Hessian stage. The default NULL runs serially; parallelism must be requested explicitly.

verbose

Logical; print progress (default FALSE).

object

A fitted INLAvaan object, or an inlavaan_loo result for summary().

Value

An object of class inlavaan_loo: a list with elements

per_unit

Data frame of pointwise results, one row per unit:

unit

Case number for LOSO, cluster position for LOCO.

group

Group membership (multigroup fits only).

nobs

1 for LOSO, the cluster size for LOCO.

l_star

Unit log-likelihood at the summary.

score_norm

Norm of the unit score \(s_u\).

lpd_1, lpd_2

Pointwise log predictive density, at first and second order.

log_cpo_1, log_cpo_2

Pointwise LOO contributions, at first and second order.

det_term

\(\tfrac12 \log |I + \Omega H_u|\), the determinant term of the second-order score.

k_max, k_min, k_sum

Leverage diagnostics (see Details).

k_ssq

\(\mathrm{tr}[(\Omega H_u)^2]\), consumed by the closed-form waic() penalty.

ok

Whether the second-order \(\log \mathrm{CPO}\) exists.

estimates

Matrix with rows elpd_loo, p_loo, looic and columns Estimate, SE, at the highest order available to each.

elpd_1, elpd_2, se_1, se_2, p_loo_1, p_loo_2

First- and second-order aggregates; the second-order ones are NA when any \(\log \mathrm{CPO}_u^{(2)}\) does not exist.

elpd_gap, pd_trace

The two sides of the curvature check (see Details); both are NA at first order and partial totals under a units subset.

type, flavour, n_units, n_groups, n_ok, n_lpd_ok, second_order, use_second, theta_overridden

Metadata.

summary() is an alias for print(): it prints the same output and returns the result invisibly.

add_loo() returns a copy of object with the LOO and WAIC results stored alongside the fit (the input object is unchanged); reassign it, e.g. fit <- add_loo(fit). Only the default LOO is stored, so the stored results always match loo(fit) and waic(fit).

Details

For a unit \(u\) (a subject for LOSO, a cluster for LOCO) with log-likelihood contribution \(\ell_u(\theta)\), score \(s_u\) and Hessian \(H_u\) evaluated at the posterior summary \((\theta^*, \Omega)\), the log conditional predictive ordinate is approximated to first and second order by $$\log \mathrm{CPO}_u^{(1)} = \ell_u - \tfrac{1}{2} s_u' \Omega s_u,$$ $$\log \mathrm{CPO}_u^{(2)} = \ell_u - \tfrac{1}{2} s_u' (\Omega^{-1} + H_u)^{-1} s_u + \tfrac{1}{2} \log |I + \Omega H_u|.$$ The reported elpd_loo is the sum of the second-order terms (first-order when second_order = FALSE), with standard error \(\sqrt{n \, \mathrm{var}(\log \mathrm{CPO}_u)}\) and looic \(= -2 \, \mathrm{elpd}\). p_loo is the loo package's effective number of parameters, \(p_{\mathrm{loo}} = \sum_u (\mathrm{lpd}_u - \log \mathrm{CPO}_u)\), where \(\mathrm{lpd}_u\) is the analogous Taylor approximation of the full-posterior pointwise log predictive density – the same definition loo::loo() reports, and not the \(p_D\) of the DIC.

\(\log \mathrm{CPO}_u^{(2)}\) exists exactly when \(\Omega^{-1} + H_u\) is positive definite (recorded in per_unit$ok) and \(\mathrm{lpd}_u^{(2)}\) exactly when \(\Omega^{-1} - H_u\) is. A unit failing the former drops every estimate to first order over all units (with a warning), while one failing only the latter contributes its first-order difference to p_loo.

The leverages k_max, k_min, and k_sum read these conditions off the spectrum of \(-\Omega H_u\) (k_max < 1, k_min > -1), with k_sum summing across units to the trace form of \(p_D\). p_loo (cross-product form) and \(p_D\) (second-derivative form) agree only in the correct-specification limit, and printing a result reports the first-to-second-order gap against its limit \(p_D/2\) as a free check on the Taylor truncation. Because the first-order elpd overstates the truth by \(\tfrac12 p_D\) in the limit, keep second_order = TRUE whenever models of different dimension are compared.

type = "auto" resolves to the marginal per-cluster "loco" for two-level fits and per-subject "loso" otherwise. Forcing "loso" on a two-level model scores the conditional predictive of Merkle, Furr & Rabe-Hesketh (2019) instead, and warns. Multigroup units are scored against their own group's implied moments and identified by case number, so compare() pairs them across fits.

The score follows the fitted likelihood's treatment of exogenous covariates, i.e. joint under fixed.x = FALSE, conditional under fixed.x = TRUE (recorded as "joint" or "conditional" in the result's flavour field), and the two flavours are never comparable (compare() refuses to mix them).

Supplying theta/Omega evaluates the LOO at an arbitrary Gaussian posterior summary (a singular Omega is restricted to its non-degenerate block), the building block for refit-free submodel scoring. Sigma is the deprecated former name of Omega, still accepted through ... with a warning.

The LOO is stored with the fit when inlavaan()'s test includes "loo" or "waic" (e.g. test = "full"), which also stores the WAIC (see waic()) from the same Taylor pass, or afterwards with add_loo(); loo(fit) with default arguments then returns the stored result. Under the default test = "standard" nothing is stored and loo(fit) computes it on demand.

The default cores = NULL runs serially, and cores > 1 parallelises the Hessian stage. Supported models are continuous-indicator models fitted with the ML estimator, single- or two-level, single-group or multigroup (multigroup two-level models are not supported yet).

References

Alhyari, M., Jamil, H., Montcho, H., & Rue, H. (2026). Deterministic leave-one-cluster-out cross-validation for multilevel Bayesian structural equation models. arXiv. (Preprint forthcoming; placeholder.)

Merkle, E. C., Furr, D., & Rabe-Hesketh, S. (2019). Bayesian comparison of latent variable models: Conditional versus marginal likelihoods. Psychometrika, 84(3), 802–829. https://doi.org/10.1007/s11336-019-09679-0

See also

fitmeasures(), compare(), inlavaan()

Examples

# \donttest{
HS.model <- "
  visual  =~ x1 + x2 + x3
  textual =~ x4 + x5 + x6
  speed   =~ x7 + x8 + x9
"
utils::data("HolzingerSwineford1939", package = "lavaan")
fit <- acfa(HS.model, HolzingerSwineford1939, meanstructure = TRUE)
#> ℹ Mode finding and Hessian computation.
#> ✔ Posterior mode and Hessian. [174ms]
#> 
#> ℹ Performing VB correction.
#> ✔ VB correction; mean |δ| = 0.133σ. [189ms]
#> 
#> ⠙ Fitting 0/30 skew-normal marginals.
#> ⠹ Fitting 15/30 skew-normal marginals.
#> ✔ Fit 30/30 skew-normal marginals. [888ms]
#> 
#> ⠙ Posterior sampling and summarising.
#> ✔ Summarise 1000 posterior draws. [531ms]
#> 
#> ℹ Fit measures: PPP, DIC.

# Leave-one-subject-out (LOSO) from the single fit -- no refitting
res <- loo(fit)
res
#> ── Leave-one-subject-out ───────────────────────── 301 subjects, second-order ──
#> 
#>          Estimate   SE
#> elpd_loo  -3769.1 43.0
#> p_loo        32.5  2.2
#> looic      7538.2 86.0
#> 
#> ── Curvature check ─────────────────────────────────────────────────────────────
#> 
#>   first-to-second-order gap        15.4
#>   pD/2 (trace)                     14.5
#>   excess over pD/2 (trace)        +6.1%
#> 
#> ℹ The gap approaches pD/2 (trace) from above. A large excess says the
#>   second-order expansion has not settled over the sample.
head(res$per_unit)
#>   unit nobs    l_star score_norm     lpd_1     lpd_2 log_cpo_1 log_cpo_2
#> 1    1    1 -17.30919   6.646703 -17.19476 -17.24151 -17.42362 -17.47403
#> 2    2    1 -13.78169   5.442203 -13.69904 -13.76970 -13.86434 -13.93963
#> 3    3    1 -11.11040   3.789743 -11.07449 -11.16498 -11.14631 -11.23920
#> 4    4    1 -10.25262   2.580054 -10.23857 -10.28526 -10.26667 -10.31439
#> 5    5    1 -10.70128   2.974402 -10.68783 -10.73479 -10.71473 -10.76276
#> 6    6    1 -13.37858   5.008821 -13.31650 -13.40400 -13.44066 -13.53238
#>      det_term      k_max        k_min      k_sum       k_ssq   ok
#> 1 -0.04877732 0.03773964 -0.041266991 0.09404961 0.006990423 TRUE
#> 2 -0.07270997 0.05910469 -0.049278162 0.14094433 0.008831504 TRUE
#> 3 -0.09217179 0.04153084 -0.009791359 0.18193722 0.004724893 TRUE
#> 4 -0.04763225 0.02969695 -0.010269099 0.09422303 0.002056503 TRUE
#> 5 -0.04796440 0.03026567 -0.011352189 0.09484979 0.002130082 TRUE
#> 6 -0.08977773 0.06593768 -0.023150767 0.17537797 0.008112617 TRUE

# Score a submodel without refitting: condition the Laplace summary on the
# visual ~~ speed covariance being zero, then evaluate at that summary
int <- get_inlavaan_internal(fit)
theta <- int$theta_star
Omega <- int$Sigma_theta
p <- which(names(coef(fit)) == "visual~~speed")
theta_c <- theta - Omega[, p] * (theta[p] / Omega[p, p])
Omega_c <- Omega - tcrossprod(Omega[, p]) / Omega[p, p]
loo(fit, theta = theta_c, Omega = Omega_c)
#> ── Leave-one-subject-out ───────────────────────── 301 subjects, second-order ──
#> 
#>          Estimate   SE
#> elpd_loo  -3786.1 45.0
#> p_loo        34.1  2.5
#> looic      7572.1 90.0
#> 
#> ℹ Evaluated at a user-supplied (theta, Sigma) summary.
#> 
#> ── Curvature check ─────────────────────────────────────────────────────────────
#> 
#>   first-to-second-order gap        16.3
#>   pD/2 (trace)                     15.3
#>   excess over pD/2 (trace)        +6.5%
#> 
#> ℹ The gap approaches pD/2 (trace) from above. A large excess says the
#>   second-order expansion has not settled over the sample.

# Two-level models are scored per cluster (LOCO) automatically
utils::data("Demo.twolevel", package = "lavaan")
model2l <- "
  level: 1
    fw =~ y1 + y2 + y3
    fw ~ x1 + x2 + x3
  level: 2
    fb =~ y1 + y2 + y3
    fb ~ w1 + w2
"
fit2l <- asem(model2l, Demo.twolevel, cluster = "cluster",
              meanstructure = TRUE, fixed.x = FALSE)
#> ℹ Mode finding and Hessian computation.
#> ✔ Posterior mode and Hessian. [895ms]
#> 
#> ℹ Performing VB correction.
#> ✔ VB correction; mean |δ| = 0.050σ. [843ms]
#> 
#> ⠙ Fitting 0/34 skew-normal marginals.
#> ⠹ Fitting 16/34 skew-normal marginals.
#> ⠸ Fitting 31/34 skew-normal marginals.
#> ✔ Fit 34/34 skew-normal marginals. [6.5s]
#> 
#> ⠙ Posterior sampling and summarising.
#> ✔ Summarise 1000 posterior draws. [1.2s]
#> 
#> ℹ Fit measures: PPP, DIC.
loo(fit2l)
#> ── Leave-one-cluster-out ───────────────────────── 200 clusters, second-order ──
#> 
#>          Estimate     SE
#> elpd_loo -23344.2  731.4
#> p_loo        34.5    2.1
#> looic     46688.4 1462.9
#> 
#> ── Curvature check ─────────────────────────────────────────────────────────────
#> 
#>   first-to-second-order gap        17.4
#>   pD/2 (trace)                     16.7
#>   excess over pD/2 (trace)        +4.3%
#> 
#> ℹ The gap approaches pD/2 (trace) from above. A large excess says the
#>   second-order expansion has not settled over the sample.
# }