fit per group needs more work

jasp-stats · Feb 23, 2024 · 6b9c520 · 6b9c520
1 parent b9ab5c3
commit 6b9c520
Show file tree

Hide file tree

Showing 2 changed files with 28 additions and 69 deletions.
diff --git a/R/sem.R b/R/sem.R
@@ -18,6 +18,8 @@
 SEMInternal <- function(jaspResults, dataset, options, ...) {
   jaspResults$addCitation("Rosseel, Y. (2012). lavaan: An R Package for Structural Equation Modeling. Journal of Statistical Software, 48(2), 1-36. URL http://www.jstatsoft.org/v48/i02/")
 
+  sink(file = "~/Downloads/log.txt")
+  on.exit(sink(NULL))
   # Read dataset
   options <- .semPrepOpts(options)
 
@@ -585,76 +587,49 @@ checkLavaanModel <- function(model, availableVars) {
 
     groupNames <- semResults[[1]]@Data@group.label
     models <- rep(rownames(lrt$value), each = length(groupNames))
-    lavopts <- .semOptionsToLavOptions(options, dataset)
-    lavopts[["group"]] <- NULL
-    results_grouped <- list()
-    k <- 1
-    for (i in 1:length(options[["models"]])) {
-      syntax   <- lavaan::parTable(semResults[[i]])
-      for (group in seq_along(groupNames)) {
-        lav_args <- lavopts
-        syntax_group <- syntax[syntax$group == group,]
-        equalityConstraints <- syntax[syntax$op == "==", ]
-        if (nrow(equalityConstraints) > 0) {
-          for (j in 1:nrow(equalityConstraints)) {
-            if (equalityConstraints[j, "lhs"] %in% syntax_group[, "plabel"]) {
-              syntax_group <- rbind(syntax_group, equalityConstraints[j,])
-            }
-          }
-        }
-        syntax_group[["group"]] <- syntax_group[["block"]] <- rep(1, nrow(syntax_group))
-        lav_args[["model"]] <- syntax_group
-        dataset_group <- dataset[dataset[[options[["group"]]]] == groupNames[group],]
-        lav_args[["data"]]  <- dataset_group
+    modelDfs <- unlist(lapply(semResults, function(x) {lavaan::lavInspect(x, what = "test")[[testName]]$df}))
+    print(modelDfs)
 
-        fit_group <- try(do.call(lavaan::lavaan, lav_args))
+    # so the models from the lrt are ordered so the model with fewer dfs is up, we should check that and
+    # adjust accorindlgy
+    # also need to check what happens with the likelihood and a different test?
+    # shouldnt the likleihood come before the test
+    # so does the LRT for models then depend on the test in each model?
 
-        if (isTryError(fit_group)) {
-          errormsg <- gettextf("The model fit by group is unavailable for the specified model: %s", options[["models"]][[i]][["name"]])
-          fittab$setError(errormsg)
-          break
-        }
-        results_grouped[[k]] <- fit_group
-        k = k + 1
-      }
-    }
+    # also see if unlist is the way to go, probably not:
 
-    chiSq <- unlist(lapply(semResults, function(x) {lavaan::lavInspect(x, what = "test")[[testName]]$stat.group}))
-    dfs <- vapply(results_grouped, function(x) {round(lavaan::lavInspect(x, what = "test")[[testName]]$df, 3)}, 0)
 
-    aic <- vapply(results_grouped, AIC, 0)
-    bic <- vapply(results_grouped, BIC, 0)
-    Ns  <- vapply(results_grouped, lavaan::lavInspect, 0, what = "ntotal")
+    chiSq <- lapply(semResults, function(x) {lavaan::lavInspect(x, what = "test")[[testName]]$stat.group})
+    logLGroup <- unlist(lapply(semResults, function(x) x@loglik$loglik.group))
+    npar <- unlist(lapply(semResults, function(x) x@loglik$npar))
+    aics <- -2 * logLGroup + 2 * npar
+    Ns  <- unlist(lapply(semResults, function(x) x@Data@nobs))
+    bics <- -2 * logLGroup + npar * log(Ns)
 
     dtFillGroup <- data.frame(matrix(ncol = 0, nrow = length(models)))
 
     dtFillGroup[["Model"]]    <- models
     dtFillGroup[["group"]]    <- rep(groupNames, length(rownames(lrt$value)))
-    dtFillGroup[["AIC"]]      <- aic
-    dtFillGroup[["BIC"]]      <- bic
+    dtFillGroup[["AIC"]]      <- aics
+    dtFillGroup[["BIC"]]      <- bics
     dtFillGroup[["N"]]        <- Ns
     dtFillGroup[["Chisq"]]    <- chiSq
-    dtFillGroup[["Df"]]       <- dfs
-    dtFillGroup[["PrChisq"]]  <- pchisq(q = chiSq, df = dfs, lower.tail = FALSE)
+    dtFillGroup[["Df"]]       <- NA
+    dtFillGroup[["PrChisq"]]  <- NA
 
+    # we want the LRT for multiple models
     if (length(semResults) > 1) {
-      dchisq   <- rep(NA, length(groupNames))
-      ddf      <- rep(NA, length(groupNames))
-      dPrChisq <- rep(NA, length(groupNames))
-      for(i in 1:(length(chiSq)-length(groupNames))) {
-        lrt <- lavaan::lavTestLRT(results_grouped[[i]], results_grouped[[i+length(groupNames)]])
-        dchisq <- c(dchisq, ifelse(lrt[2, "Chisq diff"] < 1e-05, 0, lrt[2, "Chisq diff"]))
-        ddf    <- c(ddf, lrt[2, "Df diff"])
-        dPrChisq <- c(dPrChisq, lrt[2, "Pr(>Chisq)"])
-      }
-      dtFillGroup[["dchisq"]] <- dchisq
-      dtFillGroup[["ddf"]] <- ddf
-      dtFillGroup[["dPrChisq"]] <- dPrChisq
+
+      # lrts <- -2 * log()
+
+      # dtFillGroup[["dchisq"]] <- dchisq
+      # dtFillGroup[["ddf"]] <- ddf
+      # dtFillGroup[["dPrChisq"]] <- dPrChisq
 
     }
     dtFill[["group"]] <- gettext("all")
     dtFill <- dtFill[, c(1, ncol(dtFill), 2:(ncol(dtFill)-1))]
-    dtFill <- rbind(dtFill, NA, dtFillGroup)
+    dtFill <- rbind(dtFill, dtFillGroup)
 
   }
 

diff --git a/tests/testthat/test-sem.R b/tests/testthat/test-sem.R
@@ -149,22 +149,6 @@ test_that("Multigroup, multimodel SEM works", {
                                       3, 2), "Model fit table")
 
 
-
-  grouptab <- results[["results"]][["modelContainer"]][["collection"]][["modelContainer_grouptab"]][["data"]]
-  expect_equal_tables(grouptab, list(1554.81154499834, 1622.4700973294, 51.6035277328317, 35, "default",
-                                     37, 0.0348827227017026, "", "", "", 1, 1634.45537215556, 1703.23399086407,
-                                     34.0761536514708, 35, "default", 38, 0.512542051610495, "",
-                                     "", "", 2, 1549.89265811707, 1612.7184567102, 53.4279668719255,
-                                     38, "constrained", 37, 0.0495817824685969, 0.647670229628093,
-                                     1.65194318186774, 3, 1, 1634.45537215827, 1703.23399086677,
-                                     34.5269699000817, 35, "constrained", 38, 0.490791865857937,
-                                     "", 0, 0, 2, 1549.89265811344, 1612.71845670656, 54.5525537466226,
-                                     38, "more constrained", 37, 0.0399564191888851, "", 0, 0, 1,
-                                     1631.17917549795, 1695.04503572728, 38.1908170729371, 38, "more constrained",
-                                     38, 0.460812777513164, 0.110723819122933, 6.0182887677031, 3,
-                                     2), "Model fit by group table")
-
-
   rsquared <- results[["results"]][["modelContainer"]][["collection"]][["modelContainer_rsquared"]][["data"]]
   expect_equal_tables(rsquared, list(1, "x1", 0.883076871616545, 0.88344099961725, 0.883440941017356,
                                      1, "x2", 0.993698159869737, 0.993307380054294, 0.993308312239008,