lavTestLRT(): add RMSEA if a robust test statistic is used (but only if data is complete)

Author: yrosseel

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: lavaan
 Title: Latent Variable Analysis
-Version: 0.6-20.2325
+Version: 0.6-20.2326
 Authors@R: c(person(given   = "Yves", family = "Rosseel",
                     role    = c("aut", "cre"),
                     email   = "Yves.Rosseel@UGent.be",

R/lav_options.R

@@ -367,7 +367,8 @@ lav_options_set <- function(opt = NULL) {                     # nolint
     ))) {
       lav_msg_stop(gettextf(
         "missing=%1$s is not allowed for estimator %2$s",
-        dQuote(opt$missing), dQuote(lav_options_estimatorgroup(opt$estimator))))
+        dQuote(opt$missing),
+        dQuote(lav_options_estimatorgroup(opt$estimator.orig))))
     }
   } else if (opt$missing == "ml.x") {
     if (opt$.categorical) {
@@ -381,7 +382,8 @@ lav_options_set <- function(opt = NULL) {                     # nolint
     ))) {
       lav_msg_stop(gettextf(
         "missing=%1$s is not allowed for estimator %2$s",
-        dQuote(opt$missing), dQuote(lav_options_estimatorgroup(opt$estimator))))
+        dQuote(opt$missing),
+        dQuote(lav_options_estimatorgroup(opt$estimator.orig))))
     }
   } else if (opt$missing == "two.stage") {
     if (opt$.categorical) {
@@ -395,7 +397,8 @@ lav_options_set <- function(opt = NULL) {                     # nolint
     ))) {
       lav_msg_stop(gettextf(
         "missing=%1$s is not allowed for estimator %2$s",
-        dQuote(opt$missing), dQuote(lav_options_estimatorgroup(opt$estimator))))
+        dQuote(opt$missing),
+        dQuote(lav_options_estimatorgroup(opt$estimator.orig))))
     }
   } else if (opt$missing == "robust.two.stage") {
     if (opt$.categorical) {
@@ -409,7 +412,8 @@ lav_options_set <- function(opt = NULL) {                     # nolint
     ))) {
       lav_msg_stop(gettextf(
         "missing=%1$s is not allowed for estimator %2$s",
-        dQuote(opt$missing), dQuote(lav_options_estimatorgroup(opt$estimator))))
+        dQuote(opt$missing),
+        dQuote(lav_options_estimatorgroup(opt$estimator.orig))))
     }
   } else if (opt$missing == "doubly.robust") {
     if (opt$estimator != "pml") {

R/lav_sam_step1_local.R

@@ -371,7 +371,7 @@ lav_sam_step1_local <- function(STEP1 = NULL, FIT = NULL, Y = NULL,
         # EETA is constrained somehow
         stop("not ready yet")
         # EETA[[b]] <- lav_sam_eeta_con(YBAR = YBAR, LAMBDA = LAMBDA[[b]],
-        #                               THETA = THETA[[b]], 
+        #                               THETA = THETA[[b]],
         #                               L.veta = L.veta[[b]])
 
       }

R/lav_sam_utils.R

@@ -951,8 +951,8 @@ lav_sam_veta_partable <- function(lavobject, block = 1L) {
   VETA
 }
 
-lav_sam_veta_con <- function(S = NULL, LAMBDA = NULL, THETA = NULL, 
-                             L.veta = NULL, local.M.method = "ML", 
+lav_sam_veta_con <- function(S = NULL, LAMBDA = NULL, THETA = NULL,
+                             L.veta = NULL, local.M.method = "ML",
                              tol = 1e-07, max.iter = 100L) {
 
   # first do GLS/ULS

R/lav_test_LRT.R

@@ -11,7 +11,8 @@
 #     method = "Satorra.Bentler.2010"
 #     method = "mean.var.adjusted.PLRT"
 #
-# - 0.6-13: RMSEA.D (also known as 'RDR') is added to the table (unless scaled)
+# - 0.6-13: RMSEA.D (also known as 'RDR') is added to the table (also if scaled,
+#           since 0.6-20)
 # - 0.6-13: fix multiple-group UG^2 bug in Satorra.2000 (reported by
 #           Gronneberg, Foldnes and Moss)
 #
@@ -291,14 +292,6 @@ lavTestLRT <- function(object, ..., method = "default", test = "default",
                            "satorra.2000")) {
     c.delta <- rep(as.numeric(NA), length(STAT))
   }
-  # new in 0.6-13
-  if (!scaled) {
-    RMSEA.delta <- c(NA, lav_fit_rmsea(
-      X2 = STAT.delta[-1],
-      df = Df.delta[-1],
-      N = ntotal, G = ngroups
-    ))
-  }
 
   # check for negative values in STAT.delta
   # but with a tolerance (0.6-12)!
@@ -366,9 +359,28 @@ lavTestLRT <- function(object, ..., method = "default", test = "default",
     }
   }
 
+  # unname
+  STAT.delta <- unname(STAT.delta)
+  Df.delta <- unname(Df.delta)
+  if (scaled) {
+    c.delta <- unname(c.delta)
+  }
+
   # Pvalue
   Pvalue.delta <- pchisq(STAT.delta, Df.delta, lower.tail = FALSE)
 
+  # new in 0.6-13: RMSEA (RMSEA.D or RDR)
+  if (object@Options$missing == "listwise") {
+    RMSEA.delta <- c(NA, lav_fit_rmsea(
+      X2 = STAT.delta[-1],
+      df = Df.delta[-1],
+      N = ntotal,
+      G = ngroups,
+      c.hat = ifelse (scaled, c.delta[-1], rep(1, length(STAT.delta) - 1L))
+    ))
+  }
+
+  # AIC/BIC
   aic <- bic <- rep(NA, length(mods))
   if (estimator == "ML") {
     aic <- sapply(mods, FUN = AIC)
@@ -391,33 +403,32 @@ lavTestLRT <- function(object, ..., method = "default", test = "default",
       row.names = names(mods),
       check.names = FALSE
     )
+  } else if (object@Options$missing == "listwise") {
+    val <- data.frame(
+      Df = Df,
+      AIC = aic,
+      BIC = bic,
+      Chisq = STAT,
+      "Chisq diff" = STAT.delta,
+      "RMSEA" = RMSEA.delta,
+      "Df diff" = Df.delta,
+      "Pr(>Chisq)" = Pvalue.delta,
+      row.names = names(mods),
+      check.names = FALSE
+    )
   } else {
-    if (scaled) {
-      val <- data.frame(
-        Df = Df,
-        AIC = aic,
-        BIC = bic,
-        Chisq = STAT,
-        "Chisq diff" = STAT.delta,
-        "Df diff" = Df.delta,
-        "Pr(>Chisq)" = Pvalue.delta,
-        row.names = names(mods),
-        check.names = FALSE
-      )
-    } else {
-      val <- data.frame(
-        Df = Df,
-        AIC = aic,
-        BIC = bic,
-        Chisq = STAT,
-        "Chisq diff" = STAT.delta,
-        "RMSEA" = RMSEA.delta,
-        "Df diff" = Df.delta,
-        "Pr(>Chisq)" = Pvalue.delta,
-        row.names = names(mods),
-        check.names = FALSE
-      )
-    }
+    val <- data.frame(
+      Df = Df,
+      AIC = aic,
+      BIC = bic,
+      Chisq = STAT,
+      "Chisq diff" = STAT.delta,
+      #"RMSEA" = RMSEA.delta, # if missing, not yet...
+      "Df diff" = Df.delta,
+      "Pr(>Chisq)" = Pvalue.delta,
+      row.names = names(mods),
+      check.names = FALSE
+    )
   }
 
   # catch Df.delta == 0 cases (reported by Florian Zsok in Zurich)

R/xxx_sam.R

@@ -189,7 +189,7 @@ sam <- function(model = NULL,
   # check for gamma.unbiased
   if (is.null(dotdotdot$gamma.unbiased)) {
      # put in TRUE in dotdotdot # lavaan default is still FALSE
-     dotdotdot$gamma.unbiased <- TRUE  
+     dotdotdot$gamma.unbiased <- TRUE
   }
 
 
@@ -312,7 +312,8 @@ sam <- function(model = NULL,
     STEP1 <- lav_sam_step1_local(
       STEP1 = STEP1, FIT = FIT,
       sam.method = sam.method,
-      local.options = local.options
+      local.options = local.options,
+      return.cov.iveta2 = (se %in% c("local", "local.nt"))
     )
   }