Skip to content

Commit

Permalink
added fit by group table
Browse files Browse the repository at this point in the history
  • Loading branch information
@LSLindeloo
LSLindeloo committed Feb 17, 2023
1 parent 0f0ce71 commit 6daf3fe
Show file tree
Hide file tree
Showing 2 changed files with 147 additions and 16 deletions.
137 changes: 127 additions & 10 deletions R/sem.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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"))
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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")
Expand All @@ -462,33 +462,76 @@ 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("&#967;&sup2;"), 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!
semResults <- .semComputeResults(modelContainer, dataset, options)

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 {
Ns <- vapply(semResults, lavaan::lavInspect, 0, what = "ntotal")
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)"]]
Expand Down Expand Up @@ -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) {
Expand Down
26 changes: 20 additions & 6 deletions tests/testthat/test-sem.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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"]]
Expand Down

0 comments on commit 6daf3fe

Please sign in to comment.