Nothing
R/lav_sam_utils.R
In lavaan: Latent Variable Analysis
Defines functions
lav_sam_veta_con
lav_sam_veta_pt
lav_sam_mmlist_merge_blocks
lav_sam_condx_perm_idx
lav_sam_wls_obs_labels
lav_sam_get_mmlist
lav_sam_get_cov_ybar
lav_sam_table
lav_sam_trunc_bias_se
lav_sam_trunc_bias
lav_sam_global_test
lav_sam_struc_fit
lav_sam_struc_fit_object
lav_sam_step3_joint
lav_sam_fs2
lav_sam_eeta2
lav_sam_veta2
lav_sam_veta2_floor_debias
lav_sam_fs_missing
lav_sam_eeta
lav_sam_veta
lav_sam_veta1_floor_debias
lav_sam_lambda_floor_debias
lav_sam_tmat
lav_sam_mapping_mat_tmat
lav_sam_mapping_mat
lav_sam_struc_object
lav_sam_local_flag
lav_sam_lvnames2
lav_sam_alpha_n1
# utility functions for the sam() function
# YR 4 April 2023
# clamped small-sample correction factor alpha/(n - 1), forced into [0, 1]
lav_sam_alpha_n1 <- function(alpha_correction = 0L, n = 20L) {
alpha_n1 <- alpha_correction / (n - 1)
if (!is.finite(alpha_n1)) {
alpha_n1 <- 0
}
min(max(alpha_n1, 0), 1)
}
# row/column index pairs and labels for the elements of (eta* %x% eta*),
# where lv_names already includes the leading "..int.." intercept entry;
# the first nfac labels are the plain (first-order) lv names
lav_sam_lvnames2 <- function(lv_names) {
nfac <- length(lv_names)
idx1 <- rep(seq_len(nfac), each = nfac)
idx2 <- rep(seq_len(nfac), times = nfac)
names2 <- paste(lv_names[idx1], lv_names[idx2], sep = ":")
names2[seq_len(nfac)] <- lv_names
list(idx1 = idx1, idx2 = idx2, names = names2)
}
# is this a sam object created with sam.method = "local", "fsr" or "cfsr"?
lav_sam_local_flag <- function(object) {
sam_method <- object@internal$sam.method
!is.null(sam_method) && sam_method %in% c("local", "fsr", "cfsr")
}
# the stored step-2 structural fit of a local SAM object (with the two-step
# corrected test/baseline injected, if available); NULL for global sam
# objects, and for sam objects created by older versions of lavaan
lav_sam_struc_object <- function(object) {
object@internal$sam.struc.fit.object
}
# construct 'mapping matrix' M using either "ML", "GLS" or "ULS" method
#
# by construction, M %*% LAMBDA = I (the identity matrix)
lav_sam_mapping_mat <- function(mm_lambda = NULL, mm_theta = NULL,
s = NULL, s_inv = NULL,
method = "ML") {
method <- toupper(method)
# catch empty columns in LAMBDA (eg higher order!)
lambda_orig <- mm_lambda
empty_idx <- which(apply(mm_lambda, 2L, function(x) all(x == 0)))
if (length(empty_idx) > 0L) {
mm_lambda <- lambda_orig[, -empty_idx, drop = FALSE]
}
# ULS
# M == solve( t(LAMBDA) %*% LAMBDA ) %*% t(LAMBDA)
# == MASS:::ginv(LAMBDA)
if (method == "ULS") {
m <- try(tcrossprod(solve(crossprod(mm_lambda)), mm_lambda),
silent = TRUE
)
if (inherits(m, "try-error")) {
lav_msg_warn(gettext(
"cannot invert crossprod(LAMBDA); using generalized inverse"))
m <- MASS::ginv(mm_lambda)
}
# GLS
# M == solve( t(LAMBDA) %*% S.inv %*% LAMBDA ) %*% t(LAMBDA) %*% S.inv
} else if (method == "GLS") {
if (is.null(s_inv)) {
s_inv <- try(solve(s), silent = TRUE)
}
if (inherits(s_inv, "try-error")) {
lav_msg_warn(gettext("S is not invertible; switching to ULS method"))
m <- lav_sam_mapping_mat(mm_lambda = mm_lambda, method = "ULS")
} else {
t_lsinv <- t(mm_lambda) %*% s_inv
t_lsinv_l <- t_lsinv %*% mm_lambda
m <- try(solve(t_lsinv_l, t_lsinv), silent = TRUE)
if (inherits(m, "try-error")) {
lav_msg_warn(gettext("problem constructing mapping matrix;
switching to generalized inverse"))
m <- MASS::ginv(t_lsinv_l) %*% t_lsinv
}
}
# ML
# M == solve(t(LAMBDA) %*% THETA.inv %*% LAMBDA) %*% t(LAMBDA) %*% THETA.inv
} else if (method == "ML") {
# Problem: if THETA has zero elements on the diagonal, we cannot invert
# As we do not have access to Sigma(.inv), we cannot
# use the trick as in lavPredict() (where we replace THETA.inv by
# Sigma.inv)
# old method (<0.6-16): remove rows/cols with zero values
# on the diagonal of THETA, invert the submatrix and
# set the '0' diagonal elements to one. This resulted in (somewhat)
# distorted results.
# new in 0.6-16: we use the Wall & Amemiya (2000) method using
# the so-called 'T' transformation
# new in 0.6-18: if we cannot use the marker method (eg growth models),
# use the 'old' method anyway: remove zero rows/cols and invert submatrix
zero_theta_idx <- which(abs(diag(mm_theta)) < 1e-4) # be conservative
if (length(zero_theta_idx) == 0L) {
# ok, no zero diagonal elements: try to invert THETA
if (lav_mat_is_diagonal(mm_theta)) {
theta_inv <- diag(1 / diag(mm_theta), nrow = nrow(mm_theta))
} else {
theta_inv <- try(solve(mm_theta), silent = TRUE)
if (inherits(theta_inv, "try-error")) {
theta_inv <- NULL
}
}
} else {
# see if we can use marker method
marker_idx <- lav_utils_get_marker(mm_lambda = mm_lambda, std_lv = TRUE)
if (any(is.na(marker_idx))) {
theta_inv <- try(lav_mat_sym_inverse(mm_theta), silent = TRUE)
if (inherits(theta_inv, "try-error")) {
theta_inv <- NULL
} else {
diag(theta_inv)[zero_theta_idx] <- 1
}
} else {
# try tmat method later on
theta_inv <- NULL
}
}
# could we invert THETA?
if (!is.null(theta_inv)) {
# ha, all is good; compute M the usual way
t_lti <- t(mm_lambda) %*% theta_inv
t_lti_l <- t_lti %*% mm_lambda
m <- try(solve(t_lti_l, t_lti), silent = TRUE)
if (inherits(m, "try-error")) {
lav_msg_warn(gettext(
"problem constructing ML mapping matrix; switching to ULS"))
m <- lav_sam_mapping_mat(mm_lambda = mm_lambda, method = "ULS")
}
} else {
# use W&A2000's method using the 'T' transformation
m <- try(lav_sam_mapping_mat_tmat(
mm_lambda = mm_lambda,
mm_theta = mm_theta
), silent = TRUE)
if (inherits(m, "try-error")) {
lav_msg_warn(gettext(
"problem constructing ML mapping matrix; switching to ULS"))
m <- lav_sam_mapping_mat(mm_lambda = mm_lambda, method = "ULS")
}
}
} # ML
# empty.idx?
if (length(empty_idx) > 0L) {
m_full <- m
m <- matrix(0, nrow = ncol(lambda_orig), ncol = ncol(m_full))
m[-empty_idx, ] <- m_full
}
m
}
# use 'T' transformation to create the 'Bartlett/ML' mapping matrix
# see Wall & Amemiya (2000) eq (7)
# see also Fuller 1987 page 357 (where T is called H), and page 364
# although Fuller/W&A always assumed that THETA is diagonal,
# their method seems to work equally well for non-diagonal THETA
#
# in our implementation:
# - we do NOT reorder the rows of LAMBDA
# - if std.lv = TRUE, we first rescale to 'create' marker indicators
# and then rescale back at the end
#
lav_sam_mapping_mat_tmat <- function(mm_lambda = NULL,
mm_theta = NULL,
marker_idx = NULL,
std_lv = NULL) {
mm_lambda <- as.matrix.default(mm_lambda)
# catch empty columns in LAMBDA (eg higher order!)
lambda_orig <- mm_lambda
empty_idx <- which(apply(mm_lambda, 2L, function(x) all(x == 0)))
if (length(empty_idx) > 0L) {
mm_lambda <- lambda_orig[, -empty_idx, drop = FALSE]
}
# nvar <- nrow(mm_lambda)
# nfac <- ncol(mm_lambda)
# do we have marker.idx?
if (is.null(marker_idx)) {
# 'marker' indicator has a single non-zero element in a row
marker_idx <- lav_utils_get_marker(mm_lambda = mm_lambda, std_lv = TRUE)
if (any(is.na(marker_idx))) {
lav_msg_stop(gettext("no clear markers in LAMBDA matrix"))
}
}
# std.lv TRUE or FALSE?
if (is.null(std_lv)) {
std_lv <- FALSE
if (any(diag(mm_lambda[marker_idx, , drop = FALSE]) != 1)) {
std_lv <- TRUE
}
}
# if std.lv = TRUE, rescale
if (std_lv) {
marker <- mm_lambda[marker_idx, , drop = FALSE]
marker_inv <- 1 / diag(marker)
mm_lambda <- t(t(mm_lambda) * marker_inv)
}
# compute 'T' matrix
tmat <- lav_sam_tmat(
mm_lambda = mm_lambda, mm_theta = mm_theta,
marker_idx = marker_idx
)
# ML mapping matrix
m <- tmat[marker_idx, , drop = FALSE]
if (std_lv) {
m <- m * marker_inv
}
# empty.idx?
if (length(empty_idx) > 0L) {
m_full <- m
m <- matrix(0, nrow = ncol(lambda_orig), ncol = ncol(m_full))
m[-empty_idx, ] <- m_full
}
m
}
# create 'T' matrix (tmat) for T-transformation
#
# Notes: - here we assume that LAMBDA has unity markers (no std.lv = TRUE)
# - TMAT is NOT symmetric!
# - Yc %*% t(TMAT) transforms the data in such a way that we get:
# 1) Bartlett factor scores in the marker columns
# 2) 'V' values in the non-marker columns, where:
# V = Yc - Yc[,marker.idx] %*% t(LAMBDA)
#
lav_sam_tmat <- function(mm_lambda = NULL,
mm_theta = NULL,
marker_idx = NULL) {
mm_lambda <- as.matrix.default(mm_lambda)
nvar <- nrow(mm_lambda)
nfac <- ncol(mm_lambda)
# do we have marker.idx?
if (is.null(marker_idx)) {
# 'marker' indicator has a single 1 element in a row
marker_idx <- lav_utils_get_marker(mm_lambda = mm_lambda, std_lv = FALSE)
if (any(is.na(marker_idx))) {
lav_msg_stop(gettext("no clear markers in LAMBDA matrix"))
}
}
# construct 'm_c' matrix
c2 <- diag(nvar)
c2[, marker_idx] <- -1 * mm_lambda
m_c <- c2[-marker_idx, , drop = FALSE]
# compute Sigma.ve and Sigma.vv
sigma_ve <- m_c %*% mm_theta
# Sigma.vv <- m_c %*% THETA %*% t(m_c)
sigma_vv <- sigma_ve %*% t(m_c)
# construct 'Gamma' (=gamma_1) (and gamma2) matrix
# Gamma <- (t(Sigma.ve) %*% solve(Sigma.vv))[marker.idx,, drop = FALSE]
gamma_1 <- try(t(solve(sigma_vv, sigma_ve)[, marker_idx, drop = FALSE]),
silent = TRUE
)
if (inherits(gamma_1, "try-error")) {
tmp <- t(sigma_ve) %*% MASS::ginv(sigma_vv)
gamma_1 <- tmp[marker_idx, , drop = FALSE]
}
gamma2 <- matrix(0, nfac, nvar)
gamma2[, -marker_idx] <- gamma_1
gamma2[, marker_idx] <- diag(nfac)
# transformation matrix 'T' (we call it here 'Tmat')
tmat <- matrix(0, nvar, nvar)
tmat[-marker_idx, ] <- m_c
tmat[marker_idx, ] <- -gamma2 %*% c2
tmat
}
# Gated, bootstrap-debiased truncation floor for the lambda corrections --
# the shared core behind lav_sam_veta1_floor_debias() (first order, see
# lav_sam_veta()) and lav_sam_veta2_floor_debias() (second order, see
# lav_sam_veta2()).
#
# 'score' holds the casewise contributions whose (n-divisor) covariance is
# the left-hand side of the pencil: the Bartlett factor scores
# f_i = M (y_i - ybar) for the first order (Var(f) = MSM), and the
# second-order products fs2_i (kept columns) for the second order
# (Var(fs2) = var_fs2); 'err' is the right-hand side of the pencil (MTM
# resp. var.error), and lambda its smallest generalized root.
#
# The sample lambda is downward-biased -- mildly in well-separated
# settings, severely when the bottom of the spectrum is clustered (highly
# correlated factors, weak reliability, many product terms). The
# historical rules truncate with floors (1/(n-1) resp. p2/(n-1)) that
# vanish with n, so in difficult settings the accuracy of the structural
# estimates can deteriorate with increasing n. This rule is designed to
# be a strict no-op whenever the sample is comfortably away from the
# boundary, and to repair the difficult cases:
#
# 1. GATE: six deterministic interleaved half-samples give a cheap noise
# estimate se(lambda) = sd(lambda_half)/sqrt(2). If
# lambda - 1 >= max(2 * se(lambda), floor_default)
# the sample margin is comfortably positive: return NULL, and the
# caller applies the historical rule -- which then performs the full
# correction, bit-for-bit identical to the historical behavior.
# 2. Otherwise (the danger zone), estimate the bias and the sampling
# noise of lambda by a casewise bootstrap of the score rows (B = 200
# resamples; a fixed internal seed makes the resampling plan
# deterministic, and the user's .Random.seed is saved and restored):
# bias = max(0, mean(lambda_boot) - lambda)
# se_b = sd(lambda_boot)
# and return the floor
# max(min(lambda - 1 + bias, lambda - 1 + 3 * se_b),
# 2 * se_b, floor_default)
# i.e. the debiased margin, capped at three standard errors above the
# sample margin, and never below (twice) the noise level of lambda nor
# below the historical floor. The caller truncates iff
# lambda < 1 + floor, with lambda_star = max(0, lambda - floor).
#
# Return value: a numeric floor (danger zone only), or NULL. NULL covers
# both the open gate and any ineligible/degenerate situation (no usable
# scores, n < 4 * ncol(score), failed sub-sample roots); on NULL the
# caller falls back to the historical rule, which in the open-gate case
# applies the full correction (the gate implies lambda >= 1 + cutoff -
# 1), so the result is bit-for-bit the historical behavior.
lav_sam_lambda_floor_debias <- function(score = NULL, err = NULL,
lambda = NULL, n = NULL,
floor_default = NULL) {
if (is.null(score) || nrow(score) != n || anyNA(score) ||
n < 4L * ncol(score) || !is.finite(lambda)) {
return(NULL) # unusable: historical rule
}
idx <- seq_len(n)
mod2 <- idx %% 2L
mod4 <- idx %% 4L
half_list <- list(
which(mod2 == 1L), which(mod2 == 0L),
which(mod4 <= 1L), which(mod4 >= 2L),
which(mod4 == 0L | mod4 == 3L), which(mod4 == 1L | mod4 == 2L)
)
lambda_sub <- function(ii) {
nh <- length(ii)
vh <- var(score[ii, , drop = FALSE]) * (nh - 1) / nh
out <- try(lav_mat_sym_diff_smallest_root(vh, err), silent = TRUE)
if (inherits(out, "try-error")) {
return(as.numeric(NA))
}
out
}
lambda_half <- vapply(half_list, lambda_sub, numeric(1L))
lambda_half <- lambda_half[is.finite(lambda_half)]
if (length(lambda_half) < 4L) {
return(NULL) # degenerate: historical rule
}
se_half <- stats::sd(lambda_half) / sqrt(2)
# the gate: comfortably positive sample margin -> historical behavior
if (lambda - 1 >= max(2 * se_half, floor_default)) {
return(NULL)
}
# danger zone: casewise bootstrap for bias and noise of lambda
# (B = 200: with fewer resamples, the Monte Carlo error of se_boot -- and
# hence of the floor, which is a deterministic but arbitrary function of
# the internal seed -- is not negligible: with B = 40, lambda_star varied
# by about +/- 0.1 across seeds in small samples; the cost remains
# negligible, and the floor is frozen in the numerical jacobians)
b_reps <- 200L
seed_old <- NULL
if (exists(".Random.seed", envir = globalenv())) {
seed_old <- get(".Random.seed", envir = globalenv())
}
set.seed(272727L)
lambda_boot <- vapply(seq_len(b_reps), function(b) {
lambda_sub(sample.int(n, n, replace = TRUE))
}, numeric(1L))
if (!is.null(seed_old)) {
assign(".Random.seed", seed_old, envir = globalenv())
}
lambda_boot <- lambda_boot[is.finite(lambda_boot)]
if (length(lambda_boot) < b_reps / 2) {
return(NULL) # degenerate: historical rule
}
bias_hat <- max(0, mean(lambda_boot) - lambda)
se_boot <- stats::sd(lambda_boot)
margin <- min(lambda - 1 + bias_hat, lambda - 1 + 3 * se_boot)
out <- max(margin, 2 * se_boot, floor_default)
# branch bookkeeping for the SE machinery: when the lambda-tracking
# margin binds, lambda_star = lambda - floor is locally constant in
# lambda (kappa = 0); when a lambda-independent guard binds (2*se_boot
# or floor_default), lambda_star tracks lambda one-for-one (kappa = 1)
attr(out, "kappa") <- if (margin >= max(2 * se_boot, floor_default)) {
0
} else {
1
}
out
}
# First-order wrapper: Bartlett factor scores f_i = M (y_i - ybar), so
# that Var(f) equals MSM = M S M' (the per-subsample re-centering of
# var() makes the global centering immaterial). The caller only passes y
# in the supported setting (single-level, continuous, complete data, no
# conditional.x, no dummy/higher-order lvs); everything else -> NULL ->
# historical rule (see lav_sam_lambda_floor_debias() for the design).
lav_sam_veta1_floor_debias <- function(y = NULL, m = NULL, mtm = NULL,
lambda = NULL, n = NULL) {
if (is.null(y) || nrow(y) != n || anyNA(y)) {
return(NULL) # unusable: historical rule
}
yc <- t(t(y) - colMeans(y))
lav_sam_lambda_floor_debias(
score = yc %*% t(m), err = mtm, lambda = lambda, n = n,
floor_default = 1 / (n - 1)
)
}
# compute VETA
# - if alpha_correction == 0 -> same as local SAM (or MOC)
# - if alpha_correction == (N-1) -> same as FSR+Bartlett
lav_sam_veta <- function(m = NULL, s = NULL, mm_theta = NULL,
alpha_correction = 0L, lambda_correction = TRUE,
lambda_floor = "default", y = NULL,
n = 20L, dummy_lv_idx = integer(0L), extra = FALSE) {
# catch empty rows in M (higher-order?)
m_orig <- m
empty_idx <- which(apply(m_orig, 1L, function(x) all(x == 0)))
if (length(empty_idx) > 0L) {
m <- m_orig[-empty_idx, , drop = FALSE]
}
# MSM
msm <- m %*% s %*% t(m)
# MTM
mtm <- m %*% mm_theta %*% t(m)
# note: MTM no longer needs to be regularized (force-pd) here:
# lav_mat_sym_diff_smallest_root() handles a singular (or even zero)
# MTM natively, by regressing out its null space (which includes the
# rows/columns of any dummy lvs, and any empty theta elements)
if (all(abs(mtm) < sqrt(.Machine$double.eps))) {
# all zero? certainly no need for alpha or lambda_correction
alpha_correction <- 0L
lambda_correction <- FALSE
}
# apply small sample correction (if requested)
alpha <- alpha_correction
if (alpha_correction > 0) {
alpha_n1 <- lav_sam_alpha_n1(alpha_correction, n)
mtm <- (1 - alpha_n1) * mtm
}
lambda <- lambda_star <- +Inf
lambda_kappa <- 0
lambda_w <- NULL
lambda_floor_used <- as.numeric(NA)
if (lambda_correction && is.list(lambda_floor)) {
# internal (SE machinery only): frozen effective multiplier -- the
# numeric-jacobian channels re-run this function at perturbed inputs,
# and a kappa = 0 truncation (lambda-tracking debias margin) must keep
# lambda_star fixed at its at-the-estimates value there
lambda_star <- lambda_floor[["star"]]
veta <- msm - lambda_star * mtm
} else if (lambda_correction) {
# use Fuller (1987) approach to ensure VETA is positive
lambda <- try(lav_mat_sym_diff_smallest_root(msm, mtm),
silent = TRUE
)
if (inherits(lambda, "try-error")) {
lav_msg_warn(gettext("failed to compute lambda"))
veta <- msm - mtm # and hope for the best
} else {
# determine the truncation floor:
# NULL -> historical rule (cutoff 1 + 1/(n-1), floor 1/(n-1))
# a number -> coherent rule: truncate iff lambda < 1 + floor
floor_1 <- NULL
kappa_1 <- 1 # a fixed (numeric) floor: lambda_star tracks lambda
if (identical(lambda_floor, "debias") &&
length(empty_idx) == 0L && length(dummy_lv_idx) == 0L) {
# gated bootstrap-debiased floor; NULL when the sample margin is
# comfortably positive (then lambda >= 1 + 1/(n-1) by the gate
# construction, and the historical rule below applies the full
# correction: bit-for-bit the historical behavior) or when the
# setting is ineligible/degenerate (-> historical rule as well)
floor_1 <- lav_sam_veta1_floor_debias(
y = y, m = m, mtm = mtm, lambda = lambda, n = n
)
if (!is.null(floor_1) && !is.null(attr(floor_1, "kappa"))) {
kappa_1 <- attr(floor_1, "kappa")
}
} else if (is.numeric(lambda_floor)) {
floor_1 <- lambda_floor
}
if (!is.null(floor_1)) {
lambda_floor_used <- as.numeric(floor_1)
if (lambda < 1 + floor_1) {
lambda_star <- max(0, lambda - floor_1)
veta <- msm - lambda_star * mtm
# clamped at zero -> VETA = MSM exactly, no lambda term left
lambda_kappa <- if (lambda_star > 0) kappa_1 else 0
} else {
veta <- msm - mtm
}
} else {
# historical rule
lambda_floor_used <- 1 / (n - 1)
cutoff <- 1 + 1 / (n - 1)
if (lambda < cutoff) {
lambda_star <- lambda - 1 / (n - 1)
veta <- msm - lambda_star * mtm
lambda_kappa <- 1 # fixed floor, no clamp: always tracks lambda
} else {
veta <- msm - mtm
}
}
# eigenvector for the SE machinery: when the truncation binds with
# lambda_star tracking lambda (kappa = 1), the jacobian of
# vech(VETA) w.r.t. vech(S) carries the extra rank-1 term
# -vech(MTM) . (d lambda / d vech(S))' with
# d lambda = t(w) dS w, w = t(M) v (v the smallest-root pencil
# eigenvector, normalized t(v) MTM v = 1)
if (extra && lambda_kappa == 1) {
# recompute of the same pencil: suppress its (duplicate) warnings
lam_v <- try(suppressWarnings(
lav_mat_sym_diff_smallest_root(msm, mtm,
vector = TRUE)),
silent = TRUE)
vv <- attr(lam_v, "v")
if (!inherits(lam_v, "try-error") && !is.null(vv) &&
all(is.finite(vv))) {
lambda_w <- drop(t(m) %*% vv)
}
# no usable direction -> lambda_w stays NULL and the analytic
# rank-1 term is simply skipped (old, conservative behavior);
# kappa is kept, so the frozen re-runs still track lambda
}
}
} else {
veta <- msm - mtm
}
# empty.idx?
if (length(empty_idx) > 0L) {
msm_full <- msm
mtm_full <- mtm
veta_full <- veta
nfac_orig <- nrow(m_orig)
msm <- mtm <- veta <- matrix(0, nrow = nfac_orig, ncol = nfac_orig)
msm[-empty_idx, -empty_idx] <- msm_full
mtm[-empty_idx, -empty_idx] <- mtm_full
veta[-empty_idx, -empty_idx] <- veta_full
}
# extra attributes?
if (extra) {
attr(veta, "lambda") <- lambda
attr(veta, "alpha") <- alpha
attr(veta, "lambda_star") <- lambda_star
attr(veta, "MSM") <- msm
attr(veta, "MTM") <- mtm
# SE machinery (see lav_sam_jacb_lambda1_g() and the frozen-floor
# re-runs in the numeric jacobian channels)
attr(veta, "lambda_kappa") <- lambda_kappa
attr(veta, "lambda_w") <- lambda_w
attr(veta, "lambda_floor_used") <- lambda_floor_used
}
veta
}
# compute EETA = E(Eta) = M %*% [YBAR - NU]
lav_sam_eeta <- function(m = NULL, ybar = NULL, mm_nu = NULL) {
m %*% (ybar - mm_nu)
}
# Per-missing-pattern Bartlett factor scores and their measurement-error
# covariances, for the lv-interaction (casewise) machinery under
# missing = "ml"/"ml.x".
#
# For each missing pattern p (observed variables o), the pattern-specific
# mapping matrix M_p is built from the observed sub-matrices (Lambda_o,
# Theta_o), so that M_p Lambda_o = I and hence E(f_i | eta_i) = eta_i for
# every scoreable pattern; the measurement-error covariance of the scores,
# Var(f_i | eta_i) = M_p Theta_o M_p' =: B_p, now varies over the patterns.
# The second-order (VETA2/Gamma) machinery aggregates over the patterns; see
# lav_sam_veta2() for the exact aggregation (B.bar / Q / D).
#
# A pattern is 'scoreable' if every factor keeps at least one observed
# indicator (no all-zero columns in Lambda_o), Lambda_o has full column rank,
# and the mapping matrix could be constructed with no all-zero rows. Cases in
# unscoreable patterns get ok = FALSE and are excluded from the casewise
# machinery -- their information still enters the first-order moments through
# the EM (h1) statistics.
#
# y: raw data (rows = cases, columns = ov in model order), with NAs
# mm_lambda/mm_theta/mm_nu: measurement matrices WITHOUT the lv-interaction
# columns (and in the unstandardized metric; any std.lv rescaling is
# applied by the caller, to scores and B_p alike)
# s: covariance matrix of the ov (EM/h1 based; only used for method = "GLS")
# mp: optional precomputed missing-pattern structure (lav_data_mi_patterns(y))
#
# Returns list(fs = n x nfac score matrix (NA rows where ok = FALSE),
# ok = logical(n), pat_of_case = integer(n) index into bp,
# bp = list of B_p over the scoreable patterns,
# wp = pattern proportions n_p / sum(n_p) over scoreable cases)
lav_sam_fs_missing <- function(y = NULL, mm_lambda = NULL, mm_theta = NULL,
mm_nu = NULL, s = NULL, method = "ML",
mp = NULL) {
n <- nrow(y)
nfac <- ncol(mm_lambda)
if (is.null(mp)) {
mp <- lav_data_mi_patterns(y)
}
fs <- matrix(as.numeric(NA), nrow = n, ncol = nfac)
ok <- rep(FALSE, n)
pat_of_case <- rep(NA_integer_, n)
bp <- list()
wp_n <- integer(0L)
nu <- drop(mm_nu)
for (p in seq_len(mp$npatterns)) {
obs <- mp$pat[p, ]
lam_o <- mm_lambda[obs, , drop = FALSE]
# each factor needs at least one observed indicator, and the observed
# sub-LAMBDA must keep full column rank
if (any(apply(lam_o, 2L, function(x) all(x == 0))) ||
qr(lam_o)$rank < nfac) {
next
}
th_o <- mm_theta[obs, obs, drop = FALSE]
s_o <- NULL
if (!is.null(s)) {
s_o <- s[obs, obs, drop = FALSE]
}
m_p <- tryCatch(
lav_sam_mapping_mat(
mm_lambda = lam_o, mm_theta = th_o,
s = s_o, method = method
),
error = function(e) NULL
)
if (is.null(m_p) || anyNA(m_p) ||
any(apply(m_p, 1L, function(x) all(x == 0)))) {
next
}
case_idx <- mp$case.idx[[p]]
yc_o <- t(t(y[case_idx, obs, drop = FALSE]) - nu[obs])
this_id <- length(bp) + 1L
bp[[this_id]] <- m_p %*% th_o %*% t(m_p)
wp_n <- c(wp_n, length(case_idx))
fs[case_idx, ] <- yc_o %*% t(m_p)
ok[case_idx] <- TRUE
pat_of_case[case_idx] <- this_id
}
list(
fs = fs, ok = ok, pat_of_case = pat_of_case,
bp = bp, wp = wp_n / sum(wp_n)
)
}
# Second-order wrapper: the scores are the (kept columns of the)
# second-order factor-score products fs2, whose covariance is var_fs2;
# the error side of the pencil is var.error. With missing data the
# scores are pattern-based (heterogeneous error contributions), so we
# fall back to the historical rule (see lav_sam_lambda_floor_debias()
# for the design). Note that the second-order sample lambda is much more
# severely biased than its first-order counterpart: the pencil is
# high-dimensional (all latent variables plus all product terms) and
# var_fs2 is built from fourth-order moments of the factor scores.
lav_sam_veta2_floor_debias <- function(fs2 = NULL, var_error = NULL,
lambda = NULL, n = NULL,
mi_flag = FALSE) {
if (mi_flag) {
return(NULL) # pattern-based scores: historical rule
}
lav_sam_lambda_floor_debias(
score = fs2, err = var_error, lambda = lambda, n = n,
floor_default = ncol(fs2) / (n - 1)
)
}
# compute veta including quadratic/interaction terms
lav_sam_veta2 <- function(fs = NULL, m = NULL,
veta = NULL, eeta = NULL, mm_theta = NULL,
lv_names = NULL,
lv_int_names = NULL,
dummy_lv_names = character(0L),
alpha_correction = 0L,
lambda_correction = TRUE,
lambda_floor = "default",
lambda1 = 1,
pattern_list = NULL,
fs_outlier_idx = integer(0L),
return_fs = FALSE,
return_cov_iveta2 = TRUE,
extra = FALSE) {
# small utility function: var() divided by N
varn <- function(x, n) {
var(x, use = "pairwise.complete.obs") * (n - 1) / n
}
if (length(lv_int_names) == 0L) {
lav_msg_stop(gettext("lv.int.names is empty: no lv quadratic/interaction
terms are provided"))
}
if (is.null(lv_names)) {
lv_names <- paste("eta", seq_len(ncol(fs)), sep = "")
}
# missing data (missing = "ml"): the factor scores in fs were computed per
# missing pattern (see lav_sam_fs_missing()), and their measurement-error
# covariance B_p = M_p Theta_p M_p' varies over the patterns. pattern_list
# holds bp (the B_p, in the same metric as fs), wp (pattern proportions)
# and pat_of_case (the pattern of every row of fs). Under MCAR (pattern
# independent of eta), three things change relative to the complete-data
# case:
# (1) every var.error term LINEAR in the error covariance aggregates
# with B.bar = sum_p w_p B_p;
# (2) the quadratic term Var(e* %x% e*) = (I+K)(B %x% B) aggregates as
# (I+K) Q with Q = sum_p w_p (B_p %x% B_p) (NOT B.bar %x% B.bar);
# (3) E_p(f* %x% f*) shifts by vec(B_p) across the patterns, so the raw
# Var(fs2) picks up a between-pattern term
# D = sum_p w_p vec(B_p) vec(B_p)' - vec(B.bar) vec(B.bar)'
# that must be part of var.error as well.
# With a single pattern (B_p = B.bar, Q = B %x% B, D = 0) all of this
# reduces exactly to the complete-data expressions below.
mi_flag <- !is.null(pattern_list)
# MTM (B.bar over the patterns when mi_flag)
if (mi_flag) {
bp_aug <- lapply(pattern_list$bp, function(b) lav_mat_bdiag(0, b))
wp <- pattern_list$wp
pat_of_case <- pattern_list$pat_of_case
mtm <- Reduce(`+`, Map(function(b, w) w * b, bp_aug, wp))
} else {
mtm <- m %*% mm_theta %*% t(m)
}
# note: MTM no longer needs to be regularized (force-pd) here:
# lav_mat_sym_diff_smallest_root() handles a singular var.error matrix
# natively, by regressing out its null space (which includes the
# rows/columns of any dummy lvs); this also keeps MTM identical to the
# (raw) version used in lav_sam_gamma_add()
# augment to include intercept
fs <- cbind(1, fs)
n <- nrow(fs)
if (!mi_flag) {
mtm <- lav_mat_bdiag(0, mtm)
}
veta <- lav_mat_bdiag(0, veta)
eeta <- c(1, eeta)
lv_names <- c("..int..", lv_names)
nfac <- ncol(fs)
tmp <- lav_sam_lvnames2(lv_names)
idx1 <- tmp$idx1
idx2 <- tmp$idx2
names_1 <- tmp$names
fs2 <- fs[, idx1, drop = FALSE] * fs[, idx2, drop = FALSE]
k_nfac <- lav_mat_com(nfac, nfac)
ik <- diag(nfac * nfac) + k_nfac
eeta <- as.matrix(drop(eeta))
vetak_mtm <- veta %x% mtm
# normal version (for now):
if (mi_flag) {
# quadratic error moment: (I+K) Q with Q = sum_p w_p (B_p %x% B_p)
qfull <- Reduce(`+`, Map(function(b, w) w * (b %x% b), bp_aug, wp))
gamma_me22 <- ik %*% qfull
# between-pattern term D: vec(B_p)[a] in the (idx1, idx2) enumeration of
# (eta* %x% eta*) is B_p[idx1[a], idx2[a]]
vb <- vapply(bp_aug, function(b) b[cbind(idx1, idx2)],
numeric(length(idx1)))
vbbar <- mtm[cbind(idx1, idx2)]
dfull <- vb %*% (wp * t(vb)) - tcrossprod(vbbar)
} else {
gamma_me22 <- ik %*% (mtm %x% mtm)
dfull <- 0
}
# ingredients (normal ME case)
var_fs2 <- varn(fs2, n)
var_etak_me <- (tcrossprod(eeta) %x% mtm + vetak_mtm)
var_mek_eta <- lav_mat_com_pre_post(var_etak_me)
var_me2 <- gamma_me22
cov_etak_me_mek_eta <- lav_mat_com_post(var_etak_me)
cov_mek_eta_etak_me <- t(cov_etak_me_mek_eta)
var_error <- (var_etak_me + var_mek_eta + cov_etak_me_mek_eta
+ cov_mek_eta_etak_me + var_me2 + dfull)
# select only what we need
colnames(var_fs2) <- rownames(var_fs2) <- names_1
colnames(var_error) <- rownames(var_error) <- names_1
colnames(fs2) <- names_1
lv_keep <- c(lv_names[-1], lv_int_names)
var_fs2 <- var_fs2[lv_keep, lv_keep, drop = FALSE]
var_error <- var_error[lv_keep, lv_keep, drop = FALSE]
fs_mean <- colMeans(fs2, na.rm = TRUE)
names(fs_mean) <- names_1
fs2_mean <- fs_mean # all of them
fs_mean <- fs_mean[lv_keep]
# apply small sample correction (if requested)
alpha <- alpha_correction
alpha_n1 <- 0
if (alpha_correction > 0) {
alpha_n1 <- lav_sam_alpha_n1(alpha_correction, n)
var_error <- (1 - alpha_n1) * var_error
}
lambda <- +Inf
lambda_star <- 1
lambda_kappa <- 0
lambda_v <- NULL
lambda_floor_used <- as.numeric(NA)
if (lambda_correction && is.list(lambda_floor)) {
# internal (SE machinery only): frozen effective multiplier -- the
# numeric-jacobian channels re-run this function at perturbed inputs,
# and a kappa = 0 truncation (lambda-tracking debias margin, or a
# clamp at zero) must keep lambda_star fixed at its at-the-estimates
# value there
lambda_star <- lambda_floor[["star"]]
veta2 <- var_fs2 - lambda_star * var_error
} else if (lambda_correction) {
# use Fuller (1987) approach to ensure VETA2 is positive
lambda <- try(lav_mat_sym_diff_smallest_root(
var_fs2,
var_error
), silent = TRUE)
if (inherits(lambda, "try-error")) {
lav_msg_warn(gettext("failed to compute lambda"))
veta2 <- var_fs2 - var_error # and hope for the best
} else {
# determine the truncation floor:
# NULL -> historical rule (cutoff 1 + 2/n, floor p2/(n-1))
# a number -> coherent rule: truncate iff lambda < 1 + floor
floor_2 <- NULL
kappa_2 <- 1 # a fixed (numeric) floor: lambda_star tracks lambda
if (identical(lambda_floor, "debias")) {
# gated bootstrap-debiased floor; NULL when the sample margin is
# comfortably positive (then lambda >= 1 + 2/n by the gate
# construction, and the historical rule below applies the full
# correction: bit-for-bit the historical behavior) or when the
# setting is ineligible/degenerate (-> historical rule as well)
floor_2 <- lav_sam_veta2_floor_debias(
fs2 = fs2[, lv_keep, drop = FALSE],
var_error = var_error, lambda = lambda, n = n,
mi_flag = mi_flag
)
if (!is.null(floor_2) && !is.null(attr(floor_2, "kappa"))) {
kappa_2 <- attr(floor_2, "kappa")
}
} else if (is.numeric(lambda_floor)) {
floor_2 <- lambda_floor
}
if (!is.null(floor_2)) {
lambda_floor_used <- as.numeric(floor_2)
if (lambda < 1 + floor_2) {
lambda_star <- max(c(0, lambda - floor_2))
veta2 <- var_fs2 - lambda_star * var_error
# clamped at zero -> VETA2 = var_fs2 exactly, no lambda term
lambda_kappa <- if (lambda_star > 0) kappa_2 else 0
} else {
veta2 <- var_fs2 - var_error
}
} else {
# historical rule; note that the p2/(n-1) truncation margin
# vanishes with n, while the sample lambda of this
# high-dimensional fourth-order-moment pencil is severely
# downward-biased -- in weakly-identified conditions the
# truncation then binds over a long pre-asymptotic range of n
# with an ever smaller margin, and the accuracy of the
# structural estimates can *deteriorate* with increasing n;
# lambda2_floor = "debias" is a data-driven alternative
lambda_floor_used <- ncol(var_fs2) / (n - 1)
cutoff <- 1 + 2 / n # be more conservative for VETA2
if (lambda < cutoff) {
lambda_star <- max(c(0, lambda - ncol(var_fs2) / (n - 1)))
veta2 <- var_fs2 - lambda_star * var_error
lambda_kappa <- if (lambda_star > 0) 1 else 0
} else {
veta2 <- var_fs2 - var_error
}
}
# pencil eigenvector for the SE machinery: when the truncation
# binds with lambda_star tracking lambda (kappa = 1), the casewise
# contributions below gain the influence of lambda-hat itself
# (recompute of the same pencil: suppress its duplicate warnings)
if ((return_cov_iveta2 || extra) && lambda_kappa == 1) {
lam_v <- try(suppressWarnings(
lav_mat_sym_diff_smallest_root(var_fs2, var_error,
vector = TRUE)),
silent = TRUE)
vv <- attr(lam_v, "v")
if (!inherits(lam_v, "try-error") && !is.null(vv) &&
all(is.finite(vv))) {
lambda_v <- vv
}
# no usable direction -> the IF term is skipped (old behavior)
}
}
} else {
veta2 <- var_fs2 - var_error
}
# new in 0.6-20: to compute Gamma.eta
# the casewise contributions are defined so that their average over the
# observations equals *exactly* the augmented summary statistics
# (EETA2, VETA2) that are used in the second step; this requires the
# effective *first-order* lambda1 inside the decomposition (because
# EETA2 and var.error are built from the lambda1-corrected first-order
# VETA), while the *second-order* lambda_star multiplies the whole
# correction term:
# ieeta2_i = fs2[i, ] - lambda1 * vec(MTM)
# (the terms involving the mean factor scores cancel out exactly)
# and
# iveta2_i = tcrossprod(fs2[i, ] - fs2.mean) - lambda2.eff * tmp_i
# where lambda2.eff = lambda_star * (1 - alpha/(n-1)) is the effective
# multiplier of the *unscaled* var.error (the alpha correction rescales
# var.error before the lambda step), and
# tmp_i = ((F_i - lambda1 * MTM) %x% MTM) +
# (MTM %x% (F_i - lambda1 * MTM)) +
# lav_mat_com_post((F_i - lambda1 * MTM) %x% MTM) +
# lav_mat_com_pre((F_i - lambda1 * MTM) %x% MTM) +
# (I + K) %*% (MTM %x% MTM)
# with F_i = tcrossprod(f_i); both selected in the rows/columns we keep
# note: cov(iveta2_1) is invariant to the value of lambda1 (the lambda1
# terms are constant across the observations); but lambda1 does matter
# for the jacobian of the average casewise contribution computed in
# lav_sam_gamma_add()
# rewritten 12 June 2026 to avoid the loop over the observations: every
# kept element is a linear function of the second-order factor scores
# fs2[i, ]; elementwise, for kept elements a = (r1, s1) and b = (r2, s2)
# of (eta* %x% eta*):
# tmp_i[a, b] = F_i[r1, r2] MTM[s1, s2] + MTM[r1, r2] F_i[s1, s2] +
# F_i[r1, s2] MTM[s1, r2] + F_i[s1, r2] MTM[r1, s2] +
# (1 - 2 * lambda1) * (MTM[r1, r2] MTM[s1, s2] +
# MTM[s1, r2] MTM[r1, s2])
# where F_i[r, s] = fs2[i, (r - 1) * nfac + s]
if (return_cov_iveta2) {
lambda2_eff <- lambda_star * (1 - alpha_n1)
keep_idx <- match(lv_keep, names_1)
nkeep <- length(keep_idx)
# row/column pairs of the kept elements of (eta* %x% eta*)
i1k <- idx1[keep_idx]
i2k <- idx2[keep_idx]
# casewise contributions to E(eta* %x% eta*)[keep]
if (mi_flag) {
# pattern-specific error-mean shift: vec(B_p) replaces vec(B.bar), so
# that the between-pattern shift of E_p(f* %x% f*) does not inflate
# the covariance of the contributions
part1 <- fs2[, keep_idx, drop = FALSE] -
lambda1 * t(vb)[pat_of_case, keep_idx, drop = FALSE]
} else {
part1 <- t(t(fs2[, keep_idx, drop = FALSE]) -
lambda1 * lav_mat_vec(mtm)[keep_idx])
}
# casewise contributions to vech(Var(eta* %x% eta*)[keep, keep])
# vech row (pa) and column (pb) indices over the kept elements
tmp_mat <- matrix(seq_len(nkeep), nkeep, nkeep)
pa <- tmp_mat[lower.tri(tmp_mat, diag = TRUE)]
pb <- t(tmp_mat)[lower.tri(tmp_mat, diag = TRUE)]
r1 <- i1k[pa]; s1 <- i2k[pa]
r2 <- i1k[pb]; s2 <- i2k[pb]
# fs2 columns holding F_i[r1, r2], F_i[s1, s2], F_i[r1, s2], F_i[s1, r2]
col_11 <- (r1 - 1L) * nfac + r2
col_22 <- (s1 - 1L) * nfac + s2
col_12 <- (r1 - 1L) * nfac + s2
col_21 <- (s1 - 1L) * nfac + r2
b_11 <- mtm[cbind(r1, r2)]
b_22 <- mtm[cbind(s1, s2)]
b_12 <- mtm[cbind(r1, s2)]
b_21 <- mtm[cbind(s1, r2)]
fs2kc <- t(t(fs2[, keep_idx, drop = FALSE]) - fs2_mean[keep_idx])
if (mi_flag) {
# pattern-heterogeneous quadratic terms. The F-cross terms (with B.bar)
# average to the E-terms PLUS 2*lambda1*(bb1 + bb2); the -2*lambda1
# constant cancels this, and the per-pattern products average to the
# (I+K)Q quadratic part (q12) and the between-pattern D part (dvec),
# so the average over the observations equals the aggregated var.error
# above. With a single pattern this reduces exactly to the
# (1 - 2*lambda1) * (bb1 + bb2) constant of the complete-data branch.
q12_p <- vapply(bp_aug, function(b)
b[cbind(r1, r2)] * b[cbind(s1, s2)] +
b[cbind(s1, r2)] * b[cbind(r1, s2)], numeric(length(r1)))
dvec_p <- vapply(bp_aug, function(b)
b[cbind(r1, s1)] * b[cbind(r2, s2)], numeric(length(r1)))
const_bb <- -2 * lambda1 * (b_11 * b_22 + b_21 * b_12) -
mtm[cbind(r1, s1)] * mtm[cbind(r2, s2)]
quad_case <- t(q12_p + dvec_p)[pat_of_case, , drop = FALSE] +
rep(const_bb, each = n)
} else {
quad_case <- rep((1 - 2 * lambda1) * (b_11 * b_22 + b_21 * b_12),
each = n)
}
part2 <- (fs2kc[, pa, drop = FALSE] * fs2kc[, pb, drop = FALSE] -
lambda2_eff *
(t(t(fs2[, col_11, drop = FALSE]) * b_22) +
t(t(fs2[, col_22, drop = FALSE]) * b_11) +
t(t(fs2[, col_12, drop = FALSE]) * b_21) +
t(t(fs2[, col_21, drop = FALSE]) * b_12) +
quad_case))
# sampling variability of the second-order lambda-hat itself
# (kappa = 1 truncation): d lambda = t(v) (d var_fs2 -
# lambda * d var_error) v, with v the smallest-root pencil
# eigenvector (t(v) var_error v = 1). Casewise:
# IF_i = [ (z_i'v)^2 - mean ] -
# lambda * (1 - alpha_n1) * [ u_i - mean(u) ]
# where z_i are the centered kept second-order scores, and
# u_i = t(v) tmp_i v is the casewise var.error representation: with
# V the nfac x nfac embedding of v (V[i1k, i2k] = v) and
# Vs = V + t(V), the F-linear part collapses into the single form
# u_i = < F_i, Vs B Vs >
# (all four F x B index shapes of tmp_i; the case-constant parts
# cancel in the centering; under mi the per-pattern quadratic parts
# are added per case). The vech(VETA2) contributions then gain
# - var_error[pa, pb] * IF_i, which is exactly mean-zero: the
# statistics are unchanged, only their covariance changes.
if (lambda_kappa == 1 && !is.null(lambda_v)) {
zv <- drop(fs2kc %*% lambda_v)
v2mat <- matrix(0, nrow = nfac, ncol = nfac)
v2mat[cbind(i1k, i2k)] <- lambda_v
vs_m <- v2mat + t(v2mat)
w_m <- vs_m %*% mtm %*% vs_m
u_i <- drop(fs2 %*% as.vector(t(w_m)))
if (mi_flag) {
# per-pattern quadratic part (the (I+K)Q and D terms vary over
# the patterns; globally constant parts cancel in the centering)
u_quad <- vapply(bp_aug, function(b) {
q1 <- b[i1k, i1k, drop = FALSE] * b[i2k, i2k, drop = FALSE]
q2 <- b[i2k, i1k, drop = FALSE] * b[i1k, i2k, drop = FALSE]
drop(crossprod(lambda_v, (q1 + q2) %*% lambda_v)) +
sum(lambda_v * b[cbind(i1k, i2k)])^2
}, numeric(1L))
u_i <- u_i + u_quad[pat_of_case]
}
if_lambda <- (zv * zv - mean(zv * zv)) -
lambda * (1 - alpha_n1) * (u_i - mean(u_i))
part2 <- part2 - tcrossprod(if_lambda, var_error[cbind(pa, pb)])
}
iveta2_1 <- cbind(part1, part2)
cov_iveta2 <- cov(iveta2_1) * (n - 1) / n
}
# extra attributes?
if (extra) {
attr(veta2, "lambda") <- lambda
attr(veta2, "alpha") <- alpha
attr(veta2, "lambda_star") <- lambda_star
attr(veta2, "MSM") <- var_fs2
attr(veta2, "MTM") <- var_error
attr(veta2, "FS.mean") <- fs_mean
# SE machinery (see lav_sam_gamma_add() and the frozen-floor re-runs
# in the numeric jacobian channels)
attr(veta2, "lambda_kappa") <- lambda_kappa
attr(veta2, "lambda_floor_used") <- lambda_floor_used
}
if (return_fs) {
attr(veta2, "FS") <- fs2[, lv_keep, drop = FALSE]
}
if (return_cov_iveta2) {
attr(veta2, "cov.iveta2") <- cov_iveta2
# the casewise contributions themselves (needed for the casewise
# Gamma.eta, see lav_sam_gamma_add())
attr(veta2, "iveta2") <- iveta2_1
}
veta2
}
lav_sam_eeta2 <- function(eeta = NULL, veta = NULL, lv_names = NULL,
lv_int_names = NULL) {
if (length(lv_int_names) == 0L) {
lav_msg_stop(gettext("lv.int.names is empty: no lv quadratic/interaction
terms are provided"))
}
if (is.null(lv_names)) {
lv_names <- paste("eta", seq_len(ncol(veta)), sep = "")
}
nfac <- nrow(veta)
idx1 <- rep(seq_len(nfac), each = nfac)
idx2 <- rep(seq_len(nfac), times = nfac)
names_1 <- c(lv_names, paste(lv_names[idx1], lv_names[idx2], sep = ":"))
# E(\eta %x% \eta)
eeta2 <- lav_mat_vec(veta) + eeta %x% eeta
# add 1st order
eeta2_aug <- c(eeta, eeta2)
# select only what we need
names(eeta2_aug) <- names_1
lv_keep <- c(lv_names, lv_int_names)
eeta2_aug <- eeta2_aug[lv_keep]
eeta2_aug
}
# compute var(fs2) including quadratic/interaction terms
lav_sam_fs2 <- function(fs = NULL, lv_names = NULL, lv_int_names = NULL) {
varn <- function(x, n) {
var(x) * (n - 1) / n
}
if (length(lv_int_names) == 0L) {
lav_msg_stop(gettext("lv.int.names is empty: no lv quadratic/interaction
terms are provided"))
}
if (is.null(lv_names)) {
lv_names <- paste("eta", seq_len(ncol(fs)), sep = "")
}
# augment to include intercept
fs <- cbind(1, fs)
n <- nrow(fs)
lv_names <- c("..int..", lv_names)
nfac <- ncol(fs)
idx1 <- rep(seq_len(nfac), each = nfac)
idx2 <- rep(seq_len(nfac), times = nfac)
names_1 <- paste(lv_names[idx1], lv_names[idx2], sep = ":")
fs2 <- fs[, idx1, drop = FALSE] * fs[, idx2, drop = FALSE]
var_fs2 <- varn(fs2, n)
# select only what we need
colnames(var_fs2) <- rownames(var_fs2) <- names_1
lv_main <- paste(lv_names[-1], "..int..", sep = ":")
lv_keep <- c(lv_main, lv_int_names)
var_fs2 <- var_fs2[lv_keep, lv_keep, drop = FALSE]
var_fs2
}
# create consistent lavaan object, based on (filled in) PT
lav_sam_step3_joint <- function(fit = NULL, pt_1 = NULL, sam_method = "local") {
lavoptions <- fit@Options
lavoptions_joint <- lavoptions
lavoptions_joint$optim.method <- "none"
lavoptions_joint$optim.parscale <- "none"
lavoptions_joint$start <- "default"
lavoptions_joint$optim.force.converged <- TRUE
lavoptions_joint$check.gradient <- FALSE
lavoptions_joint$check.start <- FALSE
lavoptions_joint$check.post <- FALSE
lavoptions_joint$rotation <- "none"
lavoptions_joint$se <- "none"
lavoptions_joint$store.vcov <- FALSE # we do this manually
if (sam_method %in% c("local", "fsr", "cfsr")) {
lavoptions_joint$baseline <- FALSE
lavoptions_joint$fit.by.level <- FALSE
lavoptions_joint$sample.icov <- FALSE
#lavoptions.joint$h1 <- TRUE # we need this if we re-use the sam object
lavoptions_joint$test <- "none"
lavoptions_joint$estimator <- "none"
} else {
# global: the joint computes the core (base) test(s); "yuan.chan" is not a
# core lavaan test (it is added afterwards by lav_sam_global_test()), so we
# strip it here and make sure the unscaled "standard" base is present.
joint_test <- lavoptions$test[lavoptions$test != "yuan.chan"]
if (length(joint_test) == 0L || !"standard" %in% joint_test) {
joint_test <- unique(c("standard", joint_test))
}
lavoptions_joint$test <- joint_test
lavoptions_joint$estimator <- lavoptions$estimator
}
# set ustart values
pt_1$ustart <- pt_1$est # as this is used if optim.method == "none"
# the step-0 dummy fit of the local methods skips the inverse of S
# (sample.icov = FALSE: the inverse may not exist when N < nvar), but the
# joint-model computations of some estimators need S^-1 (eg the GLS
# information/objective). Backfill it here when it is missing; if S cannot
# be inverted it simply stays missing, exactly as before.
lavsamplestats_joint <- fit@SampleStats
if (fit@Data@nlevels == 1L && !fit@Model@categorical &&
!lavsamplestats_joint@missing.flag) { # core lavaan skips icov for
# missing = "ml" on purpose
for (g in seq_len(lavsamplestats_joint@ngroups)) {
if (is.null(lavsamplestats_joint@icov[[g]]) &&
!is.null(lavsamplestats_joint@cov[[g]])) {
# quietly: the inverse is optional here (only some estimators need
# it), so a non-PD S must not surface a warning at this point
current_warn <- lav_warn()
lav_warn(FALSE)
out <- try(lav_samp_icov(
cov_1 = lavsamplestats_joint@cov[[g]], ridge = 1e-05,
x_idx = lavsamplestats_joint@x.idx[[g]],
ngroups = fit@Data@ngroups, g = g
), silent = TRUE)
lav_warn(current_warn)
if (!inherits(out, "try-error")) {
lavsamplestats_joint@icov[[g]] <- out$icov
lavsamplestats_joint@cov.log.det[[g]] <- out$cov.log.det
}
}
}
}
joint <- lavaan::lavaan(pt_1,
slot_options = lavoptions_joint,
slot_sample_stats = lavsamplestats_joint,
slot_data = fit@Data,
verbose = FALSE
)
joint
}
# the (two-step robust) STRUCTURAL fit object of a local SAM model.
#
# The structural model FIT.PA is fitted with se = "robust.sem" using the
# step-1 induced NACOV (Gamma.eta) of the latent (co)variances/means. Its
# Satorra-Bentler scaled chi-square and the robust fit measures derived from it
# are the correct 'two-step robust' fit measures.
#
# For a single group, and for multiple groups WITHOUT across-group constraints
# in the measurement model (Case A), the per-group Gamma.eta list passed to
# FIT.PA already yields the correct robust test, and we simply reproduce it.
# For multiple groups WITH across-group measurement constraints (Case B), the
# step-1 variability couples the groups, so the relevant Gamma has nonzero
# cross-group blocks (step1$Gamma.eta.full) that a per-group NACOV list cannot
# represent. We then recompute the Satorra-Bentler scaled chi-square ourselves
# -- for the model AND the baseline -- using that full Gamma.
#
# The corrected FIT.PA is stored in the @internal slot of the final sam
# object (element sam.struc.fit.object), so that fitMeasures(), lavResiduals()
# and modindices() can delegate to it (issue #517).
lav_sam_struc_fit_object <- function(fit_pa = NULL, step1 = NULL) {
# do we have a (local) Gamma.eta and a non-saturated structural model?
have_gamma <- !is.null(step1$Gamma.eta) &&
length(step1$Gamma.eta) > 0L && !is.null(step1$Gamma.eta[[1]])
df_struc <- if (length(fit_pa@test) > 0L) fit_pa@test[[1]]$df else 0L
has_baseline <- !is.null(fit_pa@baseline) &&
!is.null(fit_pa@baseline$partable)
if (!have_gamma || df_struc < 1L || !has_baseline) {
return(fit_pa)
}
out <- try(
{
fp <- fit_pa
ng <- fp@Data@ngroups
ntot <- fp@SampleStats@ntotal
fg <- unlist(fp@SampleStats@nobs) / ntot
meanstr <- fp@Model@meanstructure
# full Gamma of the stacked structural statistics (cross-group blocks
# for Case B; block-diagonal otherwise)
if (isTRUE(step1$caseB)) {
gamma_all <- ntot * step1$Gamma.eta.full
} else {
gamma_all <- lav_mat_bdiag(lapply(
seq_len(ng), function(g) step1$Gamma.eta[[g]] / fg[g]
))
}
# model: Satorra-Bentler scaled test using the full Gamma
test_model <- lav_test_sb(
lavobject = fp, test = "satorra.bentler",
gamma_full = gamma_all
)
# baseline: refit the independence model on the latent (co)variances,
# then apply the same Satorra-Bentler correction with the full Gamma.
# (The corrected adjusted tests, if any, were added afterwards by
# lav_sam_step2_corrected_test() -- they need the casewise rows of
# Gamma.eta, which this moments-only refit does not have.)
opts_base <- fp@Options
opts_base$test <- setdiff(opts_base$test, lav_sam_corrected_tests)
if (length(opts_base$test) == 0L) {
opts_base$test <- "standard"
}
fit_base <- lavaan::lavaan(
model = fp@baseline$partable,
sample.cov = step1$VETA,
sample.mean = if (meanstr) step1$EETA else NULL,
sample.nobs = as.list(unlist(fp@SampleStats@nobs)),
nacov = step1$Gamma.eta,
slot_options = opts_base
)
test_base <- lav_test_sb(
lavobject = fit_base, test = "satorra.bentler",
gamma_full = gamma_all
)
# inject the corrected scaled tests
fp@test[["satorra.bentler"]] <- test_model[["satorra.bentler"]]
base_idx <- which(vapply(fp@baseline$test,
function(x) x$test, character(1L)) == "satorra.bentler")
if (length(base_idx) > 0L) {
fp@baseline$test[[base_idx[1]]] <- test_base[["satorra.bentler"]]
}
fp@external$sam.struc.corrected <- TRUE
fp
},
silent = TRUE
)
if (inherits(out, "try-error")) {
return(fit_pa) # the untouched original
}
out
}
# (two-step robust) fit measures for the structural part of a local SAM model,
# extracted from the (corrected) structural fit object
lav_sam_struc_fit <- function(fit_pa = NULL) {
measures <- c("chisq", "df", "pvalue", "cfi", "rmsea", "srmr")
fallback <- function() {
out <- try(fitMeasures(fit_pa, measures), silent = TRUE)
if (inherits(out, "try-error")) {
out <- "(unable to obtain fit measures)"
names(out) <- "warning"
}
out
}
if (!isTRUE(fit_pa@external$sam.struc.corrected)) {
return(fallback())
}
# use the fitMeasures() labels for the SCALED/ROBUST measures, so the
# printed header makes clear these are the two-step corrected values
# (chisq.scaled / pvalue.scaled / cfi.robust / rmsea.robust), exactly as
# fitMeasures() reports them. (The naive fallback() above keeps the plain
# chisq/pvalue/cfi/rmsea labels.)
scaled_labels <- c(
"chisq.scaled", "df.scaled", "pvalue.scaled",
"cfi.robust", "rmsea.robust", "srmr"
)
fm <- try(fitMeasures(fit_pa, scaled_labels), silent = TRUE)
if (inherits(fm, "try-error")) {
return(fallback())
}
# measures that fitMeasures() does not provide for this scaled test (eg
# the robust cfi/rmsea when the scaled test is a mean-var-adjusted one)
# are reported as NA
out <- rep(as.numeric(NA), length(scaled_labels))
names(out) <- scaled_labels
out[names(fm)] <- as.numeric(fm)
out
}
# Yuan & Chan (2002) rescaled GLOBAL test statistic for sam.method = "global".
#
# The 'joint' object already carries the base statistic T = n F[S, Sigma(theta,
# gamma.hat)] -- the discrepancy of the FULL (observed-variable) model evaluated
# at the plugged-in SAM estimates -- and the correct reference df = p* - q_total
# (all measurement AND structural parameters counted) in joint@test[[1]]. What is
# missing is the scaling correction tau.hat = tr(Q Gamma)/(p* - q) of eq. (22),
# which accounts for (a) the non-normality of the data (via the ADF Gamma) and
# (b) the fact that the measurement parameters gamma were estimated SEPARATELY in
# step 1 (segregated estimation), not jointly with theta.
#
# Key simplification that lets us reuse lav_test_sb() verbatim. With
# H.theta = sigma.dot.theta (sigma.dot.theta' W sigma.dot.theta)^-1
# sigma.dot.theta' W (the W-projector onto the
# structural columns of Delta),
# D = I - sigma.dot.gamma P (P = d gamma.hat / d stats, eq. 12a),
# Yuan & Chan's U (eq. 19) gives I - sigma.dot U = (I - H.theta) D, hence
# Q = (I - sigma.dot U)' W (I - sigma.dot U) = D' U.struc D ,
# with U.struc = W - W Delta.theta (Delta.theta' W Delta.theta)^-1 Delta.theta' W
# -- exactly the Satorra-Bentler residual-weight matrix, but built from the
# STRUCTURAL columns of Delta only. By the cyclic trace property
# tr(Q Gamma) = tr( U.struc . (D Gamma D') ).
# So we hand lav_test_sb():
# - delta = Delta.theta (structural columns of the full Delta)
# - e_inv = (sum_g f_g Delta.theta_g' W_g Delta.theta_g)^-1
# - wls_v = W (the estimator weight of the full model)
# - gamma_full = D Gamma D' (the 'measurement-residualized' Gamma)
# - test_unscaled = the standard (unscaled) full-model test (T and df)
# and read off the satorra.bentler-scaled statistic, which IS eq. (22).
#
# P (the step-1 measurement influence) is the SAME jacobian already used by the
# twostep.robust SEs (lav_sam_step1_local_jac(p_only = TRUE)). Multigroup support
# (Case A block-diagonal, Case B cross-group) follows from the stacked assembly:
# D and Gamma.tilde are built in the stacked statistics space, so the cross-group
# coupling that across-group measurement constraints induce in P is carried
# through automatically.
#
# Returns a list with two components, to be stored in the joint object:
# $test -> joint@test (standard + the YC-scaled test)
# $baseline.test -> joint@baseline$test (standard + a scaled baseline test)
# The baseline (independence) model has no measurement model, so the Yuan-Chan
# correction there reduces to the ordinary Satorra-Bentler scaling; we obtain it
# by refitting the baseline with test = "satorra.bentler". Both scaled tests are
# needed so that the robust/scaled fit measures (CFI, RMSEA) can be computed. On
# any failure the original (unscaled) test list is returned unchanged.
# 'test' is the internal SB-family base used to drive lav_test_sb (eq. 22 is the
# satorra.bentler-type rescaling); the resulting entry is renamed to "yuan.chan"
# (and given a Yuan-Chan label) in both the model and the baseline test lists.
lav_sam_global_test <- function(joint = NULL, step1 = NULL, step2 = NULL,
fit = NULL, test = "satorra.bentler") {
test_orig <- joint@test
# rename a satorra.bentler entry (and its $test field) to "yuan.chan"
rename_yc <- function(test_list) {
idx <- which(names(test_list) == "satorra.bentler")
if (length(idx) == 0L) {
return(test_list)
}
e <- test_list[[idx[1]]]
e$test <- "yuan.chan"
e$label <- "Yuan-Chan (2002) correction"
names(test_list)[idx[1]] <- "yuan.chan"
test_list[[idx[1]]] <- e
test_list
}
fallback <- list(test = test_orig, baseline.test = NULL)
# conditional.x: not supported (yet). The machinery below runs, but with
# NON-NORMAL (eg skewed or binary) covariates the scaled statistic is badly
# anti-conservative (empirical rejection rates up to ~0.3 at the 5% level
# in simulations, with either random or fixed covariate values; the
# simultaneous sem() satorra.bentler test on the same data is accurate, and
# so is the sam yuan.chan test with multivariate normal covariates -- the
# interplay between the unconditionally fitted step-1 blocks and the
# conditional full-model discrepancy is not yet correctly captured by the
# D Gamma D' residualization). Report the honest unscaled test instead.
# (The corrected STRUCTURAL test of sam.method = "local" does not have
# this problem: it is accurate in the same simulations.)
if (joint@Model@conditional.x) {
lav_msg_warn(gettext(
"the Yuan & Chan (2002) scaled test statistic (test = \"yuan.chan\") is
not available (yet) when conditional_x = TRUE; the standard (unscaled)
test is reported instead."))
return(fallback)
}
# locate the standard (unscaled) full-model test = base statistic T
std_idx <- which(vapply(test_orig, function(x) x$test, character(1L)) ==
"standard")
if (length(std_idx) == 0L) {
return(fallback) # no usable base statistic (eg test = "none")
}
std <- test_orig[[std_idx[1]]]
if (is.null(std$df) || anyNA(std$df) || std$df < 1L) {
return(fallback) # saturated full model: no global test
}
out <- try(
{
ngroups <- joint@Data@ngroups
ntot <- joint@SampleStats@ntotal
fg <- unlist(joint@SampleStats@nobs) / ntot
step1_free_idx <- step1$step1.free.idx
step2_free_idx <- step2$step2.free.idx
# ingredients of the full (joint) model:
# stored NACOV, else recompute with the fit-time settings. The
# (unbiased) ADF Gamma is unavailable for the fixed.x / conditional.x
# setting (lav_samp_gamma() with unbiased = TRUE stops), so fall back
# to the *biased* ADF Gamma there (gamma.unbiased = FALSE) -- exactly
# the Gamma that simultaneous sem(test = "satorra.bentler") uses for
# such models. Asymptotically equivalent; the simple (no-x) cases keep
# the unbiased Gamma so their values are unchanged.
gamma <- tryCatch(lav_gamma_used(joint), error = function(e) NULL)
if (is.null(gamma) || is.null(gamma[[1]])) {
opts <- joint@Options
opts$gamma.unbiased <- FALSE
gamma <- lav_gamma_used(
lavdata = joint@Data,
lavoptions = opts,
lavsamplestats = joint@SampleStats,
lavh1 = joint@h1,
lavimplied = joint@implied
)
}
delta <- lavTech(joint, "Delta")
wls_v <- lavTech(joint, "WLS.V")
# P = d(gamma.mm)/d(stats), reusing the twostep.robust jacobian; align its
# rows with step1.free.idx (P rows are in ascending free-parameter order)
p <- lav_sam_step1_local_jac(step1 = step1, fit = fit, p_only = TRUE)
p_row_idx <- match(step1_free_idx, attr(p, "free.idx"))
if (anyNA(p_row_idx)) {
lav_msg_stop(gettext(
"internal error: unable to align the step 1 jacobian (P) for the
global test statistic"))
}
p <- p[p_row_idx, , drop = FALSE]
# per group: structural and measurement columns of Delta, the weight W,
# and the SB-convention Gamma (Cov of sqrt(N) * per-group statistics)
delta_theta_list <- vector("list", ngroups)
delta_gamma_list <- vector("list", ngroups)
gamma_sb_list <- vector("list", ngroups)
a_theta <- matrix(0, length(step2_free_idx), length(step2_free_idx))
for (g in seq_len(ngroups)) {
delta_theta_list[[g]] <- delta[[g]][, step2_free_idx, drop = FALSE]
delta_gamma_list[[g]] <- delta[[g]][, step1_free_idx, drop = FALSE]
gamma_sb_list[[g]] <- gamma[[g]] / fg[g]
# structural information E.theta = sum_g f_g Delta.theta_g' W_g Delta.theta_g
wd <- lav_sam_wls_delta(wls_v[[g]], delta_theta_list[[g]])
a_theta <- a_theta + fg[g] * crossprod(delta_theta_list[[g]], wd)
}
e_theta_inv <- solve(a_theta)
# D (stacked) = I - Delta.gamma P ; Gamma.tilde = D Gamma D'
delta_gamma_all <- do.call(rbind, delta_gamma_list)
nstat_all <- nrow(delta_gamma_all)
d_mat <- diag(nstat_all) - delta_gamma_all %*% p
gamma_sb_all <- lav_mat_bdiag(gamma_sb_list)
gamma_tilde <- d_mat %*% gamma_sb_all %*% t(d_mat)
# reuse lav_test_sb: structural-only Delta/E.inv/W + residualized Gamma
# (gamma_full forces the assembled 'original' path)
lav_test_sb(
lavmodel = joint@Model,
lavsamplestats = joint@SampleStats,
lavoptions = joint@Options,
lavimplied = joint@implied,
lavdata = joint@Data,
test_unscaled = std,
e_inv = e_theta_inv,
delta = delta_theta_list,
wls_v = wls_v,
m_gamma = gamma, # ignored when gamma_full is set; avoids recompute
gamma_full = gamma_tilde,
test = test
)
},
silent = TRUE
)
if (inherits(out, "try-error")) {
# the Yuan-Chan correction reuses the same step-1 jacobian (P) as the
# robust SEs, so it shares their current limitations. Equality constraints
# (within-block and across-group) and exogenous covariates (fixed.x, and
# conditional.x with the biased ADF Gamma) ARE supported. Fall back to the
# unscaled test, with a note.
lav_msg_warn(gettext(
"the Yuan & Chan (2002) scaled test statistic (test = \"yuan.chan\") is
not available for this model; the standard (unscaled) test is reported
instead."))
return(fallback) # fall back to the unscaled standard test
}
# name the scaled entry "yuan.chan" (the SAM-global analogue of the
# simultaneous "satorra.bentler"); @test = list(standard, yuan.chan)
out <- rename_yc(out)
# baseline (independence) model: no measurement model -> the Yuan-Chan
# correction reduces to the ordinary Satorra-Bentler scaling. Refit the
# baseline with the matching scaled test so the robust CFI/RMSEA can be
# computed (and so fitMeasures()/summary() do not choke on a model that has a
# scaled test while its baseline does not).
baseline_test <- NULL
if (!is.null(joint@baseline$partable)) {
refit_baseline <- function(gamma_unbiased = NULL) {
opts <- fit@Options
opts$se <- "none"
opts$test <- test
opts$baseline <- FALSE
opts$fit.by.level <- FALSE
opts$estimator <- joint@Model@estimator
if (!is.null(gamma_unbiased)) {
opts$gamma.unbiased <- gamma_unbiased
}
fit_base <- lavaan::lavaan(
model = joint@baseline$partable,
slot_data = joint@Data,
slot_sample_stats = joint@SampleStats,
slot_options = opts,
verbose = FALSE
)
rename_yc(fit_base@test)
}
baseline_test <- try(refit_baseline(), silent = TRUE)
if (inherits(baseline_test, "try-error")) {
# the unbiased ADF Gamma (the sam default) is not available for eg
# fixed.x models: retry with the biased ADF Gamma, exactly as the
# model test above does (lav_gamma_used() fallback). Without a scaled
# baseline the scaled/robust incremental fit measures would be NA.
baseline_test <- try(refit_baseline(gamma_unbiased = FALSE),
silent = TRUE)
}
if (inherits(baseline_test, "try-error")) {
baseline_test <- NULL
}
}
list(test = out, baseline.test = baseline_test)
}
# Truncation-bias report: when the (first- or second-order) lambda
# truncation engaged (lambda_star < 1), the structural summary statistics
# retain (1 - lambda_star) * MTM (resp. var.error) of extra error
# (co)variance, and the step-2 estimates are shrunken toward zero. This
# helper collects the engaged state, plus a first-order (delta-method)
# approximation of the per-coefficient shrinkage bias
# bias = d theta / d stats . (1 - lambda_star) vech(MTM)
# i.e. theta at the truncated statistics minus theta at the
# (estimand-consistent) multiplier-1 statistics, linearized at the
# truncated fit. The result is stored in the @internal slot and reported
# as a note by summary(); the estimates themselves are NOT modified.
# The magnitude is only computed in the plain setting (single level, no
# conditional.x, no categorical, no across-group Case B); otherwise only
# the engaged state is returned. Returns NULL when nothing engaged.
lav_sam_trunc_bias <- function(step1 = NULL, fit_pa = NULL) {
lam <- unlist(step1$lambda)
if (is.null(lam) || !is.numeric(lam)) {
return(NULL)
}
engaged <- is.finite(lam) & lam < 1
if (!any(engaged)) {
return(NULL)
}
vn1 <- colnames(step1$VETA[[1]])
out <- list(
lambda_star = lam, engaged = engaged,
order = if (any(grepl(":", vn1, fixed = TRUE))) 2L else 1L,
bias = NULL
)
if (isTRUE(step1$caseB) || fit_pa@Model@conditional.x ||
fit_pa@Model@categorical || fit_pa@Data@nlevels > 1L) {
return(out)
}
out$bias <- tryCatch({
delta <- lavTech(fit_pa, "delta")
wls_v <- lavTech(fit_pa, "wls.v")
if (!is.list(delta)) delta <- list(delta)
if (!is.list(wls_v)) wls_v <- list(wls_v)
ng <- length(delta)
fg <- unlist(fit_pa@SampleStats@nobs) / fit_pa@SampleStats@ntotal
a_mat <- 0
rhs <- 0
for (g in seq_len(ng)) {
# the extra error (co)variance retained by the truncation, in the
# (possibly sub-block selected) order of the step-2 statistics; note
# that the variable order of the structural fit may differ from the
# order in VETA (the sample statistics are matched by name)
vn <- fit_pa@Data@ov.names[[g]]
if (!all(vn %in% colnames(step1$VETA[[g]]))) {
vn <- colnames(step1$VETA[[g]])
}
mtm <- step1$MTM[[g]]
if (!is.null(colnames(mtm))) {
idx <- match(vn, colnames(mtm))
} else {
idx <- match(vn, step1$LV.NAMES[[g]])
}
dstats <- numeric(nrow(delta[[g]]))
if (engaged[g] && !anyNA(idx)) {
vb <- (1 - lam[g]) *
lav_mat_vech(mtm[idx, idx, drop = FALSE])
# the vech(VETA) entries are the trailing block of the step-2
# statistics (any mean/intercept entries come first)
dstats[length(dstats) - length(vb) + seq_along(vb)] <- vb
}
vd <- wls_v[[g]] %*% delta[[g]]
a_mat <- a_mat + fg[g] * crossprod(delta[[g]], vd)
rhs <- rhs + fg[g] * crossprod(vd, dstats)
}
if (fit_pa@Model@eq.constraints) {
k_mat <- fit_pa@Model@eq.constraints.K
b_full <- drop(k_mat %*% solve(
crossprod(k_mat, a_mat %*% k_mat), crossprod(k_mat, rhs)))
} else {
b_full <- drop(solve(a_mat, rhs))
}
names(b_full) <- names(coef(fit_pa))
b_full
}, error = function(e) NULL)
# the exact counterpart: by how much do the structural estimates differ
# from the ones obtained under the historical truncation rule (floor
# 1/(n-1))? The linearized bias above is blind to the nonlinearity that
# characterizes the real danger zone (an ill-conditioned predictor block,
# where the historical estimates explode); there, this shift is large,
# while it is small if the truncation engaged for a harmless reason (eg
# an almost perfectly predicted latent outcome). First-order truncation
# only; one extra (point estimates only) fit of the structural part.
out$move <- tryCatch({
if (out$order != 1L) lav_msg_stop(gettext("First-order only."))
star_used <- unlist(step1$lambda1_star_used)
floor_used <- unlist(step1$lambda1_floor_used)
nobs <- unlist(fit_pa@SampleStats@nobs)
ng <- length(step1$VETA)
veta_h <- vector("list", ng)
changed <- FALSE
for (g in seq_len(ng)) {
veta_g <- step1$VETA[[g]]
veta_h[[g]] <- veta_g
if (!engaged[g]) next
if (!is.finite(star_used[g]) || star_used[g] <= 0 ||
!is.finite(floor_used[g])) {
lav_msg_stop(gettext("Raw lambda not available.")) # clamped at zero
}
lambda_raw <- star_used[g] + floor_used[g]
floor_h <- 1 / (nobs[g] - 1)
star_h <- if (lambda_raw < 1 + floor_h) lambda_raw - floor_h else 1
if (star_h - star_used[g] < 1e-12) next # same rule was active
vn <- colnames(veta_g)
mtm <- step1$MTM[[g]]
idx <- if (!is.null(colnames(mtm))) {
match(vn, colnames(mtm))
} else {
match(vn, step1$LV.NAMES[[g]])
}
if (anyNA(idx)) lav_msg_stop(gettext("Names do not match."))
d_std <- rep(1, length(idx))
if (!is.null(step1$D.STD) && !is.null(step1$D.STD[[g]])) {
d_std <- step1$D.STD[[g]][idx]
}
veta_raw <- t(veta_g * d_std) * d_std # undo the std.lv rescaling
veta_raw <- veta_raw -
(star_h - star_used[g]) * mtm[idx, idx, drop = FALSE]
if (any(d_std != 1)) {
d_new <- d_std
d_new[d_std != 1] <- sqrt(diag(veta_raw))[d_std != 1]
veta_raw <- t(veta_raw / d_new) / d_new
}
dimnames(veta_raw) <- dimnames(veta_g)
veta_h[[g]] <- veta_raw
changed <- TRUE
}
if (!changed) {
move <- 0 * coef(fit_pa)
} else {
pt_h <- as.list(fit_pa@ParTable)
pt_h$start <- pt_h$est
pt_h$est <- pt_h$se <- NULL
opt_h <- fit_pa@Options
opt_h$se <- "none"
opt_h$test <- "none"
opt_h$baseline <- FALSE
opt_h$h1 <- FALSE
opt_h$check.post <- FALSE
opt_h$check.gradient <- FALSE
opt_h$verbose <- FALSE
opt_h$warn <- FALSE
mean_h <- NULL
if (fit_pa@Model@meanstructure) {
mean_h <- fit_pa@SampleStats@mean
for (g in seq_along(mean_h)) {
names(mean_h[[g]]) <- fit_pa@Data@ov.names[[g]]
}
}
fit_h <- suppressWarnings(lavaan::lavaan(pt_h,
sample_cov = veta_h, sample_mean = mean_h,
sample_nobs = fit_pa@SampleStats@nobs,
slot_options = opt_h, verbose = FALSE
))
if (!lavInspect(fit_h, "converged")) {
lav_msg_stop(gettext("No convergence."))
}
move <- coef(fit_pa) - coef(fit_h)
}
move
}, error = function(e) NULL)
out
}
# Truncation report, part 2 (used by summary()): put the approximate
# shrinkage bias, and the exact shift away from the historical-rule
# estimates, on the scale of the standard errors (which are only available
# in the final object). Adds to the stored sam_trunc list:
# bias_max = max |bias| among the regression coefficients (or among all
# structural parameters if there are no regressions)
# bias_se_max = max |bias| / se among the same parameters
# move_max, move_se_max = idem, for the shift ('move')
# (NA if not available). The note printed by summary() is proportional: it
# is merely informative if BOTH bias_se_max and move_se_max are available
# and at most 0.5 (a bias of half a standard error lowers the coverage of a
# 95% confidence interval to about 92%); in all other cases -- including
# the settings where the shift cannot be computed -- the full warning is
# printed.
lav_sam_trunc_bias_se <- function(object = NULL) {
sam_trunc <- object@internal$sam_trunc
if (is.null(sam_trunc) || is.null(sam_trunc$bias)) {
return(sam_trunc)
}
pt <- object@ParTable
# same labels as coef(): user labels if any, else lhs-op-rhs (+ group
# suffix); equality-constrained duplicates share label, est and se
pt_labels <- tryCatch(lav_pt_labels(pt, type = "user"),
error = function(e) NULL)
if (length(pt_labels) != length(pt$lhs)) {
pt_labels <- paste0(pt$lhs, pt$op, pt$rhs)
}
# the regression coefficients (the variance parameters absorb most of the
# retained error variance by construction)
reg_names <- unique(pt_labels[pt$op == "~" & pt$free > 0L])
se <- NULL
if (!is.null(pt$se) && !identical(object@Options$se, "none")) {
se <- pt$se[pt$free > 0L]
names(se) <- pt_labels[pt$free > 0L]
}
max_abs <- function(v) {
out <- c(as.numeric(NA), as.numeric(NA))
if (is.null(v)) {
return(out)
}
if (any(names(v) %in% reg_names)) {
v <- v[names(v) %in% reg_names]
}
v <- v[is.finite(v)]
if (length(v) == 0L) {
return(out)
}
out[1] <- max(abs(v))
if (!is.null(se)) {
ratio <- abs(v) / se[match(names(v), names(se))]
ratio <- ratio[is.finite(ratio)]
if (length(ratio) > 0L) {
out[2] <- max(ratio)
}
}
out
}
tmp <- max_abs(sam_trunc$bias)
sam_trunc$bias_max <- tmp[1]
sam_trunc$bias_se_max <- tmp[2]
tmp <- max_abs(sam_trunc$move)
sam_trunc$move_max <- tmp[1]
sam_trunc$move_se_max <- tmp[2]
sam_trunc
}
lav_sam_table <- function(joint = NULL, step1 = NULL, fit_pa = NULL,
cmd = NULL, lavoptions = NULL,
mm_args = list(), struc_args = list(),
sam_method = "local",
local_options = list(), global_options = list()) {
mm_fit <- step1$MM.FIT
sam_mm_table <- data.frame(
Block = seq_along(step1$mm.list),
Latent = sapply(mm_fit, function(x) {
paste(unique(unlist(x@pta$vnames$lv)), collapse = ",")
}),
Nind = sapply(mm_fit, function(x) {
length(unique(unlist(x@pta$vnames$ov)))
}),
Chisq = sapply(mm_fit, function(x) {
x@test[[1]]$stat
}),
Df = sapply(mm_fit, function(x) {
x@test[[1]]$df
})
)
class(sam_mm_table) <- c("lavaan.data.frame", "data.frame")
# extra info for @internal slot
sam_trunc <- NULL
if (sam_method %in% c("local", "fsr", "cfsr")) {
# the structural fit object, with the (two-step robust) corrected test
# and baseline injected; for multigroup models with across-group
# measurement constraints (Case B) these are recomputed with the full
# cross-group Gamma. The object is stored in the @internal slot, so that
# fitMeasures(), lavResiduals() and modindices() can delegate to it
# (issue #517)
sam_struc_fit_object <- lav_sam_struc_fit_object(
fit_pa = fit_pa, step1 = step1
)
# (two-step robust) fit measures for the structural part
sam_struc_fit <- lav_sam_struc_fit(fit_pa = sam_struc_fit_object)
sam_mm_rel <- step1$REL
# engaged lambda truncation? collect the state + the approximate
# shrinkage bias of the structural estimates (summary() note)
sam_trunc <- lav_sam_trunc_bias(step1 = step1, fit_pa = fit_pa)
} else {
sam_struc_fit <- paste0("no local fit measures available for",
"structural part if sam.method is global")
names(sam_struc_fit) <- "warning"
sam_struc_fit_object <- NULL
sam_mm_rel <- numeric(0L)
}
# record the EFFECTIVE structural information.bread: lav_sam_step2() may
# have resolved "default" to "observed" (the default bread for the local
# standard errors since 0.7-2), and summary() reads it from
# sam.struc.args (see lavParameterEstimates)
if (is.null(struc_args[["information.bread"]]) &&
!is.null(fit_pa@Options$information.bread) &&
fit_pa@Options$information.bread != "default") {
struc_args$information.bread <- fit_pa@Options$information.bread
}
sam_1 <- list(
sam.cmd = cmd,
sam.method = sam_method,
sam.local.options = local_options,
sam.global.options = global_options,
sam.mm.list = step1$mm.list,
sam.mm.estimator = mm_fit[[1]]@Model@estimator,
sam.mm.args = mm_args,
sam.mm.ov.names = lapply(mm_fit, function(x) {
x@pta$vnames$ov
}),
sam.mm.table = sam_mm_table,
sam.mm.rel = sam_mm_rel,
sam_trunc = sam_trunc,
sam.struc.estimator = fit_pa@Model@estimator,
sam.struc.args = struc_args,
sam.struc.fit = sam_struc_fit,
sam.struc.fit.object = sam_struc_fit_object,
sam.lavoptions = lavoptions
)
sam_1
}
lav_sam_get_cov_ybar <- function(fit = NULL, local_options = list(
m_method = "ML",
lambda_correction = TRUE,
lambda1_floor = "debias",
lambda2_floor = "debias",
alpha_correction = 0L,
twolevel_method = "h1"
)) {
# local.twolevel.method
local_twolevel_method <- tolower(local_options[["twolevel_method"]])
if (!local_twolevel_method %in% c("h1", "anova", "mean")) {
lav_msg_stop(gettext(
"local option twolevel_method should be one of h1, anova or mean."))
}
lavpta <- fit@pta
nlevels <- lavpta$nlevels
nblocks <- lavpta$nblocks
# do we need H1?
if (nlevels > 1L && local_twolevel_method == "h1") {
h1 <- lav_h1_implied_logl(
lavdata = fit@Data,
lavsamplestats = fit@SampleStats,
lavoptions = fit@Options
)
h1implied <- h1$implied
} else {
h1implied <- fit@h1$implied
}
# containers
cov_list <- vector("list", nblocks)
ybar_list <- vector("list", nblocks)
slopes_list <- vector("list", nblocks)
# label
if (nblocks > 1L) {
names(cov_list) <- fit@Data@block.label
names(ybar_list) <- fit@Data@block.label
names(slopes_list) <- fit@Data@block.label
}
# collect COV/YBAR per block
for (b in seq_len(nblocks)) {
# get sample statistics for this block
if (nlevels > 1L) {
# block order is group-major: (g1l1, g1l2, g2l1, g2l2, ...)
this_level <- (b - 1L) %% nlevels + 1L
this_group <- floor((b - 1L) / nlevels) + 1L
if (this_level > 2L) {
lav_msg_stop(gettext("level 3 not supported (yet)."))
}
if (local_twolevel_method == "h1") {
# note: the h1 estimates are already reduced to the
# level-specific variables -- no need to subset
cov_1 <- h1implied$cov[[(this_group - 1L) * nlevels + this_level]]
ybar <- h1implied$mean[[(this_group - 1L) * nlevels + this_level]]
} else {
ylp <- fit@SampleStats@YLp[[this_group]][[2]]
if (this_level == 1L) {
# "anova" and "mean" use the same within statistics
cov_1 <- ylp$Sigma.W
ybar <- ylp$Mu.W
} else if (local_twolevel_method == "anova") {
cov_1 <- ylp$Sigma.B
ybar <- ylp$Mu.B
} else { # "mean"
s_pw <- ylp$Sigma.W
nj <- ylp$s
y2 <- ylp$Y2
# grand mean
mu_y <- ylp$Mu.W + ylp$Mu.B
y2c <- t(t(y2) - mu_y) # MUST be centered
yb <- crossprod(y2c) / nrow(y2c)
cov_1 <- yb - 1 / nj * s_pw
ybar <- ylp$Mu.B
}
# reduce to the level-specific variables
ov_idx <- fit@Data@Lp[[this_group]]$ov.idx[[this_level]]
cov_1 <- cov_1[ov_idx, ov_idx, drop = FALSE]
ybar <- ybar[ov_idx]
}
# single level
} else {
this_group <- b
if (fit@Model@conditional.x) {
ybar <- h1implied$res.int[[b]]
cov_1 <- h1implied$res.cov[[b]]
slopes_list[[b]] <- h1implied$res.slopes[[b]]
} else {
ybar <- h1implied$mean[[b]] # EM version if missing="ml"
cov_1 <- h1implied$cov[[b]]
}
} # single level
cov_list[[b]] <- cov_1
ybar_list[[b]] <- ybar
}
list(COV = cov_list, YBAR = ybar_list, SLOPES = slopes_list)
}
# automatically generate mm.list; we create measurement blocks so that:
# - unlinked factors become a singleton
# - linked factor are joined into a single measurement block
# based: cross-loadings, or correlated residuals between item errors
lav_sam_get_mmlist <- function(lavobject) {
lavmodel <- lavobject@Model
lavpta <- lavobject@pta
nblocks <- lavpta$nblocks
lambda_idx <- which(names(lavmodel@GLIST) == "lambda")
mm_list <- vector("list", length = nblocks)
glist <- lavTech(lavobject, "partable")
for (b in seq_len(nblocks)) {
mm_in_block <- (seq_len(lavmodel@nmat[b]) +
cumsum(c(0, lavmodel@nmat))[b])
mlist <- glist[mm_in_block]
cc <- t(mlist$lambda) %*% mlist$theta %*% mlist$lambda
# note: CC contains dummy lv's, higher-order, etc... and they
# must be removed
# get ALL lv names (including dummy ov.x/ov.y)
ov_names <- lavmodel@dimNames[[lambda_idx[b]]][[1L]]
lv_names <- lavmodel@dimNames[[lambda_idx[b]]][[2L]]
names_1 <- lv_names
# needs to removed:
# - all lv.ind variables
# - all higher-order latent variables
# - all dummy lv's
rm_idx <- c(match(lavpta$vnames$lv.ind[[b]], lv_names),
match(lavpta$vnames$lv.interaction[[b]], lv_names),
which(lv_names %in% ov_names))
if (length(rm_idx) > 0L) {
cc <- cc[-rm_idx, -rm_idx, drop = FALSE]
names_1 <- lv_names[-rm_idx]
}
# cluster membership
membership <- lav_graph_get_connected_nodes(cc)
out <- split(names_1, names_1[membership])
names(out) <- paste("block", seq_along(out), sep = "")
mm_list[[b]] <- out
}
mm_list
}
# labels for the elements of the WLS.obs statistics vector in the
# categorical case: 1) thresholds and (negative) means, interleaved per
# variable, 2) [conditional.x only] the slopes of the y-variables on the
# exogenous covariates, column-major (lav_mat_vec(SLOPES): all y for x1,
# then x2, ...), 3) variances (numeric variables only), 4) correlations
# (vech, no diagonal); used to map the statistics of a measurement block
# into the statistics vector of the joint model
lav_sam_wls_obs_labels <- function(ov_names, th_idx, ov_names_x = NULL) {
nvar <- length(ov_names)
# variable of each th element: ordered variables have th_idx > 0; the
# th_idx == 0 elements are the (negative) means of the numeric
# variables, appearing in variable order
v_of <- th_idx
num_var_idx <- which(!seq_len(nvar) %in% th_idx)
v_of[th_idx == 0L] <- num_var_idx
th_labels <- character(length(th_idx))
count <- integer(nvar)
for (k in seq_along(th_idx)) {
v <- v_of[k]
count[v] <- count[v] + 1L
th_labels[k] <- paste0(ov_names[v], "|t", count[v])
}
# slopes (conditional.x only), column-major over the nvar x nexo matrix
slope_labels <- character(0L)
if (length(ov_names_x) > 0L) {
slope_labels <- paste0(rep(ov_names, times = length(ov_names_x)), "~",
rep(ov_names_x, each = nvar))
}
# variances (numeric variables only)
var_labels <- character(0L)
if (length(num_var_idx) > 0L) {
var_labels <- paste0(ov_names[num_var_idx], "|var")
}
# correlations: lav_mat_vech(cov, diagonal = FALSE), column-major;
# use a canonical (sorted) pair label, as the relative order of two
# variables may differ between the joint model and a measurement block
cor_labels <- character(0L)
if (nvar > 1L) {
row_idx <- matrix(seq_len(nvar), nvar, nvar)[lower.tri(diag(nvar))]
col_idx <- matrix(seq_len(nvar), nvar, nvar,
byrow = TRUE)[lower.tri(diag(nvar))]
cor_labels <- vapply(seq_along(row_idx), function(k) {
paste(sort(c(ov_names[row_idx[k]], ov_names[col_idx[k]])),
collapse = "~~")
}, character(1L))
}
c(th_labels, slope_labels, var_labels, cor_labels)
}
# conditional.x: permutation of the conditional structural-moment vector
# c( vecr(cbind(EETA, RS_eta)), vech(VETA) ) -- with the variables in
# VETA (measurement) order and the covariates in data (ov.names.x) order --
# to the structural model's internal order (target_eta = its ov.nox,
# target_x = its ov.x). This is the layout of FIT.PA's WLS.obs: the NACOV
# passed to lavaan() is consumed as-is (no name matching), so Gamma.eta must
# be permuted with this index before the step-2 handoff.
lav_sam_condx_perm_idx <- function(eta_names, x_names,
target_eta, target_x) {
ord <- match(target_eta, eta_names)
ord_x <- match(target_x, x_names)
if (length(ord) != length(eta_names) || anyNA(ord) ||
length(ord_x) != length(x_names) || anyNA(ord_x)) {
lav_msg_stop(gettext(
"internal error: unable to align the conditional structural moments
with the structural model (Gamma.eta permutation)"))
}
p <- length(eta_names)
q <- length(x_names)
mean_block <- unlist(lapply(ord, function(j) {
(j - 1L) * (q + 1L) + c(1L, 1L + ord_x)
}))
pos <- matrix(0L, p, p)
pos[lower.tri(pos, diag = TRUE)] <- seq_len(p * (p + 1L) / 2L)
pos[upper.tri(pos)] <- t(pos)[upper.tri(pos)]
sub <- pos[ord, ord, drop = FALSE]
vech_perm <- sub[lower.tri(sub, diag = TRUE)]
c(mean_block, p * (q + 1L) + vech_perm)
}
# merge the per-block mm clusters (as returned by lav_sam_get_mmlist())
# into a single mm.list, transitively joining clusters that share a latent
# variable or an observed indicator; needed for multilevel models, where
# the within and between factors of the same indicators must end up in the
# same (two-level) measurement block
lav_sam_mmlist_merge_blocks <- function(lavobject, mm_list_per_block) {
pt_1 <- lavobject@ParTable
nblocks <- lavobject@pta$nblocks
block_values <- lav_pt_block_values(pt_1)
# collect all clusters: a set of lv names + their observed indicators
clusters <- list()
for (b in seq_len(nblocks)) {
for (cl in mm_list_per_block[[b]]) {
lv <- unlist(cl)
ind_idx <- which(pt_1$op == "=~" & pt_1$block == block_values[b] &
pt_1$lhs %in% lv)
clusters <- c(clusters,
list(list(lv = lv, ov = unique(pt_1$rhs[ind_idx]))))
}
}
# merge until no two clusters share an lv or an ov
changed <- TRUE
while (changed) {
changed <- FALSE
i <- 1L
while (i < length(clusters)) {
j <- i + 1L
while (j <= length(clusters)) {
if (any(clusters[[i]]$lv %in% clusters[[j]]$lv) ||
any(clusters[[i]]$ov %in% clusters[[j]]$ov)) {
clusters[[i]]$lv <- unique(c(clusters[[i]]$lv, clusters[[j]]$lv))
clusters[[i]]$ov <- unique(c(clusters[[i]]$ov, clusters[[j]]$ov))
clusters[[j]] <- NULL
changed <- TRUE
} else {
j <- j + 1L
}
}
i <- i + 1L
}
}
lapply(clusters, "[[", "lv")
}
lav_sam_veta_pt <- function(lavobject, block = 1L) {
lavmodel <- lavobject@Model
lambda_idx <- which(names(lavmodel@GLIST) == "lambda")
glist <- lavTech(lavobject, "partable")
b <- block
mm_in_block <- (seq_len(lavmodel@nmat[b]) + cumsum(c(0, lavmodel@nmat))[b])
mlist <- glist[mm_in_block]
mm_psi <- (mlist$psi != 0) + 0L
nr <- nrow(mlist$psi)
if (!is.null(mlist$beta)) {
mm_beta <- (mlist$beta != 0) + 0L
tmp <- -mm_beta
tmp[lav_mat_diag_idx(nr)] <- 1
ib_inv <- solve(tmp)
} else {
ib_inv <- diag(nr)
}
veta <- ib_inv %*% mm_psi %*% t(ib_inv)
# get ALL lv names (including dummy ov.x/ov.y)
lv_names <- lavmodel@dimNames[[lambda_idx[b]]][[2L]]
colnames(veta) <- rownames(veta) <- lv_names
veta
}
lav_sam_veta_con <- function(s = NULL, mm_lambda = NULL, mm_theta = NULL,
l_veta = NULL, local_m_method = "ML",
tol = 1e-07, max_iter = 100L) {
# first do GLS/ULS
m <- ncol(mm_lambda)
if (local_m_method == "ULS") {
m_w <- diag(m)
} else {
m_w <- solve(s)
}
wl <- m_w %*% mm_lambda
t_lw <- t(wl)
t_lwl <- t_lw %*% mm_lambda
tmp1 <- t(l_veta) %*% (t_lwl %x% t_lwl) %*% l_veta
tmp2 <- t(l_veta) %*% lav_mat_vec(t_lw %*% (s - mm_theta) %*% wl)
out <- solve(tmp1, tmp2)
veta_init <- matrix(l_veta %*% out, m, m)
if (local_m_method != "ML") {
# we are done!
return(veta_init)
}
veta_new <- veta_init
for (i in seq_len(max_iter)) {
sigma_ml <- mm_lambda %*% veta_new %*% t(mm_lambda) + mm_theta
m_w <- solve(sigma_ml)
wl <- m_w %*% mm_lambda
t_lw <- t(wl)
t_lwl <- t_lw %*% mm_lambda
tmp1 <- t(l_veta) %*% (t_lwl %x% t_lwl) %*% l_veta
tmp2 <- t(l_veta) %*% lav_mat_vec(t_lw %*% (s - mm_theta) %*% wl)
out <- solve(tmp1, tmp2)
veta_ml <- matrix(l_veta %*% out, m, m)
rmsea <- sqrt(sum((veta_new - veta_ml)^2))
#cat("i = ", i, " rmsea = ", rmsea, "\n")
veta_new <- veta_ml
if (rmsea < tol) {
break
}
}
veta_new
}
Try the lavaan package in your browser
Any scripts or data that you put into this service are public.
lavaan documentation built on Oct. 8, 2026, 5:06 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions lav_sam_veta_con lav_sam_veta_pt lav_sam_mmlist_merge_blocks lav_sam_condx_perm_idx lav_sam_wls_obs_labels lav_sam_get_mmlist lav_sam_get_cov_ybar lav_sam_table lav_sam_trunc_bias_se lav_sam_trunc_bias lav_sam_global_test lav_sam_struc_fit lav_sam_struc_fit_object lav_sam_step3_joint lav_sam_fs2 lav_sam_eeta2 lav_sam_veta2 lav_sam_veta2_floor_debias lav_sam_fs_missing lav_sam_eeta lav_sam_veta lav_sam_veta1_floor_debias lav_sam_lambda_floor_debias lav_sam_tmat lav_sam_mapping_mat_tmat lav_sam_mapping_mat lav_sam_struc_object lav_sam_local_flag lav_sam_lvnames2 lav_sam_alpha_n1
# utility functions for the sam() function
# YR 4 April 2023
# clamped small-sample correction factor alpha/(n - 1), forced into [0, 1]
lav_sam_alpha_n1 <- function(alpha_correction = 0L, n = 20L) {
alpha_n1 <- alpha_correction / (n - 1)
if (!is.finite(alpha_n1)) {
alpha_n1 <- 0
}
min(max(alpha_n1, 0), 1)
}
# row/column index pairs and labels for the elements of (eta* %x% eta*),
# where lv_names already includes the leading "..int.." intercept entry;
# the first nfac labels are the plain (first-order) lv names
lav_sam_lvnames2 <- function(lv_names) {
nfac <- length(lv_names)
idx1 <- rep(seq_len(nfac), each = nfac)
idx2 <- rep(seq_len(nfac), times = nfac)
names2 <- paste(lv_names[idx1], lv_names[idx2], sep = ":")
names2[seq_len(nfac)] <- lv_names
list(idx1 = idx1, idx2 = idx2, names = names2)
}
# is this a sam object created with sam.method = "local", "fsr" or "cfsr"?
lav_sam_local_flag <- function(object) {
sam_method <- object@internal$sam.method
!is.null(sam_method) && sam_method %in% c("local", "fsr", "cfsr")
}
# the stored step-2 structural fit of a local SAM object (with the two-step
# corrected test/baseline injected, if available); NULL for global sam
# objects, and for sam objects created by older versions of lavaan
lav_sam_struc_object <- function(object) {
object@internal$sam.struc.fit.object
}
# construct 'mapping matrix' M using either "ML", "GLS" or "ULS" method
#
# by construction, M %*% LAMBDA = I (the identity matrix)
lav_sam_mapping_mat <- function(mm_lambda = NULL, mm_theta = NULL,
s = NULL, s_inv = NULL,
method = "ML") {
method <- toupper(method)
# catch empty columns in LAMBDA (eg higher order!)
lambda_orig <- mm_lambda
empty_idx <- which(apply(mm_lambda, 2L, function(x) all(x == 0)))
if (length(empty_idx) > 0L) {
mm_lambda <- lambda_orig[, -empty_idx, drop = FALSE]
}
# ULS
# M == solve( t(LAMBDA) %*% LAMBDA ) %*% t(LAMBDA)
# == MASS:::ginv(LAMBDA)
if (method == "ULS") {
m <- try(tcrossprod(solve(crossprod(mm_lambda)), mm_lambda),
silent = TRUE
)
if (inherits(m, "try-error")) {
lav_msg_warn(gettext(
"cannot invert crossprod(LAMBDA); using generalized inverse"))
m <- MASS::ginv(mm_lambda)
}
# GLS
# M == solve( t(LAMBDA) %*% S.inv %*% LAMBDA ) %*% t(LAMBDA) %*% S.inv
} else if (method == "GLS") {
if (is.null(s_inv)) {
s_inv <- try(solve(s), silent = TRUE)
}
if (inherits(s_inv, "try-error")) {
lav_msg_warn(gettext("S is not invertible; switching to ULS method"))
m <- lav_sam_mapping_mat(mm_lambda = mm_lambda, method = "ULS")
} else {
t_lsinv <- t(mm_lambda) %*% s_inv
t_lsinv_l <- t_lsinv %*% mm_lambda
m <- try(solve(t_lsinv_l, t_lsinv), silent = TRUE)
if (inherits(m, "try-error")) {
lav_msg_warn(gettext("problem constructing mapping matrix;
switching to generalized inverse"))
m <- MASS::ginv(t_lsinv_l) %*% t_lsinv
}
}
# ML
# M == solve(t(LAMBDA) %*% THETA.inv %*% LAMBDA) %*% t(LAMBDA) %*% THETA.inv
} else if (method == "ML") {
# Problem: if THETA has zero elements on the diagonal, we cannot invert
# As we do not have access to Sigma(.inv), we cannot
# use the trick as in lavPredict() (where we replace THETA.inv by
# Sigma.inv)
# old method (<0.6-16): remove rows/cols with zero values
# on the diagonal of THETA, invert the submatrix and
# set the '0' diagonal elements to one. This resulted in (somewhat)
# distorted results.
# new in 0.6-16: we use the Wall & Amemiya (2000) method using
# the so-called 'T' transformation
# new in 0.6-18: if we cannot use the marker method (eg growth models),
# use the 'old' method anyway: remove zero rows/cols and invert submatrix
zero_theta_idx <- which(abs(diag(mm_theta)) < 1e-4) # be conservative
if (length(zero_theta_idx) == 0L) {
# ok, no zero diagonal elements: try to invert THETA
if (lav_mat_is_diagonal(mm_theta)) {
theta_inv <- diag(1 / diag(mm_theta), nrow = nrow(mm_theta))
} else {
theta_inv <- try(solve(mm_theta), silent = TRUE)
if (inherits(theta_inv, "try-error")) {
theta_inv <- NULL
}
}
} else {
# see if we can use marker method
marker_idx <- lav_utils_get_marker(mm_lambda = mm_lambda, std_lv = TRUE)
if (any(is.na(marker_idx))) {
theta_inv <- try(lav_mat_sym_inverse(mm_theta), silent = TRUE)
if (inherits(theta_inv, "try-error")) {
theta_inv <- NULL
} else {
diag(theta_inv)[zero_theta_idx] <- 1
}
} else {
# try tmat method later on
theta_inv <- NULL
}
}
# could we invert THETA?
if (!is.null(theta_inv)) {
# ha, all is good; compute M the usual way
t_lti <- t(mm_lambda) %*% theta_inv
t_lti_l <- t_lti %*% mm_lambda
m <- try(solve(t_lti_l, t_lti), silent = TRUE)
if (inherits(m, "try-error")) {
lav_msg_warn(gettext(
"problem constructing ML mapping matrix; switching to ULS"))
m <- lav_sam_mapping_mat(mm_lambda = mm_lambda, method = "ULS")
}
} else {
# use W&A2000's method using the 'T' transformation
m <- try(lav_sam_mapping_mat_tmat(
mm_lambda = mm_lambda,
mm_theta = mm_theta
), silent = TRUE)
if (inherits(m, "try-error")) {
lav_msg_warn(gettext(
"problem constructing ML mapping matrix; switching to ULS"))
m <- lav_sam_mapping_mat(mm_lambda = mm_lambda, method = "ULS")
}
}
} # ML
# empty.idx?
if (length(empty_idx) > 0L) {
m_full <- m
m <- matrix(0, nrow = ncol(lambda_orig), ncol = ncol(m_full))
m[-empty_idx, ] <- m_full
}
m
}
# use 'T' transformation to create the 'Bartlett/ML' mapping matrix
# see Wall & Amemiya (2000) eq (7)
# see also Fuller 1987 page 357 (where T is called H), and page 364
# although Fuller/W&A always assumed that THETA is diagonal,
# their method seems to work equally well for non-diagonal THETA
#
# in our implementation:
# - we do NOT reorder the rows of LAMBDA
# - if std.lv = TRUE, we first rescale to 'create' marker indicators
# and then rescale back at the end
#
lav_sam_mapping_mat_tmat <- function(mm_lambda = NULL,
mm_theta = NULL,
marker_idx = NULL,
std_lv = NULL) {
mm_lambda <- as.matrix.default(mm_lambda)
# catch empty columns in LAMBDA (eg higher order!)
lambda_orig <- mm_lambda
empty_idx <- which(apply(mm_lambda, 2L, function(x) all(x == 0)))
if (length(empty_idx) > 0L) {
mm_lambda <- lambda_orig[, -empty_idx, drop = FALSE]
}
# nvar <- nrow(mm_lambda)
# nfac <- ncol(mm_lambda)
# do we have marker.idx?
if (is.null(marker_idx)) {
# 'marker' indicator has a single non-zero element in a row
marker_idx <- lav_utils_get_marker(mm_lambda = mm_lambda, std_lv = TRUE)
if (any(is.na(marker_idx))) {
lav_msg_stop(gettext("no clear markers in LAMBDA matrix"))
}
}
# std.lv TRUE or FALSE?
if (is.null(std_lv)) {
std_lv <- FALSE
if (any(diag(mm_lambda[marker_idx, , drop = FALSE]) != 1)) {
std_lv <- TRUE
}
}
# if std.lv = TRUE, rescale
if (std_lv) {
marker <- mm_lambda[marker_idx, , drop = FALSE]
marker_inv <- 1 / diag(marker)
mm_lambda <- t(t(mm_lambda) * marker_inv)
}
# compute 'T' matrix
tmat <- lav_sam_tmat(
mm_lambda = mm_lambda, mm_theta = mm_theta,
marker_idx = marker_idx
)
# ML mapping matrix
m <- tmat[marker_idx, , drop = FALSE]
if (std_lv) {
m <- m * marker_inv
}
# empty.idx?
if (length(empty_idx) > 0L) {
m_full <- m
m <- matrix(0, nrow = ncol(lambda_orig), ncol = ncol(m_full))
m[-empty_idx, ] <- m_full
}
m
}
# create 'T' matrix (tmat) for T-transformation
#
# Notes: - here we assume that LAMBDA has unity markers (no std.lv = TRUE)
# - TMAT is NOT symmetric!
# - Yc %*% t(TMAT) transforms the data in such a way that we get:
# 1) Bartlett factor scores in the marker columns
# 2) 'V' values in the non-marker columns, where:
# V = Yc - Yc[,marker.idx] %*% t(LAMBDA)
#
lav_sam_tmat <- function(mm_lambda = NULL,
mm_theta = NULL,
marker_idx = NULL) {
mm_lambda <- as.matrix.default(mm_lambda)
nvar <- nrow(mm_lambda)
nfac <- ncol(mm_lambda)
# do we have marker.idx?
if (is.null(marker_idx)) {
# 'marker' indicator has a single 1 element in a row
marker_idx <- lav_utils_get_marker(mm_lambda = mm_lambda, std_lv = FALSE)
if (any(is.na(marker_idx))) {
lav_msg_stop(gettext("no clear markers in LAMBDA matrix"))
}
}
# construct 'm_c' matrix
c2 <- diag(nvar)
c2[, marker_idx] <- -1 * mm_lambda
m_c <- c2[-marker_idx, , drop = FALSE]
# compute Sigma.ve and Sigma.vv
sigma_ve <- m_c %*% mm_theta
# Sigma.vv <- m_c %*% THETA %*% t(m_c)
sigma_vv <- sigma_ve %*% t(m_c)
# construct 'Gamma' (=gamma_1) (and gamma2) matrix
# Gamma <- (t(Sigma.ve) %*% solve(Sigma.vv))[marker.idx,, drop = FALSE]
gamma_1 <- try(t(solve(sigma_vv, sigma_ve)[, marker_idx, drop = FALSE]),
silent = TRUE
)
if (inherits(gamma_1, "try-error")) {
tmp <- t(sigma_ve) %*% MASS::ginv(sigma_vv)
gamma_1 <- tmp[marker_idx, , drop = FALSE]
}
gamma2 <- matrix(0, nfac, nvar)
gamma2[, -marker_idx] <- gamma_1
gamma2[, marker_idx] <- diag(nfac)
# transformation matrix 'T' (we call it here 'Tmat')
tmat <- matrix(0, nvar, nvar)
tmat[-marker_idx, ] <- m_c
tmat[marker_idx, ] <- -gamma2 %*% c2
tmat
}
# Gated, bootstrap-debiased truncation floor for the lambda corrections --
# the shared core behind lav_sam_veta1_floor_debias() (first order, see
# lav_sam_veta()) and lav_sam_veta2_floor_debias() (second order, see
# lav_sam_veta2()).
#
# 'score' holds the casewise contributions whose (n-divisor) covariance is
# the left-hand side of the pencil: the Bartlett factor scores
# f_i = M (y_i - ybar) for the first order (Var(f) = MSM), and the
# second-order products fs2_i (kept columns) for the second order
# (Var(fs2) = var_fs2); 'err' is the right-hand side of the pencil (MTM
# resp. var.error), and lambda its smallest generalized root.
#
# The sample lambda is downward-biased -- mildly in well-separated
# settings, severely when the bottom of the spectrum is clustered (highly
# correlated factors, weak reliability, many product terms). The
# historical rules truncate with floors (1/(n-1) resp. p2/(n-1)) that
# vanish with n, so in difficult settings the accuracy of the structural
# estimates can deteriorate with increasing n. This rule is designed to
# be a strict no-op whenever the sample is comfortably away from the
# boundary, and to repair the difficult cases:
#
# 1. GATE: six deterministic interleaved half-samples give a cheap noise
# estimate se(lambda) = sd(lambda_half)/sqrt(2). If
# lambda - 1 >= max(2 * se(lambda), floor_default)
# the sample margin is comfortably positive: return NULL, and the
# caller applies the historical rule -- which then performs the full
# correction, bit-for-bit identical to the historical behavior.
# 2. Otherwise (the danger zone), estimate the bias and the sampling
# noise of lambda by a casewise bootstrap of the score rows (B = 200
# resamples; a fixed internal seed makes the resampling plan
# deterministic, and the user's .Random.seed is saved and restored):
# bias = max(0, mean(lambda_boot) - lambda)
# se_b = sd(lambda_boot)
# and return the floor
# max(min(lambda - 1 + bias, lambda - 1 + 3 * se_b),
# 2 * se_b, floor_default)
# i.e. the debiased margin, capped at three standard errors above the
# sample margin, and never below (twice) the noise level of lambda nor
# below the historical floor. The caller truncates iff
# lambda < 1 + floor, with lambda_star = max(0, lambda - floor).
#
# Return value: a numeric floor (danger zone only), or NULL. NULL covers
# both the open gate and any ineligible/degenerate situation (no usable
# scores, n < 4 * ncol(score), failed sub-sample roots); on NULL the
# caller falls back to the historical rule, which in the open-gate case
# applies the full correction (the gate implies lambda >= 1 + cutoff -
# 1), so the result is bit-for-bit the historical behavior.
lav_sam_lambda_floor_debias <- function(score = NULL, err = NULL,
lambda = NULL, n = NULL,
floor_default = NULL) {
if (is.null(score) || nrow(score) != n || anyNA(score) ||
n < 4L * ncol(score) || !is.finite(lambda)) {
return(NULL) # unusable: historical rule
}
idx <- seq_len(n)
mod2 <- idx %% 2L
mod4 <- idx %% 4L
half_list <- list(
which(mod2 == 1L), which(mod2 == 0L),
which(mod4 <= 1L), which(mod4 >= 2L),
which(mod4 == 0L | mod4 == 3L), which(mod4 == 1L | mod4 == 2L)
)
lambda_sub <- function(ii) {
nh <- length(ii)
vh <- var(score[ii, , drop = FALSE]) * (nh - 1) / nh
out <- try(lav_mat_sym_diff_smallest_root(vh, err), silent = TRUE)
if (inherits(out, "try-error")) {
return(as.numeric(NA))
}
out
}
lambda_half <- vapply(half_list, lambda_sub, numeric(1L))
lambda_half <- lambda_half[is.finite(lambda_half)]
if (length(lambda_half) < 4L) {
return(NULL) # degenerate: historical rule
}
se_half <- stats::sd(lambda_half) / sqrt(2)
# the gate: comfortably positive sample margin -> historical behavior
if (lambda - 1 >= max(2 * se_half, floor_default)) {
return(NULL)
}
# danger zone: casewise bootstrap for bias and noise of lambda
# (B = 200: with fewer resamples, the Monte Carlo error of se_boot -- and
# hence of the floor, which is a deterministic but arbitrary function of
# the internal seed -- is not negligible: with B = 40, lambda_star varied
# by about +/- 0.1 across seeds in small samples; the cost remains
# negligible, and the floor is frozen in the numerical jacobians)
b_reps <- 200L
seed_old <- NULL
if (exists(".Random.seed", envir = globalenv())) {
seed_old <- get(".Random.seed", envir = globalenv())
}
set.seed(272727L)
lambda_boot <- vapply(seq_len(b_reps), function(b) {
lambda_sub(sample.int(n, n, replace = TRUE))
}, numeric(1L))
if (!is.null(seed_old)) {
assign(".Random.seed", seed_old, envir = globalenv())
}
lambda_boot <- lambda_boot[is.finite(lambda_boot)]
if (length(lambda_boot) < b_reps / 2) {
return(NULL) # degenerate: historical rule
}
bias_hat <- max(0, mean(lambda_boot) - lambda)
se_boot <- stats::sd(lambda_boot)
margin <- min(lambda - 1 + bias_hat, lambda - 1 + 3 * se_boot)
out <- max(margin, 2 * se_boot, floor_default)
# branch bookkeeping for the SE machinery: when the lambda-tracking
# margin binds, lambda_star = lambda - floor is locally constant in
# lambda (kappa = 0); when a lambda-independent guard binds (2*se_boot
# or floor_default), lambda_star tracks lambda one-for-one (kappa = 1)
attr(out, "kappa") <- if (margin >= max(2 * se_boot, floor_default)) {
0
} else {
1
}
out
}
# First-order wrapper: Bartlett factor scores f_i = M (y_i - ybar), so
# that Var(f) equals MSM = M S M' (the per-subsample re-centering of
# var() makes the global centering immaterial). The caller only passes y
# in the supported setting (single-level, continuous, complete data, no
# conditional.x, no dummy/higher-order lvs); everything else -> NULL ->
# historical rule (see lav_sam_lambda_floor_debias() for the design).
lav_sam_veta1_floor_debias <- function(y = NULL, m = NULL, mtm = NULL,
lambda = NULL, n = NULL) {
if (is.null(y) || nrow(y) != n || anyNA(y)) {
return(NULL) # unusable: historical rule
}
yc <- t(t(y) - colMeans(y))
lav_sam_lambda_floor_debias(
score = yc %*% t(m), err = mtm, lambda = lambda, n = n,
floor_default = 1 / (n - 1)
)
}
# compute VETA
# - if alpha_correction == 0 -> same as local SAM (or MOC)
# - if alpha_correction == (N-1) -> same as FSR+Bartlett
lav_sam_veta <- function(m = NULL, s = NULL, mm_theta = NULL,
alpha_correction = 0L, lambda_correction = TRUE,
lambda_floor = "default", y = NULL,
n = 20L, dummy_lv_idx = integer(0L), extra = FALSE) {
# catch empty rows in M (higher-order?)
m_orig <- m
empty_idx <- which(apply(m_orig, 1L, function(x) all(x == 0)))
if (length(empty_idx) > 0L) {
m <- m_orig[-empty_idx, , drop = FALSE]
}
# MSM
msm <- m %*% s %*% t(m)
# MTM
mtm <- m %*% mm_theta %*% t(m)
# note: MTM no longer needs to be regularized (force-pd) here:
# lav_mat_sym_diff_smallest_root() handles a singular (or even zero)
# MTM natively, by regressing out its null space (which includes the
# rows/columns of any dummy lvs, and any empty theta elements)
if (all(abs(mtm) < sqrt(.Machine$double.eps))) {
# all zero? certainly no need for alpha or lambda_correction
alpha_correction <- 0L
lambda_correction <- FALSE
}
# apply small sample correction (if requested)
alpha <- alpha_correction
if (alpha_correction > 0) {
alpha_n1 <- lav_sam_alpha_n1(alpha_correction, n)
mtm <- (1 - alpha_n1) * mtm
}
lambda <- lambda_star <- +Inf
lambda_kappa <- 0
lambda_w <- NULL
lambda_floor_used <- as.numeric(NA)
if (lambda_correction && is.list(lambda_floor)) {
# internal (SE machinery only): frozen effective multiplier -- the
# numeric-jacobian channels re-run this function at perturbed inputs,
# and a kappa = 0 truncation (lambda-tracking debias margin) must keep
# lambda_star fixed at its at-the-estimates value there
lambda_star <- lambda_floor[["star"]]
veta <- msm - lambda_star * mtm
} else if (lambda_correction) {
# use Fuller (1987) approach to ensure VETA is positive
lambda <- try(lav_mat_sym_diff_smallest_root(msm, mtm),
silent = TRUE
)
if (inherits(lambda, "try-error")) {
lav_msg_warn(gettext("failed to compute lambda"))
veta <- msm - mtm # and hope for the best
} else {
# determine the truncation floor:
# NULL -> historical rule (cutoff 1 + 1/(n-1), floor 1/(n-1))
# a number -> coherent rule: truncate iff lambda < 1 + floor
floor_1 <- NULL
kappa_1 <- 1 # a fixed (numeric) floor: lambda_star tracks lambda
if (identical(lambda_floor, "debias") &&
length(empty_idx) == 0L && length(dummy_lv_idx) == 0L) {
# gated bootstrap-debiased floor; NULL when the sample margin is
# comfortably positive (then lambda >= 1 + 1/(n-1) by the gate
# construction, and the historical rule below applies the full
# correction: bit-for-bit the historical behavior) or when the
# setting is ineligible/degenerate (-> historical rule as well)
floor_1 <- lav_sam_veta1_floor_debias(
y = y, m = m, mtm = mtm, lambda = lambda, n = n
)
if (!is.null(floor_1) && !is.null(attr(floor_1, "kappa"))) {
kappa_1 <- attr(floor_1, "kappa")
}
} else if (is.numeric(lambda_floor)) {
floor_1 <- lambda_floor
}
if (!is.null(floor_1)) {
lambda_floor_used <- as.numeric(floor_1)
if (lambda < 1 + floor_1) {
lambda_star <- max(0, lambda - floor_1)
veta <- msm - lambda_star * mtm
# clamped at zero -> VETA = MSM exactly, no lambda term left
lambda_kappa <- if (lambda_star > 0) kappa_1 else 0
} else {
veta <- msm - mtm
}
} else {
# historical rule
lambda_floor_used <- 1 / (n - 1)
cutoff <- 1 + 1 / (n - 1)
if (lambda < cutoff) {
lambda_star <- lambda - 1 / (n - 1)
veta <- msm - lambda_star * mtm
lambda_kappa <- 1 # fixed floor, no clamp: always tracks lambda
} else {
veta <- msm - mtm
}
}
# eigenvector for the SE machinery: when the truncation binds with
# lambda_star tracking lambda (kappa = 1), the jacobian of
# vech(VETA) w.r.t. vech(S) carries the extra rank-1 term
# -vech(MTM) . (d lambda / d vech(S))' with
# d lambda = t(w) dS w, w = t(M) v (v the smallest-root pencil
# eigenvector, normalized t(v) MTM v = 1)
if (extra && lambda_kappa == 1) {
# recompute of the same pencil: suppress its (duplicate) warnings
lam_v <- try(suppressWarnings(
lav_mat_sym_diff_smallest_root(msm, mtm,
vector = TRUE)),
silent = TRUE)
vv <- attr(lam_v, "v")
if (!inherits(lam_v, "try-error") && !is.null(vv) &&
all(is.finite(vv))) {
lambda_w <- drop(t(m) %*% vv)
}
# no usable direction -> lambda_w stays NULL and the analytic
# rank-1 term is simply skipped (old, conservative behavior);
# kappa is kept, so the frozen re-runs still track lambda
}
}
} else {
veta <- msm - mtm
}
# empty.idx?
if (length(empty_idx) > 0L) {
msm_full <- msm
mtm_full <- mtm
veta_full <- veta
nfac_orig <- nrow(m_orig)
msm <- mtm <- veta <- matrix(0, nrow = nfac_orig, ncol = nfac_orig)
msm[-empty_idx, -empty_idx] <- msm_full
mtm[-empty_idx, -empty_idx] <- mtm_full
veta[-empty_idx, -empty_idx] <- veta_full
}
# extra attributes?
if (extra) {
attr(veta, "lambda") <- lambda
attr(veta, "alpha") <- alpha
attr(veta, "lambda_star") <- lambda_star
attr(veta, "MSM") <- msm
attr(veta, "MTM") <- mtm
# SE machinery (see lav_sam_jacb_lambda1_g() and the frozen-floor
# re-runs in the numeric jacobian channels)
attr(veta, "lambda_kappa") <- lambda_kappa
attr(veta, "lambda_w") <- lambda_w
attr(veta, "lambda_floor_used") <- lambda_floor_used
}
veta
}
# compute EETA = E(Eta) = M %*% [YBAR - NU]
lav_sam_eeta <- function(m = NULL, ybar = NULL, mm_nu = NULL) {
m %*% (ybar - mm_nu)
}
# Per-missing-pattern Bartlett factor scores and their measurement-error
# covariances, for the lv-interaction (casewise) machinery under
# missing = "ml"/"ml.x".
#
# For each missing pattern p (observed variables o), the pattern-specific
# mapping matrix M_p is built from the observed sub-matrices (Lambda_o,
# Theta_o), so that M_p Lambda_o = I and hence E(f_i | eta_i) = eta_i for
# every scoreable pattern; the measurement-error covariance of the scores,
# Var(f_i | eta_i) = M_p Theta_o M_p' =: B_p, now varies over the patterns.
# The second-order (VETA2/Gamma) machinery aggregates over the patterns; see
# lav_sam_veta2() for the exact aggregation (B.bar / Q / D).
#
# A pattern is 'scoreable' if every factor keeps at least one observed
# indicator (no all-zero columns in Lambda_o), Lambda_o has full column rank,
# and the mapping matrix could be constructed with no all-zero rows. Cases in
# unscoreable patterns get ok = FALSE and are excluded from the casewise
# machinery -- their information still enters the first-order moments through
# the EM (h1) statistics.
#
# y: raw data (rows = cases, columns = ov in model order), with NAs
# mm_lambda/mm_theta/mm_nu: measurement matrices WITHOUT the lv-interaction
# columns (and in the unstandardized metric; any std.lv rescaling is
# applied by the caller, to scores and B_p alike)
# s: covariance matrix of the ov (EM/h1 based; only used for method = "GLS")
# mp: optional precomputed missing-pattern structure (lav_data_mi_patterns(y))
#
# Returns list(fs = n x nfac score matrix (NA rows where ok = FALSE),
# ok = logical(n), pat_of_case = integer(n) index into bp,
# bp = list of B_p over the scoreable patterns,
# wp = pattern proportions n_p / sum(n_p) over scoreable cases)
lav_sam_fs_missing <- function(y = NULL, mm_lambda = NULL, mm_theta = NULL,
mm_nu = NULL, s = NULL, method = "ML",
mp = NULL) {
n <- nrow(y)
nfac <- ncol(mm_lambda)
if (is.null(mp)) {
mp <- lav_data_mi_patterns(y)
}
fs <- matrix(as.numeric(NA), nrow = n, ncol = nfac)
ok <- rep(FALSE, n)
pat_of_case <- rep(NA_integer_, n)
bp <- list()
wp_n <- integer(0L)
nu <- drop(mm_nu)
for (p in seq_len(mp$npatterns)) {
obs <- mp$pat[p, ]
lam_o <- mm_lambda[obs, , drop = FALSE]
# each factor needs at least one observed indicator, and the observed
# sub-LAMBDA must keep full column rank
if (any(apply(lam_o, 2L, function(x) all(x == 0))) ||
qr(lam_o)$rank < nfac) {
next
}
th_o <- mm_theta[obs, obs, drop = FALSE]
s_o <- NULL
if (!is.null(s)) {
s_o <- s[obs, obs, drop = FALSE]
}
m_p <- tryCatch(
lav_sam_mapping_mat(
mm_lambda = lam_o, mm_theta = th_o,
s = s_o, method = method
),
error = function(e) NULL
)
if (is.null(m_p) || anyNA(m_p) ||
any(apply(m_p, 1L, function(x) all(x == 0)))) {
next
}
case_idx <- mp$case.idx[[p]]
yc_o <- t(t(y[case_idx, obs, drop = FALSE]) - nu[obs])
this_id <- length(bp) + 1L
bp[[this_id]] <- m_p %*% th_o %*% t(m_p)
wp_n <- c(wp_n, length(case_idx))
fs[case_idx, ] <- yc_o %*% t(m_p)
ok[case_idx] <- TRUE
pat_of_case[case_idx] <- this_id
}
list(
fs = fs, ok = ok, pat_of_case = pat_of_case,
bp = bp, wp = wp_n / sum(wp_n)
)
}
# Second-order wrapper: the scores are the (kept columns of the)
# second-order factor-score products fs2, whose covariance is var_fs2;
# the error side of the pencil is var.error. With missing data the
# scores are pattern-based (heterogeneous error contributions), so we
# fall back to the historical rule (see lav_sam_lambda_floor_debias()
# for the design). Note that the second-order sample lambda is much more
# severely biased than its first-order counterpart: the pencil is
# high-dimensional (all latent variables plus all product terms) and
# var_fs2 is built from fourth-order moments of the factor scores.
lav_sam_veta2_floor_debias <- function(fs2 = NULL, var_error = NULL,
lambda = NULL, n = NULL,
mi_flag = FALSE) {
if (mi_flag) {
return(NULL) # pattern-based scores: historical rule
}
lav_sam_lambda_floor_debias(
score = fs2, err = var_error, lambda = lambda, n = n,
floor_default = ncol(fs2) / (n - 1)
)
}
# compute veta including quadratic/interaction terms
lav_sam_veta2 <- function(fs = NULL, m = NULL,
veta = NULL, eeta = NULL, mm_theta = NULL,
lv_names = NULL,
lv_int_names = NULL,
dummy_lv_names = character(0L),
alpha_correction = 0L,
lambda_correction = TRUE,
lambda_floor = "default",
lambda1 = 1,
pattern_list = NULL,
fs_outlier_idx = integer(0L),
return_fs = FALSE,
return_cov_iveta2 = TRUE,
extra = FALSE) {
# small utility function: var() divided by N
varn <- function(x, n) {
var(x, use = "pairwise.complete.obs") * (n - 1) / n
}
if (length(lv_int_names) == 0L) {
lav_msg_stop(gettext("lv.int.names is empty: no lv quadratic/interaction
terms are provided"))
}
if (is.null(lv_names)) {
lv_names <- paste("eta", seq_len(ncol(fs)), sep = "")
}
# missing data (missing = "ml"): the factor scores in fs were computed per
# missing pattern (see lav_sam_fs_missing()), and their measurement-error
# covariance B_p = M_p Theta_p M_p' varies over the patterns. pattern_list
# holds bp (the B_p, in the same metric as fs), wp (pattern proportions)
# and pat_of_case (the pattern of every row of fs). Under MCAR (pattern
# independent of eta), three things change relative to the complete-data
# case:
# (1) every var.error term LINEAR in the error covariance aggregates
# with B.bar = sum_p w_p B_p;
# (2) the quadratic term Var(e* %x% e*) = (I+K)(B %x% B) aggregates as
# (I+K) Q with Q = sum_p w_p (B_p %x% B_p) (NOT B.bar %x% B.bar);
# (3) E_p(f* %x% f*) shifts by vec(B_p) across the patterns, so the raw
# Var(fs2) picks up a between-pattern term
# D = sum_p w_p vec(B_p) vec(B_p)' - vec(B.bar) vec(B.bar)'
# that must be part of var.error as well.
# With a single pattern (B_p = B.bar, Q = B %x% B, D = 0) all of this
# reduces exactly to the complete-data expressions below.
mi_flag <- !is.null(pattern_list)
# MTM (B.bar over the patterns when mi_flag)
if (mi_flag) {
bp_aug <- lapply(pattern_list$bp, function(b) lav_mat_bdiag(0, b))
wp <- pattern_list$wp
pat_of_case <- pattern_list$pat_of_case
mtm <- Reduce(`+`, Map(function(b, w) w * b, bp_aug, wp))
} else {
mtm <- m %*% mm_theta %*% t(m)
}
# note: MTM no longer needs to be regularized (force-pd) here:
# lav_mat_sym_diff_smallest_root() handles a singular var.error matrix
# natively, by regressing out its null space (which includes the
# rows/columns of any dummy lvs); this also keeps MTM identical to the
# (raw) version used in lav_sam_gamma_add()
# augment to include intercept
fs <- cbind(1, fs)
n <- nrow(fs)
if (!mi_flag) {
mtm <- lav_mat_bdiag(0, mtm)
}
veta <- lav_mat_bdiag(0, veta)
eeta <- c(1, eeta)
lv_names <- c("..int..", lv_names)
nfac <- ncol(fs)
tmp <- lav_sam_lvnames2(lv_names)
idx1 <- tmp$idx1
idx2 <- tmp$idx2
names_1 <- tmp$names
fs2 <- fs[, idx1, drop = FALSE] * fs[, idx2, drop = FALSE]
k_nfac <- lav_mat_com(nfac, nfac)
ik <- diag(nfac * nfac) + k_nfac
eeta <- as.matrix(drop(eeta))
vetak_mtm <- veta %x% mtm
# normal version (for now):
if (mi_flag) {
# quadratic error moment: (I+K) Q with Q = sum_p w_p (B_p %x% B_p)
qfull <- Reduce(`+`, Map(function(b, w) w * (b %x% b), bp_aug, wp))
gamma_me22 <- ik %*% qfull
# between-pattern term D: vec(B_p)[a] in the (idx1, idx2) enumeration of
# (eta* %x% eta*) is B_p[idx1[a], idx2[a]]
vb <- vapply(bp_aug, function(b) b[cbind(idx1, idx2)],
numeric(length(idx1)))
vbbar <- mtm[cbind(idx1, idx2)]
dfull <- vb %*% (wp * t(vb)) - tcrossprod(vbbar)
} else {
gamma_me22 <- ik %*% (mtm %x% mtm)
dfull <- 0
}
# ingredients (normal ME case)
var_fs2 <- varn(fs2, n)
var_etak_me <- (tcrossprod(eeta) %x% mtm + vetak_mtm)
var_mek_eta <- lav_mat_com_pre_post(var_etak_me)
var_me2 <- gamma_me22
cov_etak_me_mek_eta <- lav_mat_com_post(var_etak_me)
cov_mek_eta_etak_me <- t(cov_etak_me_mek_eta)
var_error <- (var_etak_me + var_mek_eta + cov_etak_me_mek_eta
+ cov_mek_eta_etak_me + var_me2 + dfull)
# select only what we need
colnames(var_fs2) <- rownames(var_fs2) <- names_1
colnames(var_error) <- rownames(var_error) <- names_1
colnames(fs2) <- names_1
lv_keep <- c(lv_names[-1], lv_int_names)
var_fs2 <- var_fs2[lv_keep, lv_keep, drop = FALSE]
var_error <- var_error[lv_keep, lv_keep, drop = FALSE]
fs_mean <- colMeans(fs2, na.rm = TRUE)
names(fs_mean) <- names_1
fs2_mean <- fs_mean # all of them
fs_mean <- fs_mean[lv_keep]
# apply small sample correction (if requested)
alpha <- alpha_correction
alpha_n1 <- 0
if (alpha_correction > 0) {
alpha_n1 <- lav_sam_alpha_n1(alpha_correction, n)
var_error <- (1 - alpha_n1) * var_error
}
lambda <- +Inf
lambda_star <- 1
lambda_kappa <- 0
lambda_v <- NULL
lambda_floor_used <- as.numeric(NA)
if (lambda_correction && is.list(lambda_floor)) {
# internal (SE machinery only): frozen effective multiplier -- the
# numeric-jacobian channels re-run this function at perturbed inputs,
# and a kappa = 0 truncation (lambda-tracking debias margin, or a
# clamp at zero) must keep lambda_star fixed at its at-the-estimates
# value there
lambda_star <- lambda_floor[["star"]]
veta2 <- var_fs2 - lambda_star * var_error
} else if (lambda_correction) {
# use Fuller (1987) approach to ensure VETA2 is positive
lambda <- try(lav_mat_sym_diff_smallest_root(
var_fs2,
var_error
), silent = TRUE)
if (inherits(lambda, "try-error")) {
lav_msg_warn(gettext("failed to compute lambda"))
veta2 <- var_fs2 - var_error # and hope for the best
} else {
# determine the truncation floor:
# NULL -> historical rule (cutoff 1 + 2/n, floor p2/(n-1))
# a number -> coherent rule: truncate iff lambda < 1 + floor
floor_2 <- NULL
kappa_2 <- 1 # a fixed (numeric) floor: lambda_star tracks lambda
if (identical(lambda_floor, "debias")) {
# gated bootstrap-debiased floor; NULL when the sample margin is
# comfortably positive (then lambda >= 1 + 2/n by the gate
# construction, and the historical rule below applies the full
# correction: bit-for-bit the historical behavior) or when the
# setting is ineligible/degenerate (-> historical rule as well)
floor_2 <- lav_sam_veta2_floor_debias(
fs2 = fs2[, lv_keep, drop = FALSE],
var_error = var_error, lambda = lambda, n = n,
mi_flag = mi_flag
)
if (!is.null(floor_2) && !is.null(attr(floor_2, "kappa"))) {
kappa_2 <- attr(floor_2, "kappa")
}
} else if (is.numeric(lambda_floor)) {
floor_2 <- lambda_floor
}
if (!is.null(floor_2)) {
lambda_floor_used <- as.numeric(floor_2)
if (lambda < 1 + floor_2) {
lambda_star <- max(c(0, lambda - floor_2))
veta2 <- var_fs2 - lambda_star * var_error
# clamped at zero -> VETA2 = var_fs2 exactly, no lambda term
lambda_kappa <- if (lambda_star > 0) kappa_2 else 0
} else {
veta2 <- var_fs2 - var_error
}
} else {
# historical rule; note that the p2/(n-1) truncation margin
# vanishes with n, while the sample lambda of this
# high-dimensional fourth-order-moment pencil is severely
# downward-biased -- in weakly-identified conditions the
# truncation then binds over a long pre-asymptotic range of n
# with an ever smaller margin, and the accuracy of the
# structural estimates can *deteriorate* with increasing n;
# lambda2_floor = "debias" is a data-driven alternative
lambda_floor_used <- ncol(var_fs2) / (n - 1)
cutoff <- 1 + 2 / n # be more conservative for VETA2
if (lambda < cutoff) {
lambda_star <- max(c(0, lambda - ncol(var_fs2) / (n - 1)))
veta2 <- var_fs2 - lambda_star * var_error
lambda_kappa <- if (lambda_star > 0) 1 else 0
} else {
veta2 <- var_fs2 - var_error
}
}
# pencil eigenvector for the SE machinery: when the truncation
# binds with lambda_star tracking lambda (kappa = 1), the casewise
# contributions below gain the influence of lambda-hat itself
# (recompute of the same pencil: suppress its duplicate warnings)
if ((return_cov_iveta2 || extra) && lambda_kappa == 1) {
lam_v <- try(suppressWarnings(
lav_mat_sym_diff_smallest_root(var_fs2, var_error,
vector = TRUE)),
silent = TRUE)
vv <- attr(lam_v, "v")
if (!inherits(lam_v, "try-error") && !is.null(vv) &&
all(is.finite(vv))) {
lambda_v <- vv
}
# no usable direction -> the IF term is skipped (old behavior)
}
}
} else {
veta2 <- var_fs2 - var_error
}
# new in 0.6-20: to compute Gamma.eta
# the casewise contributions are defined so that their average over the
# observations equals *exactly* the augmented summary statistics
# (EETA2, VETA2) that are used in the second step; this requires the
# effective *first-order* lambda1 inside the decomposition (because
# EETA2 and var.error are built from the lambda1-corrected first-order
# VETA), while the *second-order* lambda_star multiplies the whole
# correction term:
# ieeta2_i = fs2[i, ] - lambda1 * vec(MTM)
# (the terms involving the mean factor scores cancel out exactly)
# and
# iveta2_i = tcrossprod(fs2[i, ] - fs2.mean) - lambda2.eff * tmp_i
# where lambda2.eff = lambda_star * (1 - alpha/(n-1)) is the effective
# multiplier of the *unscaled* var.error (the alpha correction rescales
# var.error before the lambda step), and
# tmp_i = ((F_i - lambda1 * MTM) %x% MTM) +
# (MTM %x% (F_i - lambda1 * MTM)) +
# lav_mat_com_post((F_i - lambda1 * MTM) %x% MTM) +
# lav_mat_com_pre((F_i - lambda1 * MTM) %x% MTM) +
# (I + K) %*% (MTM %x% MTM)
# with F_i = tcrossprod(f_i); both selected in the rows/columns we keep
# note: cov(iveta2_1) is invariant to the value of lambda1 (the lambda1
# terms are constant across the observations); but lambda1 does matter
# for the jacobian of the average casewise contribution computed in
# lav_sam_gamma_add()
# rewritten 12 June 2026 to avoid the loop over the observations: every
# kept element is a linear function of the second-order factor scores
# fs2[i, ]; elementwise, for kept elements a = (r1, s1) and b = (r2, s2)
# of (eta* %x% eta*):
# tmp_i[a, b] = F_i[r1, r2] MTM[s1, s2] + MTM[r1, r2] F_i[s1, s2] +
# F_i[r1, s2] MTM[s1, r2] + F_i[s1, r2] MTM[r1, s2] +
# (1 - 2 * lambda1) * (MTM[r1, r2] MTM[s1, s2] +
# MTM[s1, r2] MTM[r1, s2])
# where F_i[r, s] = fs2[i, (r - 1) * nfac + s]
if (return_cov_iveta2) {
lambda2_eff <- lambda_star * (1 - alpha_n1)
keep_idx <- match(lv_keep, names_1)
nkeep <- length(keep_idx)
# row/column pairs of the kept elements of (eta* %x% eta*)
i1k <- idx1[keep_idx]
i2k <- idx2[keep_idx]
# casewise contributions to E(eta* %x% eta*)[keep]
if (mi_flag) {
# pattern-specific error-mean shift: vec(B_p) replaces vec(B.bar), so
# that the between-pattern shift of E_p(f* %x% f*) does not inflate
# the covariance of the contributions
part1 <- fs2[, keep_idx, drop = FALSE] -
lambda1 * t(vb)[pat_of_case, keep_idx, drop = FALSE]
} else {
part1 <- t(t(fs2[, keep_idx, drop = FALSE]) -
lambda1 * lav_mat_vec(mtm)[keep_idx])
}
# casewise contributions to vech(Var(eta* %x% eta*)[keep, keep])
# vech row (pa) and column (pb) indices over the kept elements
tmp_mat <- matrix(seq_len(nkeep), nkeep, nkeep)
pa <- tmp_mat[lower.tri(tmp_mat, diag = TRUE)]
pb <- t(tmp_mat)[lower.tri(tmp_mat, diag = TRUE)]
r1 <- i1k[pa]; s1 <- i2k[pa]
r2 <- i1k[pb]; s2 <- i2k[pb]
# fs2 columns holding F_i[r1, r2], F_i[s1, s2], F_i[r1, s2], F_i[s1, r2]
col_11 <- (r1 - 1L) * nfac + r2
col_22 <- (s1 - 1L) * nfac + s2
col_12 <- (r1 - 1L) * nfac + s2
col_21 <- (s1 - 1L) * nfac + r2
b_11 <- mtm[cbind(r1, r2)]
b_22 <- mtm[cbind(s1, s2)]
b_12 <- mtm[cbind(r1, s2)]
b_21 <- mtm[cbind(s1, r2)]
fs2kc <- t(t(fs2[, keep_idx, drop = FALSE]) - fs2_mean[keep_idx])
if (mi_flag) {
# pattern-heterogeneous quadratic terms. The F-cross terms (with B.bar)
# average to the E-terms PLUS 2*lambda1*(bb1 + bb2); the -2*lambda1
# constant cancels this, and the per-pattern products average to the
# (I+K)Q quadratic part (q12) and the between-pattern D part (dvec),
# so the average over the observations equals the aggregated var.error
# above. With a single pattern this reduces exactly to the
# (1 - 2*lambda1) * (bb1 + bb2) constant of the complete-data branch.
q12_p <- vapply(bp_aug, function(b)
b[cbind(r1, r2)] * b[cbind(s1, s2)] +
b[cbind(s1, r2)] * b[cbind(r1, s2)], numeric(length(r1)))
dvec_p <- vapply(bp_aug, function(b)
b[cbind(r1, s1)] * b[cbind(r2, s2)], numeric(length(r1)))
const_bb <- -2 * lambda1 * (b_11 * b_22 + b_21 * b_12) -
mtm[cbind(r1, s1)] * mtm[cbind(r2, s2)]
quad_case <- t(q12_p + dvec_p)[pat_of_case, , drop = FALSE] +
rep(const_bb, each = n)
} else {
quad_case <- rep((1 - 2 * lambda1) * (b_11 * b_22 + b_21 * b_12),
each = n)
}
part2 <- (fs2kc[, pa, drop = FALSE] * fs2kc[, pb, drop = FALSE] -
lambda2_eff *
(t(t(fs2[, col_11, drop = FALSE]) * b_22) +
t(t(fs2[, col_22, drop = FALSE]) * b_11) +
t(t(fs2[, col_12, drop = FALSE]) * b_21) +
t(t(fs2[, col_21, drop = FALSE]) * b_12) +
quad_case))
# sampling variability of the second-order lambda-hat itself
# (kappa = 1 truncation): d lambda = t(v) (d var_fs2 -
# lambda * d var_error) v, with v the smallest-root pencil
# eigenvector (t(v) var_error v = 1). Casewise:
# IF_i = [ (z_i'v)^2 - mean ] -
# lambda * (1 - alpha_n1) * [ u_i - mean(u) ]
# where z_i are the centered kept second-order scores, and
# u_i = t(v) tmp_i v is the casewise var.error representation: with
# V the nfac x nfac embedding of v (V[i1k, i2k] = v) and
# Vs = V + t(V), the F-linear part collapses into the single form
# u_i = < F_i, Vs B Vs >
# (all four F x B index shapes of tmp_i; the case-constant parts
# cancel in the centering; under mi the per-pattern quadratic parts
# are added per case). The vech(VETA2) contributions then gain
# - var_error[pa, pb] * IF_i, which is exactly mean-zero: the
# statistics are unchanged, only their covariance changes.
if (lambda_kappa == 1 && !is.null(lambda_v)) {
zv <- drop(fs2kc %*% lambda_v)
v2mat <- matrix(0, nrow = nfac, ncol = nfac)
v2mat[cbind(i1k, i2k)] <- lambda_v
vs_m <- v2mat + t(v2mat)
w_m <- vs_m %*% mtm %*% vs_m
u_i <- drop(fs2 %*% as.vector(t(w_m)))
if (mi_flag) {
# per-pattern quadratic part (the (I+K)Q and D terms vary over
# the patterns; globally constant parts cancel in the centering)
u_quad <- vapply(bp_aug, function(b) {
q1 <- b[i1k, i1k, drop = FALSE] * b[i2k, i2k, drop = FALSE]
q2 <- b[i2k, i1k, drop = FALSE] * b[i1k, i2k, drop = FALSE]
drop(crossprod(lambda_v, (q1 + q2) %*% lambda_v)) +
sum(lambda_v * b[cbind(i1k, i2k)])^2
}, numeric(1L))
u_i <- u_i + u_quad[pat_of_case]
}
if_lambda <- (zv * zv - mean(zv * zv)) -
lambda * (1 - alpha_n1) * (u_i - mean(u_i))
part2 <- part2 - tcrossprod(if_lambda, var_error[cbind(pa, pb)])
}
iveta2_1 <- cbind(part1, part2)
cov_iveta2 <- cov(iveta2_1) * (n - 1) / n
}
# extra attributes?
if (extra) {
attr(veta2, "lambda") <- lambda
attr(veta2, "alpha") <- alpha
attr(veta2, "lambda_star") <- lambda_star
attr(veta2, "MSM") <- var_fs2
attr(veta2, "MTM") <- var_error
attr(veta2, "FS.mean") <- fs_mean
# SE machinery (see lav_sam_gamma_add() and the frozen-floor re-runs
# in the numeric jacobian channels)
attr(veta2, "lambda_kappa") <- lambda_kappa
attr(veta2, "lambda_floor_used") <- lambda_floor_used
}
if (return_fs) {
attr(veta2, "FS") <- fs2[, lv_keep, drop = FALSE]
}
if (return_cov_iveta2) {
attr(veta2, "cov.iveta2") <- cov_iveta2
# the casewise contributions themselves (needed for the casewise
# Gamma.eta, see lav_sam_gamma_add())
attr(veta2, "iveta2") <- iveta2_1
}
veta2
}
lav_sam_eeta2 <- function(eeta = NULL, veta = NULL, lv_names = NULL,
lv_int_names = NULL) {
if (length(lv_int_names) == 0L) {
lav_msg_stop(gettext("lv.int.names is empty: no lv quadratic/interaction
terms are provided"))
}
if (is.null(lv_names)) {
lv_names <- paste("eta", seq_len(ncol(veta)), sep = "")
}
nfac <- nrow(veta)
idx1 <- rep(seq_len(nfac), each = nfac)
idx2 <- rep(seq_len(nfac), times = nfac)
names_1 <- c(lv_names, paste(lv_names[idx1], lv_names[idx2], sep = ":"))
# E(\eta %x% \eta)
eeta2 <- lav_mat_vec(veta) + eeta %x% eeta
# add 1st order
eeta2_aug <- c(eeta, eeta2)
# select only what we need
names(eeta2_aug) <- names_1
lv_keep <- c(lv_names, lv_int_names)
eeta2_aug <- eeta2_aug[lv_keep]
eeta2_aug
}
# compute var(fs2) including quadratic/interaction terms
lav_sam_fs2 <- function(fs = NULL, lv_names = NULL, lv_int_names = NULL) {
varn <- function(x, n) {
var(x) * (n - 1) / n
}
if (length(lv_int_names) == 0L) {
lav_msg_stop(gettext("lv.int.names is empty: no lv quadratic/interaction
terms are provided"))
}
if (is.null(lv_names)) {
lv_names <- paste("eta", seq_len(ncol(fs)), sep = "")
}
# augment to include intercept
fs <- cbind(1, fs)
n <- nrow(fs)
lv_names <- c("..int..", lv_names)
nfac <- ncol(fs)
idx1 <- rep(seq_len(nfac), each = nfac)
idx2 <- rep(seq_len(nfac), times = nfac)
names_1 <- paste(lv_names[idx1], lv_names[idx2], sep = ":")
fs2 <- fs[, idx1, drop = FALSE] * fs[, idx2, drop = FALSE]
var_fs2 <- varn(fs2, n)
# select only what we need
colnames(var_fs2) <- rownames(var_fs2) <- names_1
lv_main <- paste(lv_names[-1], "..int..", sep = ":")
lv_keep <- c(lv_main, lv_int_names)
var_fs2 <- var_fs2[lv_keep, lv_keep, drop = FALSE]
var_fs2
}
# create consistent lavaan object, based on (filled in) PT
lav_sam_step3_joint <- function(fit = NULL, pt_1 = NULL, sam_method = "local") {
lavoptions <- fit@Options
lavoptions_joint <- lavoptions
lavoptions_joint$optim.method <- "none"
lavoptions_joint$optim.parscale <- "none"
lavoptions_joint$start <- "default"
lavoptions_joint$optim.force.converged <- TRUE
lavoptions_joint$check.gradient <- FALSE
lavoptions_joint$check.start <- FALSE
lavoptions_joint$check.post <- FALSE
lavoptions_joint$rotation <- "none"
lavoptions_joint$se <- "none"
lavoptions_joint$store.vcov <- FALSE # we do this manually
if (sam_method %in% c("local", "fsr", "cfsr")) {
lavoptions_joint$baseline <- FALSE
lavoptions_joint$fit.by.level <- FALSE
lavoptions_joint$sample.icov <- FALSE
#lavoptions.joint$h1 <- TRUE # we need this if we re-use the sam object
lavoptions_joint$test <- "none"
lavoptions_joint$estimator <- "none"
} else {
# global: the joint computes the core (base) test(s); "yuan.chan" is not a
# core lavaan test (it is added afterwards by lav_sam_global_test()), so we
# strip it here and make sure the unscaled "standard" base is present.
joint_test <- lavoptions$test[lavoptions$test != "yuan.chan"]
if (length(joint_test) == 0L || !"standard" %in% joint_test) {
joint_test <- unique(c("standard", joint_test))
}
lavoptions_joint$test <- joint_test
lavoptions_joint$estimator <- lavoptions$estimator
}
# set ustart values
pt_1$ustart <- pt_1$est # as this is used if optim.method == "none"
# the step-0 dummy fit of the local methods skips the inverse of S
# (sample.icov = FALSE: the inverse may not exist when N < nvar), but the
# joint-model computations of some estimators need S^-1 (eg the GLS
# information/objective). Backfill it here when it is missing; if S cannot
# be inverted it simply stays missing, exactly as before.
lavsamplestats_joint <- fit@SampleStats
if (fit@Data@nlevels == 1L && !fit@Model@categorical &&
!lavsamplestats_joint@missing.flag) { # core lavaan skips icov for
# missing = "ml" on purpose
for (g in seq_len(lavsamplestats_joint@ngroups)) {
if (is.null(lavsamplestats_joint@icov[[g]]) &&
!is.null(lavsamplestats_joint@cov[[g]])) {
# quietly: the inverse is optional here (only some estimators need
# it), so a non-PD S must not surface a warning at this point
current_warn <- lav_warn()
lav_warn(FALSE)
out <- try(lav_samp_icov(
cov_1 = lavsamplestats_joint@cov[[g]], ridge = 1e-05,
x_idx = lavsamplestats_joint@x.idx[[g]],
ngroups = fit@Data@ngroups, g = g
), silent = TRUE)
lav_warn(current_warn)
if (!inherits(out, "try-error")) {
lavsamplestats_joint@icov[[g]] <- out$icov
lavsamplestats_joint@cov.log.det[[g]] <- out$cov.log.det
}
}
}
}
joint <- lavaan::lavaan(pt_1,
slot_options = lavoptions_joint,
slot_sample_stats = lavsamplestats_joint,
slot_data = fit@Data,
verbose = FALSE
)
joint
}
# the (two-step robust) STRUCTURAL fit object of a local SAM model.
#
# The structural model FIT.PA is fitted with se = "robust.sem" using the
# step-1 induced NACOV (Gamma.eta) of the latent (co)variances/means. Its
# Satorra-Bentler scaled chi-square and the robust fit measures derived from it
# are the correct 'two-step robust' fit measures.
#
# For a single group, and for multiple groups WITHOUT across-group constraints
# in the measurement model (Case A), the per-group Gamma.eta list passed to
# FIT.PA already yields the correct robust test, and we simply reproduce it.
# For multiple groups WITH across-group measurement constraints (Case B), the
# step-1 variability couples the groups, so the relevant Gamma has nonzero
# cross-group blocks (step1$Gamma.eta.full) that a per-group NACOV list cannot
# represent. We then recompute the Satorra-Bentler scaled chi-square ourselves
# -- for the model AND the baseline -- using that full Gamma.
#
# The corrected FIT.PA is stored in the @internal slot of the final sam
# object (element sam.struc.fit.object), so that fitMeasures(), lavResiduals()
# and modindices() can delegate to it (issue #517).
lav_sam_struc_fit_object <- function(fit_pa = NULL, step1 = NULL) {
# do we have a (local) Gamma.eta and a non-saturated structural model?
have_gamma <- !is.null(step1$Gamma.eta) &&
length(step1$Gamma.eta) > 0L && !is.null(step1$Gamma.eta[[1]])
df_struc <- if (length(fit_pa@test) > 0L) fit_pa@test[[1]]$df else 0L
has_baseline <- !is.null(fit_pa@baseline) &&
!is.null(fit_pa@baseline$partable)
if (!have_gamma || df_struc < 1L || !has_baseline) {
return(fit_pa)
}
out <- try(
{
fp <- fit_pa
ng <- fp@Data@ngroups
ntot <- fp@SampleStats@ntotal
fg <- unlist(fp@SampleStats@nobs) / ntot
meanstr <- fp@Model@meanstructure
# full Gamma of the stacked structural statistics (cross-group blocks
# for Case B; block-diagonal otherwise)
if (isTRUE(step1$caseB)) {
gamma_all <- ntot * step1$Gamma.eta.full
} else {
gamma_all <- lav_mat_bdiag(lapply(
seq_len(ng), function(g) step1$Gamma.eta[[g]] / fg[g]
))
}
# model: Satorra-Bentler scaled test using the full Gamma
test_model <- lav_test_sb(
lavobject = fp, test = "satorra.bentler",
gamma_full = gamma_all
)
# baseline: refit the independence model on the latent (co)variances,
# then apply the same Satorra-Bentler correction with the full Gamma.
# (The corrected adjusted tests, if any, were added afterwards by
# lav_sam_step2_corrected_test() -- they need the casewise rows of
# Gamma.eta, which this moments-only refit does not have.)
opts_base <- fp@Options
opts_base$test <- setdiff(opts_base$test, lav_sam_corrected_tests)
if (length(opts_base$test) == 0L) {
opts_base$test <- "standard"
}
fit_base <- lavaan::lavaan(
model = fp@baseline$partable,
sample.cov = step1$VETA,
sample.mean = if (meanstr) step1$EETA else NULL,
sample.nobs = as.list(unlist(fp@SampleStats@nobs)),
nacov = step1$Gamma.eta,
slot_options = opts_base
)
test_base <- lav_test_sb(
lavobject = fit_base, test = "satorra.bentler",
gamma_full = gamma_all
)
# inject the corrected scaled tests
fp@test[["satorra.bentler"]] <- test_model[["satorra.bentler"]]
base_idx <- which(vapply(fp@baseline$test,
function(x) x$test, character(1L)) == "satorra.bentler")
if (length(base_idx) > 0L) {
fp@baseline$test[[base_idx[1]]] <- test_base[["satorra.bentler"]]
}
fp@external$sam.struc.corrected <- TRUE
fp
},
silent = TRUE
)
if (inherits(out, "try-error")) {
return(fit_pa) # the untouched original
}
out
}
# (two-step robust) fit measures for the structural part of a local SAM model,
# extracted from the (corrected) structural fit object
lav_sam_struc_fit <- function(fit_pa = NULL) {
measures <- c("chisq", "df", "pvalue", "cfi", "rmsea", "srmr")
fallback <- function() {
out <- try(fitMeasures(fit_pa, measures), silent = TRUE)
if (inherits(out, "try-error")) {
out <- "(unable to obtain fit measures)"
names(out) <- "warning"
}
out
}
if (!isTRUE(fit_pa@external$sam.struc.corrected)) {
return(fallback())
}
# use the fitMeasures() labels for the SCALED/ROBUST measures, so the
# printed header makes clear these are the two-step corrected values
# (chisq.scaled / pvalue.scaled / cfi.robust / rmsea.robust), exactly as
# fitMeasures() reports them. (The naive fallback() above keeps the plain
# chisq/pvalue/cfi/rmsea labels.)
scaled_labels <- c(
"chisq.scaled", "df.scaled", "pvalue.scaled",
"cfi.robust", "rmsea.robust", "srmr"
)
fm <- try(fitMeasures(fit_pa, scaled_labels), silent = TRUE)
if (inherits(fm, "try-error")) {
return(fallback())
}
# measures that fitMeasures() does not provide for this scaled test (eg
# the robust cfi/rmsea when the scaled test is a mean-var-adjusted one)
# are reported as NA
out <- rep(as.numeric(NA), length(scaled_labels))
names(out) <- scaled_labels
out[names(fm)] <- as.numeric(fm)
out
}
# Yuan & Chan (2002) rescaled GLOBAL test statistic for sam.method = "global".
#
# The 'joint' object already carries the base statistic T = n F[S, Sigma(theta,
# gamma.hat)] -- the discrepancy of the FULL (observed-variable) model evaluated
# at the plugged-in SAM estimates -- and the correct reference df = p* - q_total
# (all measurement AND structural parameters counted) in joint@test[[1]]. What is
# missing is the scaling correction tau.hat = tr(Q Gamma)/(p* - q) of eq. (22),
# which accounts for (a) the non-normality of the data (via the ADF Gamma) and
# (b) the fact that the measurement parameters gamma were estimated SEPARATELY in
# step 1 (segregated estimation), not jointly with theta.
#
# Key simplification that lets us reuse lav_test_sb() verbatim. With
# H.theta = sigma.dot.theta (sigma.dot.theta' W sigma.dot.theta)^-1
# sigma.dot.theta' W (the W-projector onto the
# structural columns of Delta),
# D = I - sigma.dot.gamma P (P = d gamma.hat / d stats, eq. 12a),
# Yuan & Chan's U (eq. 19) gives I - sigma.dot U = (I - H.theta) D, hence
# Q = (I - sigma.dot U)' W (I - sigma.dot U) = D' U.struc D ,
# with U.struc = W - W Delta.theta (Delta.theta' W Delta.theta)^-1 Delta.theta' W
# -- exactly the Satorra-Bentler residual-weight matrix, but built from the
# STRUCTURAL columns of Delta only. By the cyclic trace property
# tr(Q Gamma) = tr( U.struc . (D Gamma D') ).
# So we hand lav_test_sb():
# - delta = Delta.theta (structural columns of the full Delta)
# - e_inv = (sum_g f_g Delta.theta_g' W_g Delta.theta_g)^-1
# - wls_v = W (the estimator weight of the full model)
# - gamma_full = D Gamma D' (the 'measurement-residualized' Gamma)
# - test_unscaled = the standard (unscaled) full-model test (T and df)
# and read off the satorra.bentler-scaled statistic, which IS eq. (22).
#
# P (the step-1 measurement influence) is the SAME jacobian already used by the
# twostep.robust SEs (lav_sam_step1_local_jac(p_only = TRUE)). Multigroup support
# (Case A block-diagonal, Case B cross-group) follows from the stacked assembly:
# D and Gamma.tilde are built in the stacked statistics space, so the cross-group
# coupling that across-group measurement constraints induce in P is carried
# through automatically.
#
# Returns a list with two components, to be stored in the joint object:
# $test -> joint@test (standard + the YC-scaled test)
# $baseline.test -> joint@baseline$test (standard + a scaled baseline test)
# The baseline (independence) model has no measurement model, so the Yuan-Chan
# correction there reduces to the ordinary Satorra-Bentler scaling; we obtain it
# by refitting the baseline with test = "satorra.bentler". Both scaled tests are
# needed so that the robust/scaled fit measures (CFI, RMSEA) can be computed. On
# any failure the original (unscaled) test list is returned unchanged.
# 'test' is the internal SB-family base used to drive lav_test_sb (eq. 22 is the
# satorra.bentler-type rescaling); the resulting entry is renamed to "yuan.chan"
# (and given a Yuan-Chan label) in both the model and the baseline test lists.
lav_sam_global_test <- function(joint = NULL, step1 = NULL, step2 = NULL,
fit = NULL, test = "satorra.bentler") {
test_orig <- joint@test
# rename a satorra.bentler entry (and its $test field) to "yuan.chan"
rename_yc <- function(test_list) {
idx <- which(names(test_list) == "satorra.bentler")
if (length(idx) == 0L) {
return(test_list)
}
e <- test_list[[idx[1]]]
e$test <- "yuan.chan"
e$label <- "Yuan-Chan (2002) correction"
names(test_list)[idx[1]] <- "yuan.chan"
test_list[[idx[1]]] <- e
test_list
}
fallback <- list(test = test_orig, baseline.test = NULL)
# conditional.x: not supported (yet). The machinery below runs, but with
# NON-NORMAL (eg skewed or binary) covariates the scaled statistic is badly
# anti-conservative (empirical rejection rates up to ~0.3 at the 5% level
# in simulations, with either random or fixed covariate values; the
# simultaneous sem() satorra.bentler test on the same data is accurate, and
# so is the sam yuan.chan test with multivariate normal covariates -- the
# interplay between the unconditionally fitted step-1 blocks and the
# conditional full-model discrepancy is not yet correctly captured by the
# D Gamma D' residualization). Report the honest unscaled test instead.
# (The corrected STRUCTURAL test of sam.method = "local" does not have
# this problem: it is accurate in the same simulations.)
if (joint@Model@conditional.x) {
lav_msg_warn(gettext(
"the Yuan & Chan (2002) scaled test statistic (test = \"yuan.chan\") is
not available (yet) when conditional_x = TRUE; the standard (unscaled)
test is reported instead."))
return(fallback)
}
# locate the standard (unscaled) full-model test = base statistic T
std_idx <- which(vapply(test_orig, function(x) x$test, character(1L)) ==
"standard")
if (length(std_idx) == 0L) {
return(fallback) # no usable base statistic (eg test = "none")
}
std <- test_orig[[std_idx[1]]]
if (is.null(std$df) || anyNA(std$df) || std$df < 1L) {
return(fallback) # saturated full model: no global test
}
out <- try(
{
ngroups <- joint@Data@ngroups
ntot <- joint@SampleStats@ntotal
fg <- unlist(joint@SampleStats@nobs) / ntot
step1_free_idx <- step1$step1.free.idx
step2_free_idx <- step2$step2.free.idx
# ingredients of the full (joint) model:
# stored NACOV, else recompute with the fit-time settings. The
# (unbiased) ADF Gamma is unavailable for the fixed.x / conditional.x
# setting (lav_samp_gamma() with unbiased = TRUE stops), so fall back
# to the *biased* ADF Gamma there (gamma.unbiased = FALSE) -- exactly
# the Gamma that simultaneous sem(test = "satorra.bentler") uses for
# such models. Asymptotically equivalent; the simple (no-x) cases keep
# the unbiased Gamma so their values are unchanged.
gamma <- tryCatch(lav_gamma_used(joint), error = function(e) NULL)
if (is.null(gamma) || is.null(gamma[[1]])) {
opts <- joint@Options
opts$gamma.unbiased <- FALSE
gamma <- lav_gamma_used(
lavdata = joint@Data,
lavoptions = opts,
lavsamplestats = joint@SampleStats,
lavh1 = joint@h1,
lavimplied = joint@implied
)
}
delta <- lavTech(joint, "Delta")
wls_v <- lavTech(joint, "WLS.V")
# P = d(gamma.mm)/d(stats), reusing the twostep.robust jacobian; align its
# rows with step1.free.idx (P rows are in ascending free-parameter order)
p <- lav_sam_step1_local_jac(step1 = step1, fit = fit, p_only = TRUE)
p_row_idx <- match(step1_free_idx, attr(p, "free.idx"))
if (anyNA(p_row_idx)) {
lav_msg_stop(gettext(
"internal error: unable to align the step 1 jacobian (P) for the
global test statistic"))
}
p <- p[p_row_idx, , drop = FALSE]
# per group: structural and measurement columns of Delta, the weight W,
# and the SB-convention Gamma (Cov of sqrt(N) * per-group statistics)
delta_theta_list <- vector("list", ngroups)
delta_gamma_list <- vector("list", ngroups)
gamma_sb_list <- vector("list", ngroups)
a_theta <- matrix(0, length(step2_free_idx), length(step2_free_idx))
for (g in seq_len(ngroups)) {
delta_theta_list[[g]] <- delta[[g]][, step2_free_idx, drop = FALSE]
delta_gamma_list[[g]] <- delta[[g]][, step1_free_idx, drop = FALSE]
gamma_sb_list[[g]] <- gamma[[g]] / fg[g]
# structural information E.theta = sum_g f_g Delta.theta_g' W_g Delta.theta_g
wd <- lav_sam_wls_delta(wls_v[[g]], delta_theta_list[[g]])
a_theta <- a_theta + fg[g] * crossprod(delta_theta_list[[g]], wd)
}
e_theta_inv <- solve(a_theta)
# D (stacked) = I - Delta.gamma P ; Gamma.tilde = D Gamma D'
delta_gamma_all <- do.call(rbind, delta_gamma_list)
nstat_all <- nrow(delta_gamma_all)
d_mat <- diag(nstat_all) - delta_gamma_all %*% p
gamma_sb_all <- lav_mat_bdiag(gamma_sb_list)
gamma_tilde <- d_mat %*% gamma_sb_all %*% t(d_mat)
# reuse lav_test_sb: structural-only Delta/E.inv/W + residualized Gamma
# (gamma_full forces the assembled 'original' path)
lav_test_sb(
lavmodel = joint@Model,
lavsamplestats = joint@SampleStats,
lavoptions = joint@Options,
lavimplied = joint@implied,
lavdata = joint@Data,
test_unscaled = std,
e_inv = e_theta_inv,
delta = delta_theta_list,
wls_v = wls_v,
m_gamma = gamma, # ignored when gamma_full is set; avoids recompute
gamma_full = gamma_tilde,
test = test
)
},
silent = TRUE
)
if (inherits(out, "try-error")) {
# the Yuan-Chan correction reuses the same step-1 jacobian (P) as the
# robust SEs, so it shares their current limitations. Equality constraints
# (within-block and across-group) and exogenous covariates (fixed.x, and
# conditional.x with the biased ADF Gamma) ARE supported. Fall back to the
# unscaled test, with a note.
lav_msg_warn(gettext(
"the Yuan & Chan (2002) scaled test statistic (test = \"yuan.chan\") is
not available for this model; the standard (unscaled) test is reported
instead."))
return(fallback) # fall back to the unscaled standard test
}
# name the scaled entry "yuan.chan" (the SAM-global analogue of the
# simultaneous "satorra.bentler"); @test = list(standard, yuan.chan)
out <- rename_yc(out)
# baseline (independence) model: no measurement model -> the Yuan-Chan
# correction reduces to the ordinary Satorra-Bentler scaling. Refit the
# baseline with the matching scaled test so the robust CFI/RMSEA can be
# computed (and so fitMeasures()/summary() do not choke on a model that has a
# scaled test while its baseline does not).
baseline_test <- NULL
if (!is.null(joint@baseline$partable)) {
refit_baseline <- function(gamma_unbiased = NULL) {
opts <- fit@Options
opts$se <- "none"
opts$test <- test
opts$baseline <- FALSE
opts$fit.by.level <- FALSE
opts$estimator <- joint@Model@estimator
if (!is.null(gamma_unbiased)) {
opts$gamma.unbiased <- gamma_unbiased
}
fit_base <- lavaan::lavaan(
model = joint@baseline$partable,
slot_data = joint@Data,
slot_sample_stats = joint@SampleStats,
slot_options = opts,
verbose = FALSE
)
rename_yc(fit_base@test)
}
baseline_test <- try(refit_baseline(), silent = TRUE)
if (inherits(baseline_test, "try-error")) {
# the unbiased ADF Gamma (the sam default) is not available for eg
# fixed.x models: retry with the biased ADF Gamma, exactly as the
# model test above does (lav_gamma_used() fallback). Without a scaled
# baseline the scaled/robust incremental fit measures would be NA.
baseline_test <- try(refit_baseline(gamma_unbiased = FALSE),
silent = TRUE)
}
if (inherits(baseline_test, "try-error")) {
baseline_test <- NULL
}
}
list(test = out, baseline.test = baseline_test)
}
# Truncation-bias report: when the (first- or second-order) lambda
# truncation engaged (lambda_star < 1), the structural summary statistics
# retain (1 - lambda_star) * MTM (resp. var.error) of extra error
# (co)variance, and the step-2 estimates are shrunken toward zero. This
# helper collects the engaged state, plus a first-order (delta-method)
# approximation of the per-coefficient shrinkage bias
# bias = d theta / d stats . (1 - lambda_star) vech(MTM)
# i.e. theta at the truncated statistics minus theta at the
# (estimand-consistent) multiplier-1 statistics, linearized at the
# truncated fit. The result is stored in the @internal slot and reported
# as a note by summary(); the estimates themselves are NOT modified.
# The magnitude is only computed in the plain setting (single level, no
# conditional.x, no categorical, no across-group Case B); otherwise only
# the engaged state is returned. Returns NULL when nothing engaged.
lav_sam_trunc_bias <- function(step1 = NULL, fit_pa = NULL) {
lam <- unlist(step1$lambda)
if (is.null(lam) || !is.numeric(lam)) {
return(NULL)
}
engaged <- is.finite(lam) & lam < 1
if (!any(engaged)) {
return(NULL)
}
vn1 <- colnames(step1$VETA[[1]])
out <- list(
lambda_star = lam, engaged = engaged,
order = if (any(grepl(":", vn1, fixed = TRUE))) 2L else 1L,
bias = NULL
)
if (isTRUE(step1$caseB) || fit_pa@Model@conditional.x ||
fit_pa@Model@categorical || fit_pa@Data@nlevels > 1L) {
return(out)
}
out$bias <- tryCatch({
delta <- lavTech(fit_pa, "delta")
wls_v <- lavTech(fit_pa, "wls.v")
if (!is.list(delta)) delta <- list(delta)
if (!is.list(wls_v)) wls_v <- list(wls_v)
ng <- length(delta)
fg <- unlist(fit_pa@SampleStats@nobs) / fit_pa@SampleStats@ntotal
a_mat <- 0
rhs <- 0
for (g in seq_len(ng)) {
# the extra error (co)variance retained by the truncation, in the
# (possibly sub-block selected) order of the step-2 statistics; note
# that the variable order of the structural fit may differ from the
# order in VETA (the sample statistics are matched by name)
vn <- fit_pa@Data@ov.names[[g]]
if (!all(vn %in% colnames(step1$VETA[[g]]))) {
vn <- colnames(step1$VETA[[g]])
}
mtm <- step1$MTM[[g]]
if (!is.null(colnames(mtm))) {
idx <- match(vn, colnames(mtm))
} else {
idx <- match(vn, step1$LV.NAMES[[g]])
}
dstats <- numeric(nrow(delta[[g]]))
if (engaged[g] && !anyNA(idx)) {
vb <- (1 - lam[g]) *
lav_mat_vech(mtm[idx, idx, drop = FALSE])
# the vech(VETA) entries are the trailing block of the step-2
# statistics (any mean/intercept entries come first)
dstats[length(dstats) - length(vb) + seq_along(vb)] <- vb
}
vd <- wls_v[[g]] %*% delta[[g]]
a_mat <- a_mat + fg[g] * crossprod(delta[[g]], vd)
rhs <- rhs + fg[g] * crossprod(vd, dstats)
}
if (fit_pa@Model@eq.constraints) {
k_mat <- fit_pa@Model@eq.constraints.K
b_full <- drop(k_mat %*% solve(
crossprod(k_mat, a_mat %*% k_mat), crossprod(k_mat, rhs)))
} else {
b_full <- drop(solve(a_mat, rhs))
}
names(b_full) <- names(coef(fit_pa))
b_full
}, error = function(e) NULL)
# the exact counterpart: by how much do the structural estimates differ
# from the ones obtained under the historical truncation rule (floor
# 1/(n-1))? The linearized bias above is blind to the nonlinearity that
# characterizes the real danger zone (an ill-conditioned predictor block,
# where the historical estimates explode); there, this shift is large,
# while it is small if the truncation engaged for a harmless reason (eg
# an almost perfectly predicted latent outcome). First-order truncation
# only; one extra (point estimates only) fit of the structural part.
out$move <- tryCatch({
if (out$order != 1L) lav_msg_stop(gettext("First-order only."))
star_used <- unlist(step1$lambda1_star_used)
floor_used <- unlist(step1$lambda1_floor_used)
nobs <- unlist(fit_pa@SampleStats@nobs)
ng <- length(step1$VETA)
veta_h <- vector("list", ng)
changed <- FALSE
for (g in seq_len(ng)) {
veta_g <- step1$VETA[[g]]
veta_h[[g]] <- veta_g
if (!engaged[g]) next
if (!is.finite(star_used[g]) || star_used[g] <= 0 ||
!is.finite(floor_used[g])) {
lav_msg_stop(gettext("Raw lambda not available.")) # clamped at zero
}
lambda_raw <- star_used[g] + floor_used[g]
floor_h <- 1 / (nobs[g] - 1)
star_h <- if (lambda_raw < 1 + floor_h) lambda_raw - floor_h else 1
if (star_h - star_used[g] < 1e-12) next # same rule was active
vn <- colnames(veta_g)
mtm <- step1$MTM[[g]]
idx <- if (!is.null(colnames(mtm))) {
match(vn, colnames(mtm))
} else {
match(vn, step1$LV.NAMES[[g]])
}
if (anyNA(idx)) lav_msg_stop(gettext("Names do not match."))
d_std <- rep(1, length(idx))
if (!is.null(step1$D.STD) && !is.null(step1$D.STD[[g]])) {
d_std <- step1$D.STD[[g]][idx]
}
veta_raw <- t(veta_g * d_std) * d_std # undo the std.lv rescaling
veta_raw <- veta_raw -
(star_h - star_used[g]) * mtm[idx, idx, drop = FALSE]
if (any(d_std != 1)) {
d_new <- d_std
d_new[d_std != 1] <- sqrt(diag(veta_raw))[d_std != 1]
veta_raw <- t(veta_raw / d_new) / d_new
}
dimnames(veta_raw) <- dimnames(veta_g)
veta_h[[g]] <- veta_raw
changed <- TRUE
}
if (!changed) {
move <- 0 * coef(fit_pa)
} else {
pt_h <- as.list(fit_pa@ParTable)
pt_h$start <- pt_h$est
pt_h$est <- pt_h$se <- NULL
opt_h <- fit_pa@Options
opt_h$se <- "none"
opt_h$test <- "none"
opt_h$baseline <- FALSE
opt_h$h1 <- FALSE
opt_h$check.post <- FALSE
opt_h$check.gradient <- FALSE
opt_h$verbose <- FALSE
opt_h$warn <- FALSE
mean_h <- NULL
if (fit_pa@Model@meanstructure) {
mean_h <- fit_pa@SampleStats@mean
for (g in seq_along(mean_h)) {
names(mean_h[[g]]) <- fit_pa@Data@ov.names[[g]]
}
}
fit_h <- suppressWarnings(lavaan::lavaan(pt_h,
sample_cov = veta_h, sample_mean = mean_h,
sample_nobs = fit_pa@SampleStats@nobs,
slot_options = opt_h, verbose = FALSE
))
if (!lavInspect(fit_h, "converged")) {
lav_msg_stop(gettext("No convergence."))
}
move <- coef(fit_pa) - coef(fit_h)
}
move
}, error = function(e) NULL)
out
}
# Truncation report, part 2 (used by summary()): put the approximate
# shrinkage bias, and the exact shift away from the historical-rule
# estimates, on the scale of the standard errors (which are only available
# in the final object). Adds to the stored sam_trunc list:
# bias_max = max |bias| among the regression coefficients (or among all
# structural parameters if there are no regressions)
# bias_se_max = max |bias| / se among the same parameters
# move_max, move_se_max = idem, for the shift ('move')
# (NA if not available). The note printed by summary() is proportional: it
# is merely informative if BOTH bias_se_max and move_se_max are available
# and at most 0.5 (a bias of half a standard error lowers the coverage of a
# 95% confidence interval to about 92%); in all other cases -- including
# the settings where the shift cannot be computed -- the full warning is
# printed.
lav_sam_trunc_bias_se <- function(object = NULL) {
sam_trunc <- object@internal$sam_trunc
if (is.null(sam_trunc) || is.null(sam_trunc$bias)) {
return(sam_trunc)
}
pt <- object@ParTable
# same labels as coef(): user labels if any, else lhs-op-rhs (+ group
# suffix); equality-constrained duplicates share label, est and se
pt_labels <- tryCatch(lav_pt_labels(pt, type = "user"),
error = function(e) NULL)
if (length(pt_labels) != length(pt$lhs)) {
pt_labels <- paste0(pt$lhs, pt$op, pt$rhs)
}
# the regression coefficients (the variance parameters absorb most of the
# retained error variance by construction)
reg_names <- unique(pt_labels[pt$op == "~" & pt$free > 0L])
se <- NULL
if (!is.null(pt$se) && !identical(object@Options$se, "none")) {
se <- pt$se[pt$free > 0L]
names(se) <- pt_labels[pt$free > 0L]
}
max_abs <- function(v) {
out <- c(as.numeric(NA), as.numeric(NA))
if (is.null(v)) {
return(out)
}
if (any(names(v) %in% reg_names)) {
v <- v[names(v) %in% reg_names]
}
v <- v[is.finite(v)]
if (length(v) == 0L) {
return(out)
}
out[1] <- max(abs(v))
if (!is.null(se)) {
ratio <- abs(v) / se[match(names(v), names(se))]
ratio <- ratio[is.finite(ratio)]
if (length(ratio) > 0L) {
out[2] <- max(ratio)
}
}
out
}
tmp <- max_abs(sam_trunc$bias)
sam_trunc$bias_max <- tmp[1]
sam_trunc$bias_se_max <- tmp[2]
tmp <- max_abs(sam_trunc$move)
sam_trunc$move_max <- tmp[1]
sam_trunc$move_se_max <- tmp[2]
sam_trunc
}
lav_sam_table <- function(joint = NULL, step1 = NULL, fit_pa = NULL,
cmd = NULL, lavoptions = NULL,
mm_args = list(), struc_args = list(),
sam_method = "local",
local_options = list(), global_options = list()) {
mm_fit <- step1$MM.FIT
sam_mm_table <- data.frame(
Block = seq_along(step1$mm.list),
Latent = sapply(mm_fit, function(x) {
paste(unique(unlist(x@pta$vnames$lv)), collapse = ",")
}),
Nind = sapply(mm_fit, function(x) {
length(unique(unlist(x@pta$vnames$ov)))
}),
Chisq = sapply(mm_fit, function(x) {
x@test[[1]]$stat
}),
Df = sapply(mm_fit, function(x) {
x@test[[1]]$df
})
)
class(sam_mm_table) <- c("lavaan.data.frame", "data.frame")
# extra info for @internal slot
sam_trunc <- NULL
if (sam_method %in% c("local", "fsr", "cfsr")) {
# the structural fit object, with the (two-step robust) corrected test
# and baseline injected; for multigroup models with across-group
# measurement constraints (Case B) these are recomputed with the full
# cross-group Gamma. The object is stored in the @internal slot, so that
# fitMeasures(), lavResiduals() and modindices() can delegate to it
# (issue #517)
sam_struc_fit_object <- lav_sam_struc_fit_object(
fit_pa = fit_pa, step1 = step1
)
# (two-step robust) fit measures for the structural part
sam_struc_fit <- lav_sam_struc_fit(fit_pa = sam_struc_fit_object)
sam_mm_rel <- step1$REL
# engaged lambda truncation? collect the state + the approximate
# shrinkage bias of the structural estimates (summary() note)
sam_trunc <- lav_sam_trunc_bias(step1 = step1, fit_pa = fit_pa)
} else {
sam_struc_fit <- paste0("no local fit measures available for",
"structural part if sam.method is global")
names(sam_struc_fit) <- "warning"
sam_struc_fit_object <- NULL
sam_mm_rel <- numeric(0L)
}
# record the EFFECTIVE structural information.bread: lav_sam_step2() may
# have resolved "default" to "observed" (the default bread for the local
# standard errors since 0.7-2), and summary() reads it from
# sam.struc.args (see lavParameterEstimates)
if (is.null(struc_args[["information.bread"]]) &&
!is.null(fit_pa@Options$information.bread) &&
fit_pa@Options$information.bread != "default") {
struc_args$information.bread <- fit_pa@Options$information.bread
}
sam_1 <- list(
sam.cmd = cmd,
sam.method = sam_method,
sam.local.options = local_options,
sam.global.options = global_options,
sam.mm.list = step1$mm.list,
sam.mm.estimator = mm_fit[[1]]@Model@estimator,
sam.mm.args = mm_args,
sam.mm.ov.names = lapply(mm_fit, function(x) {
x@pta$vnames$ov
}),
sam.mm.table = sam_mm_table,
sam.mm.rel = sam_mm_rel,
sam_trunc = sam_trunc,
sam.struc.estimator = fit_pa@Model@estimator,
sam.struc.args = struc_args,
sam.struc.fit = sam_struc_fit,
sam.struc.fit.object = sam_struc_fit_object,
sam.lavoptions = lavoptions
)
sam_1
}
lav_sam_get_cov_ybar <- function(fit = NULL, local_options = list(
m_method = "ML",
lambda_correction = TRUE,
lambda1_floor = "debias",
lambda2_floor = "debias",
alpha_correction = 0L,
twolevel_method = "h1"
)) {
# local.twolevel.method
local_twolevel_method <- tolower(local_options[["twolevel_method"]])
if (!local_twolevel_method %in% c("h1", "anova", "mean")) {
lav_msg_stop(gettext(
"local option twolevel_method should be one of h1, anova or mean."))
}
lavpta <- fit@pta
nlevels <- lavpta$nlevels
nblocks <- lavpta$nblocks
# do we need H1?
if (nlevels > 1L && local_twolevel_method == "h1") {
h1 <- lav_h1_implied_logl(
lavdata = fit@Data,
lavsamplestats = fit@SampleStats,
lavoptions = fit@Options
)
h1implied <- h1$implied
} else {
h1implied <- fit@h1$implied
}
# containers
cov_list <- vector("list", nblocks)
ybar_list <- vector("list", nblocks)
slopes_list <- vector("list", nblocks)
# label
if (nblocks > 1L) {
names(cov_list) <- fit@Data@block.label
names(ybar_list) <- fit@Data@block.label
names(slopes_list) <- fit@Data@block.label
}
# collect COV/YBAR per block
for (b in seq_len(nblocks)) {
# get sample statistics for this block
if (nlevels > 1L) {
# block order is group-major: (g1l1, g1l2, g2l1, g2l2, ...)
this_level <- (b - 1L) %% nlevels + 1L
this_group <- floor((b - 1L) / nlevels) + 1L
if (this_level > 2L) {
lav_msg_stop(gettext("level 3 not supported (yet)."))
}
if (local_twolevel_method == "h1") {
# note: the h1 estimates are already reduced to the
# level-specific variables -- no need to subset
cov_1 <- h1implied$cov[[(this_group - 1L) * nlevels + this_level]]
ybar <- h1implied$mean[[(this_group - 1L) * nlevels + this_level]]
} else {
ylp <- fit@SampleStats@YLp[[this_group]][[2]]
if (this_level == 1L) {
# "anova" and "mean" use the same within statistics
cov_1 <- ylp$Sigma.W
ybar <- ylp$Mu.W
} else if (local_twolevel_method == "anova") {
cov_1 <- ylp$Sigma.B
ybar <- ylp$Mu.B
} else { # "mean"
s_pw <- ylp$Sigma.W
nj <- ylp$s
y2 <- ylp$Y2
# grand mean
mu_y <- ylp$Mu.W + ylp$Mu.B
y2c <- t(t(y2) - mu_y) # MUST be centered
yb <- crossprod(y2c) / nrow(y2c)
cov_1 <- yb - 1 / nj * s_pw
ybar <- ylp$Mu.B
}
# reduce to the level-specific variables
ov_idx <- fit@Data@Lp[[this_group]]$ov.idx[[this_level]]
cov_1 <- cov_1[ov_idx, ov_idx, drop = FALSE]
ybar <- ybar[ov_idx]
}
# single level
} else {
this_group <- b
if (fit@Model@conditional.x) {
ybar <- h1implied$res.int[[b]]
cov_1 <- h1implied$res.cov[[b]]
slopes_list[[b]] <- h1implied$res.slopes[[b]]
} else {
ybar <- h1implied$mean[[b]] # EM version if missing="ml"
cov_1 <- h1implied$cov[[b]]
}
} # single level
cov_list[[b]] <- cov_1
ybar_list[[b]] <- ybar
}
list(COV = cov_list, YBAR = ybar_list, SLOPES = slopes_list)
}
# automatically generate mm.list; we create measurement blocks so that:
# - unlinked factors become a singleton
# - linked factor are joined into a single measurement block
# based: cross-loadings, or correlated residuals between item errors
lav_sam_get_mmlist <- function(lavobject) {
lavmodel <- lavobject@Model
lavpta <- lavobject@pta
nblocks <- lavpta$nblocks
lambda_idx <- which(names(lavmodel@GLIST) == "lambda")
mm_list <- vector("list", length = nblocks)
glist <- lavTech(lavobject, "partable")
for (b in seq_len(nblocks)) {
mm_in_block <- (seq_len(lavmodel@nmat[b]) +
cumsum(c(0, lavmodel@nmat))[b])
mlist <- glist[mm_in_block]
cc <- t(mlist$lambda) %*% mlist$theta %*% mlist$lambda
# note: CC contains dummy lv's, higher-order, etc... and they
# must be removed
# get ALL lv names (including dummy ov.x/ov.y)
ov_names <- lavmodel@dimNames[[lambda_idx[b]]][[1L]]
lv_names <- lavmodel@dimNames[[lambda_idx[b]]][[2L]]
names_1 <- lv_names
# needs to removed:
# - all lv.ind variables
# - all higher-order latent variables
# - all dummy lv's
rm_idx <- c(match(lavpta$vnames$lv.ind[[b]], lv_names),
match(lavpta$vnames$lv.interaction[[b]], lv_names),
which(lv_names %in% ov_names))
if (length(rm_idx) > 0L) {
cc <- cc[-rm_idx, -rm_idx, drop = FALSE]
names_1 <- lv_names[-rm_idx]
}
# cluster membership
membership <- lav_graph_get_connected_nodes(cc)
out <- split(names_1, names_1[membership])
names(out) <- paste("block", seq_along(out), sep = "")
mm_list[[b]] <- out
}
mm_list
}
# labels for the elements of the WLS.obs statistics vector in the
# categorical case: 1) thresholds and (negative) means, interleaved per
# variable, 2) [conditional.x only] the slopes of the y-variables on the
# exogenous covariates, column-major (lav_mat_vec(SLOPES): all y for x1,
# then x2, ...), 3) variances (numeric variables only), 4) correlations
# (vech, no diagonal); used to map the statistics of a measurement block
# into the statistics vector of the joint model
lav_sam_wls_obs_labels <- function(ov_names, th_idx, ov_names_x = NULL) {
nvar <- length(ov_names)
# variable of each th element: ordered variables have th_idx > 0; the
# th_idx == 0 elements are the (negative) means of the numeric
# variables, appearing in variable order
v_of <- th_idx
num_var_idx <- which(!seq_len(nvar) %in% th_idx)
v_of[th_idx == 0L] <- num_var_idx
th_labels <- character(length(th_idx))
count <- integer(nvar)
for (k in seq_along(th_idx)) {
v <- v_of[k]
count[v] <- count[v] + 1L
th_labels[k] <- paste0(ov_names[v], "|t", count[v])
}
# slopes (conditional.x only), column-major over the nvar x nexo matrix
slope_labels <- character(0L)
if (length(ov_names_x) > 0L) {
slope_labels <- paste0(rep(ov_names, times = length(ov_names_x)), "~",
rep(ov_names_x, each = nvar))
}
# variances (numeric variables only)
var_labels <- character(0L)
if (length(num_var_idx) > 0L) {
var_labels <- paste0(ov_names[num_var_idx], "|var")
}
# correlations: lav_mat_vech(cov, diagonal = FALSE), column-major;
# use a canonical (sorted) pair label, as the relative order of two
# variables may differ between the joint model and a measurement block
cor_labels <- character(0L)
if (nvar > 1L) {
row_idx <- matrix(seq_len(nvar), nvar, nvar)[lower.tri(diag(nvar))]
col_idx <- matrix(seq_len(nvar), nvar, nvar,
byrow = TRUE)[lower.tri(diag(nvar))]
cor_labels <- vapply(seq_along(row_idx), function(k) {
paste(sort(c(ov_names[row_idx[k]], ov_names[col_idx[k]])),
collapse = "~~")
}, character(1L))
}
c(th_labels, slope_labels, var_labels, cor_labels)
}
# conditional.x: permutation of the conditional structural-moment vector
# c( vecr(cbind(EETA, RS_eta)), vech(VETA) ) -- with the variables in
# VETA (measurement) order and the covariates in data (ov.names.x) order --
# to the structural model's internal order (target_eta = its ov.nox,
# target_x = its ov.x). This is the layout of FIT.PA's WLS.obs: the NACOV
# passed to lavaan() is consumed as-is (no name matching), so Gamma.eta must
# be permuted with this index before the step-2 handoff.
lav_sam_condx_perm_idx <- function(eta_names, x_names,
target_eta, target_x) {
ord <- match(target_eta, eta_names)
ord_x <- match(target_x, x_names)
if (length(ord) != length(eta_names) || anyNA(ord) ||
length(ord_x) != length(x_names) || anyNA(ord_x)) {
lav_msg_stop(gettext(
"internal error: unable to align the conditional structural moments
with the structural model (Gamma.eta permutation)"))
}
p <- length(eta_names)
q <- length(x_names)
mean_block <- unlist(lapply(ord, function(j) {
(j - 1L) * (q + 1L) + c(1L, 1L + ord_x)
}))
pos <- matrix(0L, p, p)
pos[lower.tri(pos, diag = TRUE)] <- seq_len(p * (p + 1L) / 2L)
pos[upper.tri(pos)] <- t(pos)[upper.tri(pos)]
sub <- pos[ord, ord, drop = FALSE]
vech_perm <- sub[lower.tri(sub, diag = TRUE)]
c(mean_block, p * (q + 1L) + vech_perm)
}
# merge the per-block mm clusters (as returned by lav_sam_get_mmlist())
# into a single mm.list, transitively joining clusters that share a latent
# variable or an observed indicator; needed for multilevel models, where
# the within and between factors of the same indicators must end up in the
# same (two-level) measurement block
lav_sam_mmlist_merge_blocks <- function(lavobject, mm_list_per_block) {
pt_1 <- lavobject@ParTable
nblocks <- lavobject@pta$nblocks
block_values <- lav_pt_block_values(pt_1)
# collect all clusters: a set of lv names + their observed indicators
clusters <- list()
for (b in seq_len(nblocks)) {
for (cl in mm_list_per_block[[b]]) {
lv <- unlist(cl)
ind_idx <- which(pt_1$op == "=~" & pt_1$block == block_values[b] &
pt_1$lhs %in% lv)
clusters <- c(clusters,
list(list(lv = lv, ov = unique(pt_1$rhs[ind_idx]))))
}
}
# merge until no two clusters share an lv or an ov
changed <- TRUE
while (changed) {
changed <- FALSE
i <- 1L
while (i < length(clusters)) {
j <- i + 1L
while (j <= length(clusters)) {
if (any(clusters[[i]]$lv %in% clusters[[j]]$lv) ||
any(clusters[[i]]$ov %in% clusters[[j]]$ov)) {
clusters[[i]]$lv <- unique(c(clusters[[i]]$lv, clusters[[j]]$lv))
clusters[[i]]$ov <- unique(c(clusters[[i]]$ov, clusters[[j]]$ov))
clusters[[j]] <- NULL
changed <- TRUE
} else {
j <- j + 1L
}
}
i <- i + 1L
}
}
lapply(clusters, "[[", "lv")
}
lav_sam_veta_pt <- function(lavobject, block = 1L) {
lavmodel <- lavobject@Model
lambda_idx <- which(names(lavmodel@GLIST) == "lambda")
glist <- lavTech(lavobject, "partable")
b <- block
mm_in_block <- (seq_len(lavmodel@nmat[b]) + cumsum(c(0, lavmodel@nmat))[b])
mlist <- glist[mm_in_block]
mm_psi <- (mlist$psi != 0) + 0L
nr <- nrow(mlist$psi)
if (!is.null(mlist$beta)) {
mm_beta <- (mlist$beta != 0) + 0L
tmp <- -mm_beta
tmp[lav_mat_diag_idx(nr)] <- 1
ib_inv <- solve(tmp)
} else {
ib_inv <- diag(nr)
}
veta <- ib_inv %*% mm_psi %*% t(ib_inv)
# get ALL lv names (including dummy ov.x/ov.y)
lv_names <- lavmodel@dimNames[[lambda_idx[b]]][[2L]]
colnames(veta) <- rownames(veta) <- lv_names
veta
}
lav_sam_veta_con <- function(s = NULL, mm_lambda = NULL, mm_theta = NULL,
l_veta = NULL, local_m_method = "ML",
tol = 1e-07, max_iter = 100L) {
# first do GLS/ULS
m <- ncol(mm_lambda)
if (local_m_method == "ULS") {
m_w <- diag(m)
} else {
m_w <- solve(s)
}
wl <- m_w %*% mm_lambda
t_lw <- t(wl)
t_lwl <- t_lw %*% mm_lambda
tmp1 <- t(l_veta) %*% (t_lwl %x% t_lwl) %*% l_veta
tmp2 <- t(l_veta) %*% lav_mat_vec(t_lw %*% (s - mm_theta) %*% wl)
out <- solve(tmp1, tmp2)
veta_init <- matrix(l_veta %*% out, m, m)
if (local_m_method != "ML") {
# we are done!
return(veta_init)
}
veta_new <- veta_init
for (i in seq_len(max_iter)) {
sigma_ml <- mm_lambda %*% veta_new %*% t(mm_lambda) + mm_theta
m_w <- solve(sigma_ml)
wl <- m_w %*% mm_lambda
t_lw <- t(wl)
t_lwl <- t_lw %*% mm_lambda
tmp1 <- t(l_veta) %*% (t_lwl %x% t_lwl) %*% l_veta
tmp2 <- t(l_veta) %*% lav_mat_vec(t_lw %*% (s - mm_theta) %*% wl)
out <- solve(tmp1, tmp2)
veta_ml <- matrix(l_veta %*% out, m, m)
rmsea <- sqrt(sum((veta_new - veta_ml)^2))
#cat("i = ", i, " rmsea = ", rmsea, "\n")
veta_new <- veta_ml
if (rmsea < tol) {
break
}
}
veta_new
}
Try the lavaan package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.