From b163308b507d8ede93c7438f5e67c5c8d4655b67 Mon Sep 17 00:00:00 2001 From: Yves Rosseel Date: Fri, 14 Feb 2025 17:00:42 +0100 Subject: [PATCH] add "lv.ho" in lavNames() to list higher-order lv's; prepare sam() for higher-order latent factors --- DESCRIPTION | 2 +- R/lav_partable_vnames.R | 27 +++++++++++++++++ R/lav_sam_utils.R | 65 +++++++++++++++++++++++++++++++++++++---- R/lav_start.R | 3 +- 4 files changed, 90 insertions(+), 7 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index 617da0e8..62882541 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: lavaan Title: Latent Variable Analysis -Version: 0.6-20.2264 +Version: 0.6-20.2265 Authors@R: c(person(given = "Yves", family = "Rosseel", role = c("aut", "cre"), email = "Yves.Rosseel@UGent.be", diff --git a/R/lav_partable_vnames.R b/R/lav_partable_vnames.R index 4c771ece..30e29297 100644 --- a/R/lav_partable_vnames.R +++ b/R/lav_partable_vnames.R @@ -36,6 +36,7 @@ lavaanNames <- lavNames # nolint # or returns a list per block # LDW 30/1/24: 'block' argument not explicitly tested !? # LDW 29/2/24: ov.order = "data" via attribute "ovda" +# - YR 11/02/25: add type = "lv.ho" for higher-order latent variables lav_partable_vnames <- function(partable, type = NULL, ..., # nolint force.warn = FALSE, ov.x.fatal = FALSE) { # This function derives the names of some types of variable (as specified @@ -98,6 +99,7 @@ lav_partable_vnames <- function(partable, type = NULL, ..., # nolint "lv.efa", # latent variables involved in efa "lv.rv", # random slopes, random variables "lv.ind", # latent indicators (higher-order cfa) + "lv.ho", # higher-order latent variables "lv.marker", # marker indicator per lv "eqs.y", # y's in regression @@ -311,6 +313,19 @@ lav_partable_vnames <- function(partable, type = NULL, ..., # nolint next } + # higher-order latent variables + if ("lv.ho" == type) { + out.ind <- unique(partable$rhs[block.ind & + partable$op == "=~" & + partable$rhs %in% lv.names]) + out <- unique(partable$lhs[block.ind & + partable$op == "=~" & + partable$lhs %in% lv.names & + partable$rhs %in% out.ind]) + return.value$lv.ho[[b]] <- out + next + } + # eqs.y if (!(any(type == c("lv", "lv.regular")))) { eqs.y <- unique(partable$lhs[block.ind & @@ -834,6 +849,18 @@ lav_partable_vnames <- function(partable, type = NULL, ..., # nolint return.value$lv.ind[[b]] <- out } + # higher-order latent variables + if (any("lv.ho" == type)) { + out.ind <- unique(partable$rhs[block.ind & + partable$op == "=~" & + partable$rhs %in% lv.names]) + out <- unique(partable$lhs[block.ind & + partable$op == "=~" & + partable$lhs %in% lv.names & + partable$rhs %in% out.ind]) + return.value$lv.ho[[b]] <- out + } + # eqs.y eqs.y <- unique(partable$lhs[block.ind & partable$op == "~"]) diff --git a/R/lav_sam_utils.R b/R/lav_sam_utils.R index 3477aa8e..c97adcef 100644 --- a/R/lav_sam_utils.R +++ b/R/lav_sam_utils.R @@ -11,6 +11,13 @@ lav_sam_mapping_matrix <- function(LAMBDA = NULL, THETA = NULL, method <- toupper(method) + # catch empty columns in LAMBDA (eg higher order!) + LAMBDA.orig <- LAMBDA + empty.idx <- which(apply(LAMBDA, 2L, function(x) all(x == 0))) + if (length(empty.idx) > 0L) { + LAMBDA <- LAMBDA.orig[, -empty.idx, drop = FALSE] + } + # ULS # M == solve( t(LAMBDA) %*% LAMBDA ) %*% t(LAMBDA) # == MASS:::ginv(LAMBDA) @@ -119,6 +126,13 @@ lav_sam_mapping_matrix <- function(LAMBDA = NULL, THETA = NULL, } } # 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 } @@ -138,7 +152,16 @@ lav_sam_mapping_matrix_tmat <- function(LAMBDA = NULL, THETA = NULL, marker.idx = NULL, std.lv = NULL) { + LAMBDA <- as.matrix.default(LAMBDA) + + # catch empty columns in LAMBDA (eg higher order!) + LAMBDA.orig <- LAMBDA + empty.idx <- which(apply(LAMBDA, 2L, function(x) all(x == 0))) + if (length(empty.idx) > 0L) { + LAMBDA <- LAMBDA.orig[, -empty.idx, drop = FALSE] + } + nvar <- nrow(LAMBDA) nfac <- ncol(LAMBDA) @@ -179,6 +202,13 @@ lav_sam_mapping_matrix_tmat <- function(LAMBDA = NULL, 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 } @@ -245,21 +275,33 @@ lav_sam_tmat <- function(LAMBDA = NULL, lav_sam_veta <- function(M = NULL, S = NULL, THETA = NULL, alpha.correction = 0L, lambda.correction = TRUE, 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 %*% THETA %*% t(M) + # empty theta elements? + empty.theta.idx <- which(diag(MTM) == 0) + # new in 0.6-16: make sure MTM is pd # (otherwise lav_matrix_symmetric_diff_smallest_root will fail) - if (length(dummy.lv.idx) > 0L) { - MTM.nodummy <- MTM[-dummy.lv.idx, -dummy.lv.idx, drop = FALSE] - MTM.nodummy <- zapsmall(lav_matrix_symmetric_force_pd( - MTM.nodummy, + theta.rm.idx <- unique(c(dummy.lv.idx, empty.theta.idx)) + if (length(theta.rm.idx) > 0L) { + MTM.small <- MTM[-theta.rm.idx, -theta.rm.idx, drop = FALSE] + MTM.small <- zapsmall(lav_matrix_symmetric_force_pd( + MTM.small, tol = 1e-04 )) - MTM[-dummy.lv.idx, -dummy.lv.idx] <- MTM.nodummy + MTM[-theta.rm.idx, -theta.rm.idx] <- MTM.small } else { MTM <- zapsmall(lav_matrix_symmetric_force_pd(MTM, tol = 1e-04)) } @@ -299,6 +341,19 @@ lav_sam_veta <- function(M = NULL, S = NULL, THETA = NULL, 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 diff --git a/R/lav_start.R b/R/lav_start.R index e9767d24..3700b05b 100644 --- a/R/lav_start.R +++ b/R/lav_start.R @@ -304,7 +304,8 @@ lav_start <- function(start.method = "default", } # check for negative marker - if (!std.lv && lavpartable$ustart[lambda.idx[1]] < 0) { + if (!std.lv && !is.na(lavpartable$ustart[lambda.idx[1]]) && + lavpartable$ustart[lambda.idx[1]] < 0) { tmp <- -1 * tmp }