From 6daf3fe21057e29b4f633a85af983ef504466176 Mon Sep 17 00:00:00 2001 From: Lorenzo Lindeloo <94008175+LSLindeloo@users.noreply.github.com> Date: Fri, 17 Feb 2023 12:44:42 +0100 Subject: [PATCH] added fit by group table --- R/sem.R | 137 +++++++++++++++++++++++++++++++++++--- tests/testthat/test-sem.R | 26 ++++++-- 2 files changed, 147 insertions(+), 16 deletions(-) diff --git a/R/sem.R b/R/sem.R index fd2c12c1..5e65000e 100644 --- a/R/sem.R +++ b/R/sem.R @@ -203,7 +203,7 @@ checkLavaanModel <- function(model, availableVars) { modelContainer <- createJaspContainer() modelContainer$dependOn(c("samplingWeights", "meanStructure", "manifestInterceptFixedToZero", "latentInterceptFixedToZero", "exogenousCovariateFixed", "orthogonal", "factorScaling", "residualSingleIndicatorOmitted", "residualVariance", "exogenousLatentCorrelation", - "dependentCorrelation", "threshold", "scalingParameter", "efaConstrained", "standardizedVariable", "naAction", "estimator", "test", + "dependentCorrelation", "threshold", "scalingParameter", "efaConstrained", "standardizedVariable", "naAction", "estimator", "modelTest", "errorCalculationMethod", "informationMatrix", "emulation", "group", "equalLoading", "equalIntercept", "equalResidual", "equalResidualCovariance", "equalMean", "equalThreshold", "equalRegression", "equalVariance", "equalLatentCovariance", "dataType", "sampleSize", "freeParameters")) @@ -360,11 +360,11 @@ checkLavaanModel <- function(model, availableVars) { lavopts[["estimator"]] <- options[["estimator"]] lavopts[["se"]] <- ifelse(options[["errorCalculationMethod"]] == "bootstrap", "standard", options[["errorCalculationMethod"]]) lavopts[["information"]] <- options[["informationMatrix"]] - lavopts[["test"]] <- ifelse(options[["modelTest"]] == "satorraBentler", "Satorra.Bentler", - ifelse(options[["modelTest"]] == "yuanBentler", "Yuan.Bentler", + lavopts[["test"]] <- ifelse(options[["modelTest"]] == "satorraBentler", "satorra.bentler", + ifelse(options[["modelTest"]] == "yuanBentler", "yuan.bentler", ifelse(options[["modelTest"]] == "meanAndVarianceAdjusted", "mean.var.adjusted", ifelse(options[["modelTest"]] == "scaledAndShifted", "scaled.shifted", - ifelse(options[["modelTest"]] == "bollenStine", "Bollen.Stine", + ifelse(options[["modelTest"]] == "bollenStine", "bollen.stine", options[["modelTest"]]))))) # group.equal options @@ -443,7 +443,7 @@ checkLavaanModel <- function(model, availableVars) { fittab$dependOn("models") fittab$position <- 0 - fittab$addColumnInfo(name = "Model", title = "", type = "string" ) + fittab$addColumnInfo(name = "Model", title = "", type = "string" , combine = TRUE) fittab$addColumnInfo(name = "AIC", title = gettext("AIC"), type = "number" ) fittab$addColumnInfo(name = "BIC", title = gettext("BIC"), type = "number" ) fittab$addColumnInfo(name = "N", title = gettext("n"), type = "integer") @@ -462,6 +462,35 @@ checkLavaanModel <- function(model, availableVars) { modelContainer[["fittab"]] <- fittab + if (options[["group"]] != "") { + grouptab <- createJaspTable(title = gettext("Model fit by group")) + grouptab$dependOn("models") + grouptab$position <- 0.05 + + grouptab$addColumnInfo(name = "Model", title = "", type = "string" , combine = TRUE) + grouptab$addColumnInfo(name = "group", title = gettext("Group"), type = "string" ) + grouptab$addColumnInfo(name = "AIC", title = gettext("AIC"), type = "number" ) + grouptab$addColumnInfo(name = "BIC", title = gettext("BIC"), type = "number" ) + grouptab$addColumnInfo(name = "N", title = gettext("n"), type = "integer") + grouptab$addColumnInfo(name = "Chisq", title = gettext("χ²"), type = "number" , + overtitle = gettext("Baseline test")) + grouptab$addColumnInfo(name = "Df", title = gettext("df"), type = "integer", + overtitle = gettext("Baseline test")) + grouptab$addColumnInfo(name = "PrChisq", title = gettext("p"), type = "pvalue", + overtitle = gettext("Baseline test")) + if (length(options[["models"]]) > 1) { + grouptab$addColumnInfo(name = "dchisq", title = gettext("\u0394\u03C7\u00B2"), type = "number" , + overtitle = gettext("Difference test")) + grouptab$addColumnInfo(name = "ddf", title = gettext("\u0394df"), type = "integer", + overtitle = gettext("Difference test")) + grouptab$addColumnInfo(name = "dPrChisq", title = gettext("p"), type = "pvalue" , + overtitle = gettext("Difference test")) + } + + modelContainer[["grouptab"]] <- grouptab + } + + if (!ready) return() # add data to the table! @@ -469,8 +498,19 @@ checkLavaanModel <- function(model, availableVars) { if (modelContainer$getError()) return() + testName <- switch(options[["modelTest"]], + "satorraBentler" = "satorra.bentler", + "yuanBentler" = "yuan.bentler", + "scaledAndShifted" = "scaled.shifted", + "meanAndVarianceAdjusted" = "mean.var.adjusted", + "bollenStine" = "bollen.stine", + options[["modelTest"]]) + if (testName == "default") + testName <- "standard" if (length(semResults) == 1) { - lrt <- .withWarnings(lavaan::lavTestLRT(semResults[[1]])[-1, ]) + lrt <- .withWarnings(lavaan::lavTestLRT(semResults[[1]], type = "Chisq")[-1, ]) + chiSq <- lavaan::lavInspect(semResults[[1]], what = "test")[[testName]]$stat + dfs <- lavaan::lavInspect(semResults[[1]], what = "test")[[testName]]$df rownames(lrt$value) <- options[["models"]][[1]][["name"]] Ns <- lavaan::lavInspect(semResults[[1]], "ntotal") } else { @@ -478,17 +518,20 @@ checkLavaanModel <- function(model, availableVars) { lrt_args <- semResults names(lrt_args) <- "object" # (the first result is object, the others ...) lrt_args[["model.names"]] <- vapply(options[["models"]], getElement, name = "name", "") + lrt_args[["type"]] <- "Chisq" lrt <- .withWarnings(do.call(lavaan::lavTestLRT, lrt_args)) lrt$value[1,5:7] <- NA + chiSq <- unlist(lapply(semResults, function(x) {lavaan::lavInspect(x, what = "test")[[testName]]$stat})) + dfs <- unlist(lapply(semResults, function(x) {round(lavaan::lavInspect(x, what = "test")[[testName]]$df, 3)})) } fittab[["Model"]] <- rownames(lrt$value) fittab[["AIC"]] <- lrt$value[["AIC"]] fittab[["BIC"]] <- lrt$value[["BIC"]] fittab[["N"]] <- Ns - fittab[["Chisq"]] <- lrt$value[["Chisq"]] - fittab[["Df"]] <- lrt$value[["Df"]] - fittab[["PrChisq"]] <- pchisq(q = lrt$value[["Chisq"]], df = lrt$value[["Df"]], lower.tail = FALSE) + fittab[["Chisq"]] <- chiSq + fittab[["Df"]] <- dfs + fittab[["PrChisq"]] <- pchisq(q = chiSq, df = dfs, lower.tail = FALSE) fittab[["dchisq"]] <- lrt$value[["Chisq diff"]] fittab[["ddf"]] <- lrt$value[["Df diff"]] fittab[["dPrChisq"]] <- lrt$value[["Pr(>Chisq)"]] @@ -526,13 +569,87 @@ checkLavaanModel <- function(model, availableVars) { "boot", gettext("bootstrap (Bollen-Stine) probability value") ) testname <- LUT[test == tolower(LUT$option), "name"][[1]] - ftext <- gettextf("Model tests based on %s.", testname) + if (length(semResults) == 1) + ftext <- gettextf("Baseline tests based on %s.", testname) + if (length(semResults) > 1) + ftext <- gettextf("Baseline tests based on %s. Difference tests based on a function of two standard test-statistics.", testname) + fittab$addFootnote(message = ftext) } if (options$estimator %in% c("dwls", "gls", "wls", "uls")) { fittab$addFootnote(message = gettext("The AIC, BIC and additional information criteria are only available with ML-type estimators")) } + if (options[["group"]] != "") { + + 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 + fit_group <- try(do.call(lavaan::lavaan, lav_args)) + if (isTryError(fit_group)) { + errormsg <- gettextf("The model fit by group is unavailable for the specified model: %s", options[["models"]][[i]][["name"]]) + grouptab$setError(errormsg) + break + } + results_grouped[[k]] <- fit_group + k = k + 1 + } + } + + 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") + + grouptab[["Model"]] <- models + grouptab[["group"]] <- rep(groupNames, length(rownames(lrt$value))) + grouptab[["AIC"]] <- aic + grouptab[["BIC"]] <- bic + grouptab[["N"]] <- Ns + grouptab[["Chisq"]] <- chiSq + grouptab[["Df"]] <- dfs + grouptab[["PrChisq"]] <- pchisq(q = chiSq, df = dfs, lower.tail = FALSE) + 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)"]) + } + grouptab[["dchisq"]] <- dchisq + grouptab[["ddf"]] <- ddf + grouptab[["dPrChisq"]] <- dPrChisq + + } + if (test != "standard") { + grouptab$addFootnote(message = ftext) + } + } } .semParameters <- function(modelContainer, dataset, options, ready) { diff --git a/tests/testthat/test-sem.R b/tests/testthat/test-sem.R index e92400b1..350bd517 100644 --- a/tests/testthat/test-sem.R +++ b/tests/testthat/test-sem.R @@ -118,12 +118,26 @@ test_that("Multigroup, multimodel SEM works", { results <- jaspTools::runAnalysis("SEM", "poldem_grouped.csv", options) fittab <- results[["results"]][["modelContainer"]][["collection"]][["modelContainer_fittab"]][["data"]] - expect_equal_tables(fittab, list(3189.26691715402, 3383.93591869107, 82.2032643259552, 70, "default", - 75, 0.150923770163707, "", "", "", 3184.34803034567, 3372.06456754211, - 83.2843775176126, 73, "constrained", 75, 0.19249046571302, 0.647754490401111, - 1.65156784895969, 3, 3181.07183366569, 3361.83590652152, 86.0081808376312, - 76, "more constrained", 75, 0.202652084622712, 0.110596895607335, - 6.02091896557529, 3), "Model fit table") + expect_equal_tables(fittab, list(3189.26691715402, 3383.93591869107, 85.6796813843025, 70, "default", + 75, 0.0980338951401389, "", "", "", 3184.34803034567, 3372.06456754211, + 87.9549367720072, 73, "constrained", 75, 0.111927575441429, + 0.647754490401136, 1.65156784895958, 3, 3181.07183366569, 3361.83590652152, + 92.7433708195597, 76, "more constrained", 75, 0.0929761753671129, + 0.110596895607244, 6.02091896557718, 3), "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"]]