Package {INLAvaan}

Type:	Package
Title:	Approximate Bayesian Latent Variable Analysis
Version:	0.3.0
Description:	Implements approximate Bayesian inference for Structural Equation Models (SEM) using a custom adaptation of the Integrated Nested Laplace Approximation (Rue et al., 2009) <doi:10.1111/j.1467-9868.2008.00700.x> as described in Jamil and Rue (2026a) <doi:10.48550/arXiv.2603.25690>. Provides a computationally efficient alternative to Markov Chain Monte Carlo (MCMC) for Bayesian estimation, allowing users to fit latent variable models using the 'lavaan' syntax. See also the companion paper on implementation and workflows, Jamil and Rue (2026b) <doi:10.48550/arXiv.2604.00671>.
License:	GPL (≥ 3)
URL:	https://inlavaan.haziqj.ml/, https://github.com/haziqj/INLAvaan
BugReports:	https://github.com/haziqj/INLAvaan/issues
Depends:	R (≥ 3.5)
Imports:	cli, graphics, lavaan (≥ 0.6-19), methods, parallel, stats, utils
Suggests:	blavaan, ggplot2, knitr, lme4, numDeriv, qrng, quarto, sn, testthat (≥ 3.0.0), ucminf
VignetteBuilder:	quarto
Config/Needs/website:	rmarkdown
Config/roxygen2/markdown:	TRUE
Config/roxygen2/version:	8.0.0
Config/testthat/edition:	3
Encoding:	UTF-8
NeedsCompilation:	yes
Packaged:	2026-07-11 11:30:51 UTC; haziqj
Author:	Haziq Jamil [aut, cre, cph], Håvard Rue [ctb] (Statistical and computational methodology), Alvin Bong [ctb] (Initial site build)
Maintainer:	Haziq Jamil <haziq.jamil@gmail.com>
Repository:	CRAN
Date/Publication:	2026-07-11 12:00:02 UTC

INLAvaan: Approximate Bayesian Latent Variable Analysis

Description

Implements approximate Bayesian inference for Structural Equation Models (SEM) using a custom adaptation of the Integrated Nested Laplace Approximation (Rue et al., 2009) doi:10.1111/j.1467-9868.2008.00700.x as described in Jamil and Rue (2026a) doi:10.48550/arXiv.2603.25690. Provides a computationally efficient alternative to Markov Chain Monte Carlo (MCMC) for Bayesian estimation, allowing users to fit latent variable models using the 'lavaan' syntax. See also the companion paper on implementation and workflows, Jamil and Rue (2026b) doi:10.48550/arXiv.2604.00671.

Main features

• acfa(): Approximate Confirmatory Factor Analysis.

• asem(): Approximate Structural Equation Modelling.

• agrowth(): Approximate Latent Growth Curve models.

Model specifications

Supports advanced 'lavaan' syntax features, including:

• Equality constraints

• Defined parameters (e.g., ⁠:=⁠ operator for indirect effects)

• Flexible prior specifications

Methods for INLAvaan objects

After fitting a model an INLAvaan object is returned. The following S4 methods are available. See INLAvaan for the class definition.

• Summaries and parameter estimates: summary(), coef(), vcov(), standardisedsolution()

• Fit assessment and model comparison: fitmeasures(), bfit_indices(), compare(), diagnostics(), timing()

• Posterior inference and simulation: predict(), sampling(), simulate()

• Visualisation: plot()

Online vignettes

The package website contains comprehensive examples covering:

• Confirmatory Factor Analysis (CFA)

• Structural Equation Models (SEM)

• Latent Growth Curve Models

• Multigroup and Invariance Testing

• Mediation Analysis

Author(s)

Maintainer: Haziq Jamil haziq.jamil@gmail.com (ORCID) [copyright holder]

Authors:

• Haziq Jamil haziq.jamil@gmail.com (ORCID) [copyright holder]

Other contributors:

• Håvard Rue (ORCID) (Statistical and computational methodology) [contributor]

• Alvin Bong (Initial site build) [contributor]

See Also

Useful links:

• https://inlavaan.haziqj.ml/

• https://github.com/haziqj/INLAvaan

• Report bugs at https://github.com/haziqj/INLAvaan/issues

Class For Representing a (Fitted) Latent Variable Model

Description

This is a class that extends the lavaan class. Several S4 methods are available.

Usage

## S4 method for signature 'INLAvaan'
coef(object, ...)

## S4 method for signature 'INLAvaan'
nobs(object, ...)

## S4 method for signature 'INLAvaan'
show(object)

## S4 method for signature 'INLAvaan'
summary(
  object,
  header = TRUE,
  fit.measures = TRUE,
  estimates = TRUE,
  ci = TRUE,
  standardized = FALSE,
  standardised = standardized,
  rsquare = FALSE,
  postmedian = FALSE,
  postmode = FALSE,
  nmad = TRUE,
  kld = FALSE,
  vb_shift = FALSE,
  priors = TRUE,
  nd = 3L,
  ...
)

Arguments

Slots

external

A list containing an inlavaan_internal object.

See Also

lavaan, inlavaan(), acfa(), asem(), agrowth()

Examples

HS.model <- "
  visual  =~ x1 + x2 + x3
  textual =~ x4 + x5 + x6
  speed   =~ x7 + x8 + x9
"
utils::data("HolzingerSwineford1939", package = "lavaan")
fit <- acfa(HS.model, HolzingerSwineford1939, std.lv = TRUE, nsamp = 100,
            test = "none", verbose = FALSE)

# Print basic info
fit

# Detailed summary
summary(fit)

# Extract coefficients
coef(fit)

Fit an Approximate Bayesian Confirmatory Factor Analysis Model

Description

Fit an Approximate Bayesian Confirmatory Factor Analysis Model

Usage

acfa(
  model,
  data,
  dp = priors_for(),
  test = "standard",
  vb_correction = TRUE,
  marginal_method = c("skewnorm", "asymgaus", "marggaus", "sampling"),
  marginal_correction = c("shortcut", "shortcut_fd", "hessian", "none"),
  nsamp = 1000,
  samp_copula = TRUE,
  sn_fit_ngrid = 21,
  sn_fit_logthresh = -6,
  sn_fit_temp = 1,
  sn_fit_sample = TRUE,
  control = list(),
  verbose = TRUE,
  debug = FALSE,
  add_priors = TRUE,
  optim_method = c("nlminb", "ucminf", "optim"),
  numerical_grad = FALSE,
  cores = NULL,
  ...
)

Arguments

Details

The acfa() function is a wrapper for the more general inlavaan() function, using the following default arguments:

• int.ov.free = TRUE

• int.lv.free = FALSE

• auto.fix.first = TRUE (unless std.lv = TRUE)

• auto.fix.single = TRUE

• auto.var = TRUE

• auto.cov.lv.x = TRUE

• auto.efa = TRUE

• auto.th = TRUE

• auto.delta = TRUE

• auto.cov.y = TRUE

For further information regarding these arguments, please refer to the lavaan::lavOptions() documentation.

Value

An S4 object of class INLAvaan which is a subclass of the lavaan class.

See Also

Typically, users will interact with the specific latent variable model functions instead, including acfa(), asem(), and agrowth().

Examples

# The famous Holzinger and Swineford (1939) example
HS.model <- "
  visual  =~ x1 + x2 + x3
  textual =~ x4 + x5 + x6
  speed   =~ x7 + x8 + x9
"
utils::data("HolzingerSwineford1939", package = "lavaan")

# Fit a CFA model with standardised latent variables
fit <- acfa(HS.model, data = HolzingerSwineford1939, std.lv = TRUE, nsamp = 100)
summary(fit)

Fit an Approximate Bayesian Growth Curve Model

Description

Fit an Approximate Bayesian Growth Curve Model

Usage

agrowth(
  model,
  data,
  dp = priors_for(),
  test = "standard",
  vb_correction = TRUE,
  marginal_method = c("skewnorm", "asymgaus", "marggaus", "sampling"),
  marginal_correction = c("shortcut", "shortcut_fd", "hessian", "none"),
  nsamp = 1000,
  samp_copula = TRUE,
  sn_fit_ngrid = 21,
  sn_fit_logthresh = -6,
  sn_fit_temp = 1,
  sn_fit_sample = TRUE,
  control = list(),
  verbose = TRUE,
  debug = FALSE,
  add_priors = TRUE,
  optim_method = c("nlminb", "ucminf", "optim"),
  numerical_grad = FALSE,
  cores = NULL,
  ...
)

Arguments

Details

The asem() function is a wrapper for the more general inlavaan() function, using the following default arguments:

• meanstructure = TRUE

• int.ov.free = FALSE

• int.lv.free = TRUE

• auto.fix.first = TRUE (unless std.lv = TRUE)

• auto.fix.single = TRUE

• auto.var = TRUE

• auto.cov.lv.x = TRUE

• auto.efa = TRUE

• auto.th = TRUE

• auto.delta = TRUE

• auto.cov.y = TRUE

Value

An S4 object of class INLAvaan which is a subclass of the lavaan class.

See Also

Typically, users will interact with the specific latent variable model functions instead, including acfa(), asem(), and agrowth().

Examples

# Linear growth model with a time-varying covariate
mod <- "
  # Intercept and slope with fixed coefficients
    i =~ 1*t1 + 1*t2 + 1*t3 + 1*t4
    s =~ 0*t1 + 1*t2 + 2*t3 + 3*t4

  # (Latent) regressions
    i ~ x1 + x2
    s ~ x1 + x2

  # Time-varying covariates
    t1 ~ c1
    t2 ~ c2
    t3 ~ c3
    t4 ~ c4
"
utils::data("Demo.growth", package = "lavaan")
str(Demo.growth)

fit <- agrowth(mod, data = Demo.growth, nsamp = 100)
summary(fit)

Convert function to single string

Description

Convert function to single string

Usage

as_fun_string(f)

Arguments

Value

A single character vector representing the function.

Examples

f <- function(x) { x^2 + 1 }
as_fun_string(f)

Fit an Approximate Bayesian Structural Equation Model

Description

Fit an Approximate Bayesian Structural Equation Model

Usage

asem(
  model,
  data,
  dp = priors_for(),
  test = "standard",
  vb_correction = TRUE,
  marginal_method = c("skewnorm", "asymgaus", "marggaus", "sampling"),
  marginal_correction = c("shortcut", "shortcut_fd", "hessian", "none"),
  nsamp = 1000,
  samp_copula = TRUE,
  sn_fit_ngrid = 21,
  sn_fit_logthresh = -6,
  sn_fit_temp = 1,
  sn_fit_sample = TRUE,
  control = list(),
  verbose = TRUE,
  debug = FALSE,
  add_priors = TRUE,
  optim_method = c("nlminb", "ucminf", "optim"),
  numerical_grad = FALSE,
  cores = NULL,
  ...
)

Arguments

Details

The asem() function is a wrapper for the more general inlavaan() function, using the following default arguments:

• int.ov.free = TRUE

• int.lv.free = FALSE

• auto.fix.first = TRUE (unless std.lv = TRUE)

• auto.fix.single = TRUE

• auto.var = TRUE

• auto.cov.lv.x = TRUE

• auto.efa = TRUE

• auto.th = TRUE

• auto.delta = TRUE

• auto.cov.y = TRUE

For further information regarding these arguments, please refer to the lavaan::lavOptions() documentation.

Value

An S4 object of class INLAvaan which is a subclass of the lavaan class.

See Also

Typically, users will interact with the specific latent variable model functions instead, including acfa(), asem(), and agrowth().

Examples

# The industrialization and Political Democracy Example from Bollen (1989), page
# 332
model <- "
  # Latent variable definitions
     ind60 =~ x1 + x2 + x3
     dem60 =~ y1 + a*y2 + b*y3 + c*y4
     dem65 =~ y5 + a*y6 + b*y7 + c*y8

  # (Latent) regressions
    dem60 ~ ind60
    dem65 ~ ind60 + dem60

  # Residual correlations
    y1 ~~ y5
    y2 ~~ y4 + y6
    y3 ~~ y7
    y4 ~~ y8
    y6 ~~ y8
"
utils::data("PoliticalDemocracy", package = "lavaan")

fit <- asem(model, PoliticalDemocracy, test = "none")
summary(fit)

Bayesian Fit Indices

Description

Compute posterior distributions of Bayesian fit indices for an INLAvaan model, analogous to blavaan::blavFitIndices().

Usage

bfit_indices(
  object,
  baseline.model = NULL,
  rescale = c("devM", "MCMC"),
  nsamp = NULL,
  samp_copula = TRUE
)

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

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

Arguments

Value

An S3 object of class "bfit_indices" containing:

indices

Named list of numeric vectors (one per posterior sample) for each computed fit index.

details

List with chisq (per-sample deviance), df, pD, rescale, and nsamp.

Use summary() to obtain a table of posterior summaries (Mean, SD, quantiles, Mode) for each index.

See Also

lavaan::fitMeasures(), blavaan::blavFitIndices(), fitmeasures(), compare()

Examples

HS.model <- "
  visual  =~ x1 + x2 + x3
  textual =~ x4 + x5 + x6
  speed   =~ x7 + x8 + x9
"
utils::data("HolzingerSwineford1939", package = "lavaan")
fit <- acfa(HS.model, HolzingerSwineford1939, std.lv = TRUE, nsamp = 100,
            verbose = FALSE)

# Absolute fit indices
bf <- bfit_indices(fit)
bf
summary(bf)

Compare Bayesian Models Fitted with INLAvaan

Description

Compare two or more Bayesian SEM fitted with INLAvaan, reporting model-fit statistics and (optionally) fit indices side by side.

Usage

compare(x, y, ..., fit.measures = NULL, loo = FALSE)

## S4 method for signature 'INLAvaan'
compare(x, y, ..., fit.measures = NULL, loo = FALSE)

## S4 method for signature 'INLAvaan'
anova(object, ...)

Arguments

Details

The first argument x serves as the baseline (null) model. All models (including the baseline) appear in the comparison table. The baseline is also passed to fitMeasures() when incremental fit indices (BCFI, BTLI, BNFI) are requested via fit.measures.

The default table always includes:

• npar: Number of free parameters.

• Marg.Loglik: Approximated marginal log-likelihood.

• logBF: Natural-log Bayes factor relative to the best model.

• DIC / pD: Deviance Information Criterion and effective number of parameters (when test != "none" was used during fitting).

Set fit.measures to a character vector of measure names (anything returned by fitMeasures()) to append extra columns. Use fit.measures = "all" to include every available measure.

Set loo = TRUE to compare models by leave-one-out cross-validation (see loo()). This appends ELPD / SE (the second-order Taylor expected log predictive density and its standard error), p_loo, and, against the best-ELPD model, the difference elpd_diff with its paired standard error se_diff computed from the pointwise contributions (the appropriate uncertainty for nested or same-data comparisons). The table is then sorted by descending ELPD. All models must be fitted to the same data with matching units; units are paired by id rather than by row order, so fits that stack groups differently – a pooled fit against a multigroup fit, or multigroup fits with different group orderings – still pair up unit by unit. For missing-data (FIML) fits, "the same data" also means the same observed entries: each unit is scored on the entries it has, so comparisons require identical missingness patterns across models. All models must also share the score flavour (see loo()): mixing fits with modelled covariates (fixed.x = FALSE, joint scores) and fixed covariates (fixed.x = TRUE, conditional scores) is refused. Joint scores additionally require identical variable sets across models, while conditional scores require only matching outcome variables – covariate sets may differ, which is the covariate-selection setting. Stored LOO results (test = "loo" or add_loo()) are reused.

anova() is disabled for INLAvaan fits – there is no direct Bayesian analogue of the classical likelihood-ratio test – and points here instead.

Value

A data frame of class compare.inlavaan_internal containing model fit statistics, sorted by descending marginal log-likelihood (or by descending ELPD when loo = TRUE).

References

https://lavaan.ugent.be/tutorial/groups.html

See Also

fitmeasures(), bfit_indices()

Examples

# Model comparison on multigroup analysis (measurement invariance)
HS.model <- "
  visual  =~ x1 + x2 + x3
  textual =~ x4 + x5 + x6
  speed   =~ x7 + x8 + x9
"
utils::data("HolzingerSwineford1939", package = "lavaan")

# Configural invariance
fit1 <- acfa(HS.model, data = HolzingerSwineford1939, group = "school")

# Weak invariance
fit2 <- acfa(
  HS.model,
  data = HolzingerSwineford1939,
  group = "school",
  group.equal = "loadings"
)

# Strong invariance
fit3 <- acfa(
  HS.model,
  data = HolzingerSwineford1939,
  group = "school",
  group.equal = c("intercepts", "loadings")
)

# Compare models (fit1 = configural = baseline, always first argument)
compare(fit1, fit2, fit3)

# With extra fit measures
compare(fit1, fit2, fit.measures = c("BRMSEA", "BMc"))

# With incremental indices (baseline = fit1, passed to fitMeasures())
compare(fit1, fit2, fit3, fit.measures = c("BCFI", "BTLI"))

Density of a Beta distribution on a bounded interval

Description

Density of a Beta distribution on a bounded interval

Usage

dbeta_box(x, shape1, shape2, a, b, log = FALSE)

Arguments

Value

A numeric vector of density values.

Examples

# Beta(2,5) on (0,100)
x <- seq(0, 100, length.out = 100)
y <- dbeta_box(x, shape1 = 2, shape2 = 5, a = 0, b = 100)
plot(x, y, type = "l", main = "Beta(2,5) on (0,100)")

# Beta(1,1) i.e. uniform on (-1, 1)
x <- seq(-1, 1, length.out = 100)
y <- dbeta_box(x, shape1 = 1, shape2 = 1, a = -1, b = 1)
plot(x, y, type = "l", main = "Beta(1,1) on (-1,1)")

Deviance for INLAvaan Models

Description

Extract the (Bayesian) deviance of a fitted INLAvaan model. Unlike lavaan, which has no deviance() method, this follows the BUGS/JAGS/Stan convention: "deviance" is -2 times the log-likelihood, summarised over the posterior.

Usage

## S3 method for class 'INLAvaan'
deviance(object, type = c("mean", "plugin"), ...)

Arguments

Details

\bar{D} and \hat{D} are the two ingredients of the Deviance Information Criterion, DIC = \bar{D} + p_D where p_D = \bar{D} - \hat{D} is the effective number of parameters. Use compare() to compare models by DIC (or Bayes factors, or LOO/WAIC) rather than comparing raw deviances directly.

Value

A length-one numeric of class inlavaan_deviance, with the effective number of parameters (pD) and DIC attached as attributes.

See Also

logLik(), compare()

Examples

HS.model <- "
  visual  =~ x1 + x2 + x3
  textual =~ x4 + x5 + x6
  speed   =~ x7 + x8 + x9
"
utils::data("HolzingerSwineford1939", package = "lavaan")
fit <- acfa(HS.model, HolzingerSwineford1939, std.lv = TRUE, nsamp = 100,
            test = "standard", verbose = FALSE)

deviance(fit)
attr(deviance(fit), "DIC")

Convergence and Approximation Diagnostics for INLAvaan Models

Description

Extract convergence and approximation-quality diagnostics from a fitted INLAvaan model.

Usage

diagnostics(object, ...)

## S4 method for signature 'INLAvaan'
diagnostics(object, type = c("global", "param"), ...)

Arguments

Details

Global diagnostics (type = "global"):

npar

Number of free parameters.

nsamp

Number of posterior samples drawn.

converged

1 if the optimiser converged, 0 otherwise.

iterations

Number of optimiser iterations.

grad_inf

L-infinity norm of the analytic gradient at the mode (max |grad|). Should be ~0 at convergence.

grad_inf_rel

Relative L-infinity norm of the analytic gradient (max |grad| / (|par| + 1e-6)).

grad_l2

L2 (Euclidean) norm of the analytic gradient at the mode.

mode_shift_max

Maximum, across parameters, of the Newton step at the reported mode in posterior-SD units (max |\Sigma_\theta grad| / se). Unlike the raw gradient norms this is scale-free: it estimates how far the reported mode sits from the true posterior mode relative to the posterior uncertainty. Should be ~0 at convergence.

hess_cond

Condition number of the Hessian (precision matrix) computed from \Sigma_\theta. Large values indicate near-singularity.

vb_kld_global

Global KL divergence from the VB mean correction (NA if VB correction was not applied).

vb_applied

1 if VB correction was applied, 0 otherwise.

kld_max

Maximum per-parameter KL divergence from the VB correction.

kld_mean

Mean per-parameter KL divergence.

nmad_max

Maximum normalised max-absolute-deviation across marginals (skew-normal method only; NA otherwise).

nmad_mean

Mean NMAD across marginals.

Per-parameter diagnostics (type = "param"): A data frame with columns:

param

Parameter name.

grad

Analytic gradient of the negative log-posterior at the mode. Should be ~0 at convergence.

grad_num

Numerical (finite-difference) gradient at the mode. Should agree with grad; large discrepancies indicate a bug in the analytic gradient.

grad_diff

Difference grad_num - grad: should be ~0.

grad_abs

Absolute analytic gradient.

grad_rel

Relative analytic gradient |grad| / (|par| + 1e-6).

mode_shift_sigma

Newton step at the reported mode in posterior-SD units. Should be ~0 at convergence.

kld

Per-parameter KL divergence from the VB correction.

vb_shift

VB correction shift (in original scale).

vb_shift_sigma

VB shift in units of posterior SD.

nmad

Normalised max-absolute-deviation of the skew-normal fit (NA when not using the skewnorm method).

Fit-time warnings: inlavaan() runs these checks once at the end of every fit and emits a single consolidated warning (condition class "inlavaan_diagnostics_warning") when any of them look off: the optimiser did not converge, mode_shift_max exceeds 0.1, any marginal has NMAD above 0.1, the VB correction shifted a posterior mean by more than 1 posterior SD, or the Hessian condition number exceeds 1e8. A healthy fit stays silent. Silence the check with suppressWarnings(), or selectively by handling the condition class.

Value

For type = "global", a named numeric vector (class "diagnostics.INLAvaan"). For type = "param", a data frame (class c("diagnostics.INLAvaan.param", "data.frame")).

See Also

timing(), fitmeasures(), plot()

Examples

HS.model <- "
  visual  =~ x1 + x2 + x3
  textual =~ x4 + x5 + x6
  speed   =~ x7 + x8 + x9
"
utils::data("HolzingerSwineford1939", package = "lavaan")
fit <- acfa(HS.model, HolzingerSwineford1939, std.lv = TRUE, nsamp = 100,
            test = "none", verbose = FALSE)

# Global convergence summary
diagnostics(fit)

# Per-parameter table
diagnostics(fit, type = "param")

The Skew Normal Distribution

Description

Density for the skew normal distribution with location xi, scale omega, and shape alpha.

Usage

dsnorm(x, xi, omega, alpha, logC = 0, log = FALSE)

Arguments

Value

A numeric vector of (log) density values.

References

https://en.wikipedia.org/wiki/Skew_normal_distribution

Examples

x <- seq(-2, 5, length.out = 100)
y <- dsnorm(x, xi = 0, omega = 1, alpha = 5)
plot(x, y, type = "l", main = "Skew Normal Density")

Fit a skew normal distribution to log-density evaluations

Description

Fit a skew normal distribution to log-density evaluations

Usage

fit_skew_normal(x, y, threshold_log_drop = -6, temp = NA)

Arguments

Details

This skew normal fitting function uses a weighted least squares approach to fit the log-density evaluations provided in y at points x. The weights are set to be the density evaluations raised to the power of the temperature parameter k. This has somewhat an interpretation of finding the skew normal fit that minimises the Kullback-Leibler divergence from the true density to it.

In R-INLA, the C code implementation from which this was translated from can be found here.

Value

A list with fitted parameters:

• xi: location parameter

• omega: scale parameter

• alpha: shape parameter

• logC: log-normalization constant

• k: temperature parameter

• rsq: R-squared of the fit

Note that logC and k are not used when fitting from a sample.

Examples

# Fit a SN curve to gamma log-density
logdens <- function(x) dgamma(x, shape = 3, rate = 1, log = TRUE)

x_grid <- seq(0.1, 8, length.out = 21)
y_log  <- sapply(x_grid, logdens)
y_log  <- y_log - max(y_log)  # normalise to have maximum at zero

res <- fit_skew_normal(x_grid, y_log, temp = 10)
unlist(res)

# Compare truth vs skew-normal approximation
x_fine <- seq(0.1, 8, length.out = 200)
y_true <- exp(logdens(x_fine))
y_sn   <- dsnorm(x_fine, xi = res$xi, omega = res$omega, alpha = res$alpha)

plot(x_fine, y_true, type = "n", xlab = "x", ylab = "Density", bty = "n")
polygon(c(x_fine, rev(x_fine)), c(y_true, rep(0, length(x_fine))),
        col = adjustcolor("#131516", 0.25), border = NA)
lines(x_fine, y_sn, col = "#00A6AA", lwd = 2)
legend("topright", legend = c("Truth", "SN Approx."),
       fill = c(adjustcolor("#131516", 0.25), NA), border = NA,
       col = c(NA, "#00A6AA"), lwd = c(NA, 2), lty = c(NA, 1),
       bty = "n")

Fit a skew normal distribution to a sample

Description

Fit a skew normal distribution to a sample

Usage

fit_skew_normal_samp(x)

Arguments

Details

Uses maximum likelihood estimation to fit a skew normal distribution to the provided numeric vector x.

Value

A list with fitted parameters:

• xi: location parameter

• omega: scale parameter

• alpha: shape parameter

• logC: log-normalization constant

• k: temperature parameter

• rsq: R-squared of the fit

Note that logC and k are not used when fitting from a sample.

Examples

x <- rnorm(100, mean = 5, sd = 1)
unlist(fit_skew_normal_samp(x))

Fit Measures for a Latent Variable Model estimated using INLA

Description

Fit Measures for a Latent Variable Model estimated using INLA

Arguments

Value

A named numeric vector of fit measures.

See Also

bfit_indices(), compare(), diagnostics()

Examples

HS.model <- "
  visual  =~ x1 + x2 + x3
  textual =~ x4 + x5 + x6
  speed   =~ x7 + x8 + x9
"
utils::data("HolzingerSwineford1939", package = "lavaan")
fit <- acfa(HS.model, HolzingerSwineford1939, std.lv = TRUE, nsamp = 100,
            verbose = FALSE)

# All available fit measures
fitMeasures(fit)

# Specific measures
fitMeasures(fit, c("npar", "DIC", "pD", "ppp"))

Model-Implied Moments for INLAvaan Models

Description

Extract the model-implied (fitted) sample statistics from a fitted INLAvaan model. As in lavaan and blavaan, the moments are the model-implied covariance matrix (and mean vector, when a mean structure is present) evaluated at the parameter estimates – here the posterior means.

Usage

## S4 method for signature 'INLAvaan'
fitted(object, type = "moments", labels = TRUE, ...)

## S4 method for signature 'INLAvaan'
fitted.values(object, type = "moments", labels = TRUE, ...)

Arguments

Details

This delegates to lavaan's own fitted() machinery, so the return structure matches lavaan exactly. Because INLAvaan stores the posterior means as the point estimates of the fitted object, the implied moments are the posterior-mean model-implied moments (mirroring blavaan).

Value

For type = "moments", a list (or list of lists, for multiple groups) with elements such as cov, mean, and th. For type = "casewise", a numeric matrix of predicted observed-variable values.

See Also

predict(), coef(), fitMeasures()

Examples

HS.model <- "
  visual  =~ x1 + x2 + x3
  textual =~ x4 + x5 + x6
  speed   =~ x7 + x8 + x9
"
utils::data("HolzingerSwineford1939", package = "lavaan")
fit <- acfa(HS.model, HolzingerSwineford1939, std.lv = TRUE, nsamp = 100,
            test = "none", verbose = FALSE)

# Model-implied covariance matrix (posterior means)
fitted(fit)

# Casewise model-predicted observed values
head(fitted(fit, type = "ov"))

Extract the Internal INLAvaan Object

Description

Returns the inlavaan_internal list stored inside a fitted INLAvaan object, optionally extracting a single named element.

Usage

get_inlavaan_internal(object, what)

Arguments

Value

The full inlavaan_internal list, or the named element when what is supplied.

See Also

INLAvaan, diagnostics(), timing()

Examples

HS.model <- "
  visual  =~ x1 + x2 + x3
  textual =~ x4 + x5 + x6
  speed   =~ x7 + x8 + x9
"
utils::data("HolzingerSwineford1939", package = "lavaan")
fit <- acfa(HS.model, HolzingerSwineford1939, std.lv = TRUE, nsamp = 100,
            test = "none", verbose = FALSE)

# Full internal object
int <- get_inlavaan_internal(fit)
names(int)

# Extract a specific element
get_inlavaan_internal(fit, "coefficients")

Fit an Approximate Bayesian Latent Variable Model

Description

This function fits a Bayesian latent variable model by approximating the posterior distributions of the model parameters using various methods, including skew-normal, asymmetric Gaussian, marginal Gaussian, or sampling-based approaches. It leverages the lavaan package for model specification and estimation.

Usage

inlavaan(
  model,
  data,
  model.type = "sem",
  dp = priors_for(),
  test = "standard",
  vb_correction = TRUE,
  marginal_method = c("skewnorm", "asymgaus", "marggaus", "sampling"),
  marginal_correction = c("shortcut", "shortcut_fd", "hessian", "none"),
  nsamp = 1000,
  samp_copula = TRUE,
  sn_fit_ngrid = 21,
  sn_fit_logthresh = -6,
  sn_fit_temp = 1,
  sn_fit_sample = TRUE,
  control = list(),
  verbose = TRUE,
  debug = FALSE,
  add_priors = TRUE,
  optim_method = c("nlminb", "ucminf", "optim"),
  numerical_grad = FALSE,
  start = NULL,
  cores = NULL,
  ...
)

Arguments

Value

An S4 object of class INLAvaan which is a subclass of the lavaan class.

See Also

Typically, users will interact with the specific latent variable model functions instead, including acfa(), asem(), and agrowth().

Examples

# The Holzinger and Swineford (1939) example
HS.model <- "
  visual  =~ x1 + x2 + x3
  textual =~ x4 + x5 + x6
  speed   =~ x7 + x8 + x9
"
utils::data("HolzingerSwineford1939", package = "lavaan")

fit <- inlavaan(
  HS.model,
  data = HolzingerSwineford1939,
  auto.var = TRUE,
  auto.fix.first = TRUE,
  auto.cov.lv.x = TRUE
)
summary(fit)

Helper function to check if two functions are the same

Description

Helper function to check if two functions are the same

Usage

is_same_function(f, g)

Arguments

Value

Logical.

Examples

f1 <- function(x) { x^2 + 1 }
f2 <- function(x) { x^2 + 1 }
is_same_function(f1, f2)  # TRUE

Log-Likelihood for INLAvaan Models

Description

Extract a log-likelihood-flavoured summary from a fitted INLAvaan model. Two distinct quantities are available, deliberately not conflated: the Bayesian marginal log-likelihood (the default) and the classical log-likelihood evaluated at the posterior mean.

Usage

## S4 method for signature 'INLAvaan'
logLik(object, type = c("marginal", "plugin"), ...)

## S4 method for signature 'INLAvaan'
AIC(object, ..., k = 2)

## S4 method for signature 'INLAvaan'
BIC(object, ...)

Arguments

Details

The marginal log-likelihood already integrates over the (Laplace- approximated) posterior, so it is not on the same scale as a classical log-likelihood and should not be passed to AIC()/BIC() – doing so would double-penalise model complexity that the evidence has already accounted for. Use compare() to compare models via Bayes factors, DIC, or LOO/WAIC. The plug-in variant exists for users who specifically want a point-estimate-based classical comparison.

AIC()/BIC() on an INLAvaan fit are themselves disabled (mirroring anova()): both are large-sample asymptotic approximations to quantities INLAvaan already computes directly – AIC approximates predictive accuracy, which loo()/waic() already estimate more rigorously; BIC approximates -2 times the log marginal likelihood, which logLik() already returns directly (up to the Laplace approximation). Point-estimate AIC/BIC remain available for reporting-convention purposes via AIC(logLik(object, type = "plugin")) / BIC(logLik(object, type = "plugin")).

Value

For type = "marginal", a length-one numeric of class inlavaan_logLik that prints with a note on its interpretation. For type = "plugin", a standard "logLik" object.

See Also

deviance(), compare(), loo(), waic()

Examples

HS.model <- "
  visual  =~ x1 + x2 + x3
  textual =~ x4 + x5 + x6
  speed   =~ x7 + x8 + x9
"
utils::data("HolzingerSwineford1939", package = "lavaan")
fit <- acfa(HS.model, HolzingerSwineford1939, std.lv = TRUE, nsamp = 100,
            test = "standard", verbose = FALSE)

# Marginal log-likelihood (log evidence)
logLik(fit)

# Classical log-likelihood at the posterior mean, AIC/BIC-compatible
ll <- logLik(fit, type = "plugin")
AIC(ll)

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

Description

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,
  Sigma = 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,
  Sigma = NULL,
  cores = NULL,
  verbose = FALSE,
  ...
)

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

Arguments

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^*, \Sigma), 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' \Sigma s_u,

\log \mathrm{CPO}_u^{(2)} = \ell_u - \tfrac{1}{2} s_u' (\Sigma^{-1} + H_u)^{-1} s_u + \tfrac{1}{2} \log |I + \Sigma 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}. The effective number of parameters is 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. Units whose second-order curvature matrix is not positive definite fall back to first order for that unit only (flagged in per_unit$ok).

The type is resolved automatically: per-cluster ("loco") when the model was fitted with a cluster argument, per-subject ("loso") otherwise. For a two-level model these are the two estimands of Merkle, Furr & Rabe-Hesketh (2019): the default per-cluster "loco" is the marginal predictive (leave-one-cluster-out – prediction for a new cluster), while type = "loso" forces the conditional predictive (leave-one-unit-out – prediction for a new observation within an observed cluster), where row i of cluster j contributes \ell_i = \ell_j(\mathrm{full}) - \ell_j(\mathrm{minus\ row\ } i), the conditional density of the row given the rest of its cluster. The two answer different questions and are easily conflated, so the per-cluster marginal is the default and type = "loso" warns. It works with and without missing data, costs one cluster evaluation per row per Hessian direction, and is best subset with units.

Multigroup models. Groups are independent, so each unit is scored against its own group's implied moments; without a mean structure the exchangeability transformation applies per group, and cross-group equality constraints (group.equal) flow through the packed parameter space automatically. The per-unit results are stacked by group (a group column records the membership), and units are identified by case number – the row number of the analysed dataset – so a unit keeps its identity across fits that assign or order groups differently (e.g. a pooled fit versus a grouped fit of the same data, which compare() pairs unit by unit). This makes compare(..., loo = TRUE) the instrument of choice for the measurement-invariance ladder: configural, metric, and scalar fits are compared on a proper predictive scale with paired standard errors.

Supplying theta and/or Sigma scores the model at an arbitrary Gaussian posterior summary instead of the fit's own, without refitting. This is the building block for refit-free model exploration: for example, conditioning the encompassing model's summary on a parameter being zero (a rank-one update of theta and Sigma) and scoring the result gives the LOO of that submodel from a single fit. INLAvaan provides only this evaluation API; search strategies are left to the user. A conditioned Sigma may be singular; the computation automatically restricts to the non-degenerate block, which is exact.

Parallelism is strictly opt-in: the default cores = NULL runs serially, and cores > 1 parallelises the Hessian stage via forking (not available on Windows).

Calling loo() never modifies the fitted object. Under the default test = "standard", inlavaan() already computes and stores the full LOO at fit time whenever the model is supported, has a mean structure, and the predicted serial cost is within a 10-second budget (measured by timing one score evaluation); test = "loo" forces the computation regardless of the budget, and fit <- add_loo(fit) stores it post hoc. A stored result is returned directly by loo(fit) when called with default arguments, and is reused by fitmeasures() and compare() without recomputation.

Exogenous covariates. The flavour of the score follows the fitted likelihood. Under fixed.x = FALSE the covariates are modelled jointly and each unit is scored by the joint predictive density of its outcomes and covariates (flavour = "joint"). Under fixed.x = TRUE (the lavaan default) the fitted likelihood is the conditional one, and each unit is scored by the predictive density of its outcomes given its covariates (flavour = "conditional"); since the conditional likelihood is exactly invariant to the frozen covariate moments, this involves no additional approximation. The two flavours estimate different quantities whose scales differ by the covariate predictive density, so a joint and a conditional elpd must never be compared (compare() refuses mixed-flavour comparisons). Conditional scores of models conditioning on different covariate sets are comparable provided the outcome variables match – the natural setting for covariate selection. Both flavours support any covariate placement: single-level covariates, and cluster-level (between) and/or within-level covariates in two-level models.

Missing data. Fits estimated by full-information maximum likelihood (missing = "ml") are scored on the observed-data predictive: each unit contributes the density of the entries it actually has, with its full row (single-level) or whole cluster (two-level) removed from the conditioning set. For single-level fits the casewise kernels operate on each unit's observed subset, grouping rows by missing pattern, so a unit with fewer observed entries contributes a smaller log-likelihood term and a smaller score and thus self-weights in the elpd. Two-level fits are scored per cluster ("loco"): each cluster contributes its observed-data marginal likelihood, evaluated by lavaan's raw-data cluster kernels (no per-cluster sufficient statistics are needed, since leave-one-cluster-out deletes the whole cluster). All carry the same missing-at-random assumption as the FIML fit itself. Because the score is the observed-entry predictive, a compare() of two missing-data fits is meaningful only when they share the same observed entries (the same data and the same holes). The two-level conditional predictive (type = "loso") is available under missing data too, on the same kernels.

Supported models: continuous-indicator models fitted with the ML estimator (including FIML, missing = "ml", single- and two-level), single-group or multigroup (multigroup two-level models are not supported yet). If the loo package is attached it masks this generic, but loo(fit) continues to dispatch correctly because the method is registered by generic name.

Value

An object of class inlavaan_loo: a list with elements

per_unit

Data frame of pointwise results: unit (case number for LOSO, cluster position for LOCO), group (multigroup fits only), nobs (1 for LOSO, the cluster size for LOCO), l_star (unit log-likelihood at the summary), score_norm, lpd_1/lpd_2 (pointwise log predictive density), log_cpo_1/log_cpo_2 (pointwise LOO contributions), det_term, and ok (second-order success flag).

estimates

Matrix with rows elpd_loo, p_loo, looic and columns Estimate, SE (headline second-order values).

elpd_1, elpd_2, se_1, se_2, p_loo_1, p_loo_2

First- and second-order aggregates.

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

Metadata; flavour records whether units were scored jointly with their covariates ("joint") or conditionally on them ("conditional", for fixed.x fits).

add_loo() returns a copy of object with the LOO result 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 result always matches loo(fit).

See Also

fitmeasures(), compare(), inlavaan()

Examples

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)

# Leave-one-subject-out (LOSO) from the single fit -- no refitting
res <- loo(fit)
res
head(res$per_unit)

# 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
Sigma <- int$Sigma_theta
p <- which(names(coef(fit)) == "visual~~speed")
theta_c <- theta - Sigma[, p] * (theta[p] / Sigma[p, p])
Sigma_c <- Sigma - tcrossprod(Sigma[, p]) / Sigma[p, p]
loo(fit, theta = theta_c, Sigma = Sigma_c)

# 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)
loo(fit2l)

Plot an INLAvaan Object

Description

Generates diagnostic plots for a fitted INLAvaan model.

Usage

## S4 method for signature 'INLAvaan,ANY'
plot(
  x,
  y,
  type = c("marg_pdf", "sn_fit", "sn_fit_log"),
  params = "all",
  nrow = NULL,
  ncol = NULL,
  use_ggplot = TRUE,
  points = FALSE,
  ...
)

Arguments

See Also

diagnostics(), summary()

Examples

HS.model <- "
  visual  =~ x1 + x2 + x3
  textual =~ x4 + x5 + x6
  speed   =~ x7 + x8 + x9
"
utils::data("HolzingerSwineford1939", package = "lavaan")
fit <- acfa(HS.model, HolzingerSwineford1939, std.lv = TRUE, nsamp = 100,
            test = "none", verbose = FALSE)

# Posterior marginal densities (default)
plot(fit)

# Skew-normal fit diagnostic for a single parameter
plot(fit, type = "sn_fit", params = "visual=~x1")

Posterior Predictions for INLAvaan Models

Description

Compute posterior predictions from a fitted INLAvaan model, including latent variable scores, predicted observed values, and imputed missing data.

Usage

## S4 method for signature 'INLAvaan'
predict(
  object,
  type = c("lv", "yhat", "ov", "ypred", "ydist", "ymis", "ovmis"),
  newdata = NULL,
  level = 1L,
  nsamp = 1000,
  ymis_only = FALSE,
  summary = FALSE,
  ...
)

Arguments

Value

A list of nsamp posterior draws, each a matrix (or data frame, for multiple groups) with rows corresponding to cases and columns to variables or latent factors. When summary = TRUE, instead returns a summary.predict.inlavaan_internal object with the posterior mean, SD, quantiles, and mode for each case/variable.

See Also

sampling(), simulate(), summary()

Examples

HS.model <- "
  visual  =~ x1 + x2 + x3
  textual =~ x4 + x5 + x6
  speed   =~ x7 + x8 + x9
"
utils::data("HolzingerSwineford1939", package = "lavaan")
fit <- acfa(HS.model, HolzingerSwineford1939, std.lv = TRUE, nsamp = 100,
            test = "none", verbose = FALSE)

# Posterior latent variable scores
lv_scores <- predict(fit)
head(lv_scores)

# Predicted observed variable means
yhat <- predict(fit, type = "yhat")
head(yhat)

# Point estimates only, skipping the manual summary() step
predict(fit, type = "yhat", summary = TRUE)

Specify priors for a SEM

Description

Specify priors for a SEM, similar to how blavaan::dpriors() works.

Usage

priors_for(...)

Arguments

Details

This function provides a convenient way to specify prior distributions for different types of parameters in a structural equation model (SEM). It uses a registry of default priors for common lavaan parameter types (e.g., loadings, regressions, residuals, etc.) and allows users to override these defaults by passing named arguments.

The parameter names, and default settings, are:

• nu = "normal(0,32)": Observed variable intercepts

• alpha = "normal(0,10)": Latent variable intercepts

• lambda = "normal(0,10)": Factor loadings

• beta = "normal(0,10)": Regression coefficients

• theta = "gamma(1,.5)[sd]": Residual precisions

• psi = "gamma(1,.5)[sd]": Latent variable precisions

• rho = "beta(1,1)": Correlations (both latent and observed)

• tau = "normal(0,1.5)": Thresholds for ordinal variables

Note that the normal distributions are parameterised using standard deviations, and not variances. For example, normal(0,10) means a normal distribution with mean 0 and standard deviation 10 (not variance 10).

Value

A named character vector of prior specifications, where names correspond to lavaan parameter types (e.g., "lambda", "beta", "theta", etc.) and values are character strings specifying the prior distribution (e.g., "normal(0,10)", "gamma(1,0.5)[sd]", "gamma(1,1)[prec]", etc.).

Scale qualifiers

For variance parameters (theta, psi), the prior distribution can be placed on a transformed scale by appending a qualifier:

• [sd]: Prior is on the standard deviation \sigma. Example: "gamma(1,0.5)[sd]" places a Gamma(1, 0.5) prior on \sigma = \sqrt{\text{variance}}.

• [prec]: Prior is on the precision \tau = 1/\sigma^2. Example: "gamma(1,1)[prec]" places a Gamma(1, 1) prior on \tau = 1/\text{variance}. This is the parameterisation used by blavaan and corresponds to an Inverse-Gamma prior on the variance.

The necessary Jacobian adjustment is applied automatically in both cases.

See Also

inlavaan(), acfa(), asem(), agrowth()

Examples

priors_for(nu = "normal(0,10)", lambda = "normal(0,1)", rho = "beta(3,3)")

# Precision-scale prior for residual variances (blavaan-style)
priors_for(theta = "gamma(1,1)[prec]")

Fast Approximation of Skew-Normal Quantile Function

Description

A fast approximation of skew-normal quantiles using the high-performance approximation algorithm from the INLA GMRFLib C source, and originally by Thomas Luu (see details for reference).

Usage

qsnorm_fast(p, xi = 0, omega = 1, alpha = 0)

Arguments

Details

This function implements a high-performance approximation for the skew-normal quantile function based on the algorithm described by Luu (2016). The method uses a domain decomposition strategy to achieve high accuracy (< 10^{-7} relative error) without iterative numerical inversion.

The domain is split into two regions:

• Tail Regions: For extreme probabilities where \vert u \vert is large, the quantile is approximated using the Lambert W-function, W(z), solving z = \Phi(q) via asymptotic expansion:

  q \approx \sqrt{2 W\left(\frac{1}{2\pi (1-p)^2}\right)}

• Central Region: For the main body of the distribution, the function uses a high-order Taylor expansion of the inverse error function around a carefully selected expansion point $x_0$:

  \Phi^{-1}(p) \approx \sum_{k=0}^5 c_k (z - x_0)^k

This approach is significantly faster than standard numerical inversion (e.g., uniroot) while maintaining sufficient precision for most statistical applications.

Value

Vector of quantiles.

References

Luu, T. (2016). Fast and accurate parallel computation of quantile functions for random number generation (Doctoral thesis). UCL (University College London). https://discovery.ucl.ac.uk/1482128/

Examples

qsnorm_fast(c(0.025, 0.5, 0.975))
qsnorm_fast(c(0.025, 0.5, 0.975), xi = 2, omega = 0.5, alpha = 1)

Model Residuals for INLAvaan Models

Description

Extract the difference between the observed and model-implied (fitted) sample statistics from a fitted INLAvaan model. As in lavaan and blavaan, residuals are computed at the parameter estimates – here the posterior means – not as a posterior distribution over residuals.

Usage

## S4 method for signature 'INLAvaan'
residuals(object, type = "raw", labels = TRUE, ...)

## S4 method for signature 'INLAvaan'
resid(object, type = "raw", ...)

Arguments

Details

This delegates to lavaan's own residuals() machinery, so the return structure matches lavaan exactly. Because INLAvaan stores the posterior means as the point estimates of the fitted object, the residuals are the observed statistics minus the posterior-mean model-implied statistics (mirroring blavaan, which likewise inherits lavaan's residuals() without overriding it).

Value

For moment-based types, a list with elements type, cov, and (when relevant) mean. For type = "casewise", a numeric matrix of observed-minus-fitted values.

See Also

fitted(), predict(), fitMeasures()

Examples

HS.model <- "
  visual  =~ x1 + x2 + x3
  textual =~ x4 + x5 + x6
  speed   =~ x7 + x8 + x9
"
utils::data("HolzingerSwineford1939", package = "lavaan")
fit <- acfa(HS.model, HolzingerSwineford1939, std.lv = TRUE, nsamp = 100,
            test = "none", verbose = FALSE)

# Raw residual covariance matrix (posterior means)
residuals(fit)

# SRMR-basis residuals
residuals(fit, type = "cor.bentler")

# Casewise observed-minus-fitted values
head(residuals(fit, type = "casewise"))

Draw Samples from the Generative Model

Description

Sample model parameters, latent variables, or observed variables from the generative model underlying a fitted INLAvaan model. By default, parameters are drawn from the posterior distribution; set prior = TRUE to draw from the prior instead (useful for prior predictive checks).

Usage

sampling(object, ...)

## S4 method for signature 'INLAvaan'
sampling(
  object,
  type = c("lavaan", "theta", "latent", "observed", "implied", "all"),
  nsamp = 1000L,
  samp_copula = TRUE,
  prior = FALSE,
  silent = FALSE,
  ...
)

Arguments

Details

Each row of the output corresponds to a fresh parameter draw: a new \boldsymbol\theta^{(s)} is sampled and then propagated through the generative chain to produce one latent vector and one observed vector. This makes sampling() ideal for prior and posterior predictive checks (e.g., density overlays, test statistic distributions).

The generative chain is:

\boldsymbol\theta^{(s)} \sim \pi(\boldsymbol\theta \mid \mathbf{y})

\boldsymbol\eta^{(s)} \sim N((\mathbf{I} - \mathbf{B})^{-1}\boldsymbol\alpha,\,\boldsymbol\Phi)

\mathbf{y}^{*(s)} \sim N(\boldsymbol\Lambda\boldsymbol\eta^{(s)} + \boldsymbol\nu,\,\boldsymbol\Theta)

If you need complete replicate datasets (many observations from a single parameter draw) — for example, for simulation-based calibration (SBC) — use simulate() instead.

This is distinct from predict(), which computes individual-specific factor scores \boldsymbol\eta \mid \mathbf{y},\boldsymbol\theta conditional on observed data.

Value

A matrix or named list, depending on type.

See Also

simulate() for generating complete replicate datasets (e.g., for SBC); predict() for individual-specific factor scores; bfit_indices() for Bayesian fit indices.

Examples

utils::data("HolzingerSwineford1939", package = "lavaan")
fit <- acfa("visual =~ x1 + x2 + x3", HolzingerSwineford1939)

# Posterior samples of lavaan-side parameters
samps <- sampling(fit, nsamp = 500)
head(samps)

# Compare copula vs Gaussian sampling
s_cop <- sampling(fit, nsamp = 500, samp_copula = TRUE)
s_gaus <- sampling(fit, nsamp = 500, samp_copula = FALSE)

# Prior predictive samples
y_prior <- sampling(fit, type = "observed", nsamp = 500, prior = TRUE)

Simulate Datasets from the Generative Model

Description

Generate complete synthetic datasets from a fitted INLAvaan model. For each simulation, a single parameter vector is drawn (from the posterior or prior), and then sample.nobs observations are generated from the model-implied distribution at that parameter value.

Usage

## S4 method for signature 'INLAvaan'
simulate(
  object,
  nsim = 1L,
  seed = NULL,
  sample.nobs = NULL,
  prior = FALSE,
  samp_copula = TRUE,
  silent = FALSE,
  ...
)

Arguments

Details

This function is designed for tasks that require full replicate datasets from a single parameter draw, such as simulation-based calibration (SBC) and posterior predictive p-values. It differs from sampling() which generates one observation per parameter draw (useful for prior/posterior predictive density overlays).

For each simulation s = 1, \ldots, S:

1. Draw \boldsymbol\theta^{(s)} from the posterior (or prior).

2. Compute the model-implied covariance \boldsymbol\Sigma(\boldsymbol\theta^{(s)}). If it is not positive-definite, reject and redraw.

3. Generate a dataset of sample.nobs rows from N(\boldsymbol\mu(\boldsymbol\theta^{(s)}),\, \boldsymbol\Sigma(\boldsymbol\theta^{(s)})).

Parameter draws reuse the same internal machinery as sampling() (sample_params_prior / sample_params_posterior), so the prior specification is consistent.

Value

A list of length nsim. Each element is a data frame with sample.nobs rows and two attributes:

• "truth" — named numeric vector of lavaan-side (x-space, constrained) parameter values used to generate the dataset.

• "truth_theta" — named numeric vector of the corresponding unconstrained (theta-space) parameter values.

See Also

sampling() for single-observation draws from the predictive distribution (prior/posterior predictive checks).

Examples

utils::data("HolzingerSwineford1939", package = "lavaan")
fit <- acfa("visual =~ x1 + x2 + x3", HolzingerSwineford1939)

# Simulate one replicate dataset from the posterior
sims <- simulate(fit, nsim = 1)
head(sims[[1]])                    # data frame
attr(sims[[1]], "truth")           # true lavaan-side (x-space) parameters
attr(sims[[1]], "truth_theta")     # corresponding unconstrained (theta-space) parameters

# Simulate from the prior (e.g., for SBC)
sims_prior <- simulate(fit, nsim = 5, prior = TRUE)
lapply(sims_prior, nrow)

Standardised solution of a latent variable model

Description

Standardised solution of a latent variable model

Usage

standardisedsolution(
  object,
  type = "std.all",
  se = TRUE,
  ci = TRUE,
  level = 0.95,
  postmedian = FALSE,
  postmode = FALSE,
  cov.std = TRUE,
  remove.eq = TRUE,
  remove.ineq = TRUE,
  remove.def = FALSE,
  nsamp = 250,
  ...
)

standardisedSolution(
  object,
  type = "std.all",
  se = TRUE,
  ci = TRUE,
  level = 0.95,
  postmedian = FALSE,
  postmode = FALSE,
  cov.std = TRUE,
  remove.eq = TRUE,
  remove.ineq = TRUE,
  remove.def = FALSE,
  nsamp = 250,
  ...
)

standardizedsolution(
  object,
  type = "std.all",
  se = TRUE,
  ci = TRUE,
  level = 0.95,
  postmedian = FALSE,
  postmode = FALSE,
  cov.std = TRUE,
  remove.eq = TRUE,
  remove.ineq = TRUE,
  remove.def = FALSE,
  nsamp = 250,
  ...
)

standardizedSolution(
  object,
  type = "std.all",
  se = TRUE,
  ci = TRUE,
  level = 0.95,
  postmedian = FALSE,
  postmode = FALSE,
  cov.std = TRUE,
  remove.eq = TRUE,
  remove.ineq = TRUE,
  remove.def = FALSE,
  nsamp = 250,
  ...
)

Arguments

Value

A data.frame containing standardised model parameters.

See Also

summary(), coef(), vcov()

Examples

HS.model <- "
  visual  =~ x1 + x2 + x3
  textual =~ x4 + x5 + x6
"
utils::data("HolzingerSwineford1939", package = "lavaan")

# Fit a CFA model with standardised latent variables
fit <- acfa(
  HS.model,
  data = HolzingerSwineford1939,
  test = "none",
  nsamp = 10,
  vb_correction = FALSE,
  verbose = FALSE
)
standardisedsolution(fit, nsamp = 10, se = FALSE, ci = FALSE)

Timing Information for INLAvaan Models

Description

Extract wall-clock timings for individual computation stages of a fitted INLAvaan model.

Usage

timing(object, ...)

## S4 method for signature 'INLAvaan'
timing(object, what = "total", ...)

Arguments

Value

A named numeric vector (class c("timing.INLAvaan", "numeric")) of elapsed times in seconds. Printing formats short durations as seconds, longer ones as minutes or hours.

See Also

diagnostics(), summary()

Examples

HS.model <- "
  visual  =~ x1 + x2 + x3
  textual =~ x4 + x5 + x6
  speed   =~ x7 + x8 + x9
"
utils::data("HolzingerSwineford1939", package = "lavaan")
fit <- acfa(HS.model, HolzingerSwineford1939, std.lv = TRUE, nsamp = 100,
            test = "none", verbose = FALSE)

# Total elapsed time
timing(fit)

# All stages
timing(fit, what = "all")

# Specific stages
timing(fit, what = c("optim", "marginals"))

Update and re-fit an INLAvaan model

Description

Re-fit an INLAvaan model with modified arguments, in the style of stats::update(). This is convenient for prior-sensitivity analyses (vary dp or the prior() modifiers), iterative respecification (swap the model or add syntax), and swapping approximation settings (e.g. marginal_method, nsamp) without retyping the whole call.

Usage

## S4 method for signature 'INLAvaan'
update(object, model, add, ..., evaluate = TRUE)

Arguments

Details

The method edits the recorded inlavaan() call, replacing model (or extending it via add) and overriding any argument supplied through ..., then re-evaluates it. Because INLAvaan extends the lavaan class, a dedicated method is required: without it, update() would dispatch to lavaan's method and silently return a frequentist lavaan fit instead of a Bayesian one.

When the parameter structure is preserved (no change to model, add, or model.type), the previous posterior mode is passed as the optimiser's starting value (a warm start), so mode-finding is typically much faster than a cold fit. This makes prior sweeps particularly cheap. The warm start only changes the optimiser's starting point, not the fitted result. Supply start explicitly to override it.

Fits produced by INLAvaan versions that did not record their call cannot be updated; re-fit once with inlavaan() (or acfa(), asem(), agrowth()) and update() will work on the new object.

Value

A re-fitted INLAvaan object, or the updated call when evaluate = FALSE.

See Also

inlavaan(), compare()

Examples

## Not run: 
model <- "visual =~ x1 + x2 + x3"
fit <- acfa(model, data = lavaan::HolzingerSwineford1939)

# Prior sensitivity: refit under a tighter loading prior (warm-started)
fit2 <- update(fit, dp = priors_for(lambda = "normal(0,1)"))

# Respecify: add a residual covariance
fit3 <- update(fit, add = "x1 ~~ x2")
compare(fit, fit3)

# Swap approximation settings only
fit4 <- update(fit, marginal_method = "sampling", nsamp = 4000)

## End(Not run)

Variance-Covariance Matrix for INLAvaan Models

Description

Extract the posterior variance-covariance matrix of model parameters from a fitted INLAvaan model.

Usage

## S4 method for signature 'INLAvaan'
vcov(object, type = c("lavaan", "theta"), ...)

Arguments

Value

A square numeric matrix.

See Also

summary(), coef(), standardisedsolution()

Examples

HS.model <- "
  visual  =~ x1 + x2 + x3
  textual =~ x4 + x5 + x6
  speed   =~ x7 + x8 + x9
"
utils::data("HolzingerSwineford1939", package = "lavaan")
fit <- acfa(HS.model, HolzingerSwineford1939, std.lv = TRUE, nsamp = 100,
            test = "none", verbose = FALSE)

# Default: posterior covariance of lavaan parameters
vcov(fit)

# Internal parameterisation (Laplace approximation)
vcov(fit, type = "theta")

Widely Applicable Information Criterion for INLAvaan Models

Description

Computes the WAIC of a fitted INLAvaan model from unit log-likelihoods evaluated over posterior draws. Single-level models are scored per subject; two-level models are scored per cluster by default, matching the units used by loo(). For a two-level model type = "loso" instead scores the conditional (leave-one-unit-out) WAIC; see Details.

Usage

waic(x, ...)

## S3 method for class 'INLAvaan'
waic(
  x,
  type = c("auto", "loso", "loco"),
  units = NULL,
  nsamp = NULL,
  cores = NULL,
  verbose = FALSE,
  ...
)

## S3 method for class 'inlavaan_internal'
waic(
  x,
  type = c("auto", "loso", "loco"),
  units = NULL,
  nsamp = NULL,
  cores = NULL,
  verbose = FALSE,
  ...
)

Arguments

Details

For each posterior draw \theta_s (drawn with the same copula sampler used for the fit's posterior summaries) and unit u, the log-likelihood \log p(y_u \mid \theta_s) is evaluated; then \mathrm{lpd}_u = \log \tfrac{1}{S}\sum_s p(y_u \mid \theta_s), p_{\mathrm{waic},u} = \mathrm{var}_s \log p(y_u \mid \theta_s), and \mathrm{elpd}_{\mathrm{waic}} = \sum_u (\mathrm{lpd}_u - p_{\mathrm{waic},u}) with \mathrm{WAIC} = -2\, \mathrm{elpd}_{\mathrm{waic}}. Unlike loo(), this is a sampling-based estimate: results vary with the random draws, and units with p_{\mathrm{waic},u} > 0.4 trigger a reliability warning (also annotated when printing). The same model restrictions as loo() apply, and so does the flavour rule: fits with fixed.x = TRUE are scored conditionally on the exogenous covariates, fits with fixed.x = FALSE jointly (see loo()).

Marginal vs conditional WAIC (two-level models). The default per-cluster scoring is the marginal WAIC, which corresponds to leave-one-cluster-out cross-validation – prediction for a new cluster. Setting type = "loso" scores the conditional WAIC, corresponding to leave-one-unit-out – prediction for a new observation within an observed cluster (each row contributes the conditional density of its observed entries given the rest of its cluster). The two answer different questions and are easily conflated (Merkle, Furr & Rabe-Hesketh, 2019); the per-cluster marginal is the usual model-comparison target, so it is the default, and type = "loso" warns. This matches loo(type = "loso") – the two compute the same estimand by sampling and by Taylor expansion.

Under the default test = "standard", inlavaan() computes the WAIC at fit time by reusing the posterior draws the fit already produced (when the model is supported and nsamp >= 100), and stores it with the fit: waic(fit) then returns the stored result when called with default arguments, and fitmeasures() reports waic, p_waic, se_waic as part of "all" for free. If the loo package is attached it masks this generic, but dispatch on INLAvaan objects continues to work.

Value

An object of class inlavaan_waic: a list with per_unit (pointwise lpd, p_waic, elpd_waic, with the same unit/group identification as loo()), estimates (matrix with rows elpd_waic, p_waic, waic and columns Estimate, SE), type, n_units, n_groups, and nsamp.

See Also

loo(), fitmeasures()

Examples

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)
waic(fit)