add "lv.ho" in lavNames() to list higher-order lv's; prepare sam() fo…

…r higher-order latent factors

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   = "[email protected]",

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 == "~"])

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

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
           }