From a640959b1dc3875d36590a78fbc87a6b0344584f Mon Sep 17 00:00:00 2001 From: Julius Pfadt Date: Thu, 15 Aug 2024 16:35:50 +0200 Subject: [PATCH] add boostrapped CIs for standardized estimates --- R/common.R | 41 +++++++++++- R/sem.R | 30 ++++++--- inst/qml/SEM.qml | 87 ++++++++++--------------- tests/testthat/test-mediationanalysis.R | 82 +++++++++++------------ tests/testthat/test-sem.R | 53 ++++++++------- 5 files changed, 166 insertions(+), 127 deletions(-) diff --git a/R/common.R b/R/common.R index 11bca257..2410ceab 100644 --- a/R/common.R +++ b/R/common.R @@ -18,7 +18,7 @@ # Function commonly used in the various procedures within the SEM module -lavBootstrap <- function(fit, samples = 1000) { +lavBootstrap <- function(fit, samples = 1000, standard = FALSE, typeStd = NULL) { # Run bootstrap, track progress with progress bar # Notes: faulty runs are simply ignored # recommended: add a warning if not all boot samples are successful @@ -38,14 +38,51 @@ lavBootstrap <- function(fit, samples = 1000) { return(lavaan::coef(lav_object)) } + + coef_with_callback_std <- function(lav_object, typeStd) { + std <- lavaan::standardizedSolution(lav_object, type = typeStd) + out <- std$est.std + + progressbarTick() + + return(out) + } + startProgressbar(samples + 1) - bootres <- lavaan::bootstrapLavaan(object = fit, R = samples, FUN = coef_with_callback) + if (!standard) { + bootres <- lavaan::bootstrapLavaan(object = fit, R = samples, FUN = coef_with_callback) + } else { + bootres <- lavaan::bootstrapLavaan(object = fit, R = samples, FUN = coef_with_callback_std, typeStd = typeStd) + } # Add the bootstrap samples to the fit object fit@boot <- list(coef = bootres) fit@Options$se <- "bootstrap" + # exclude error bootstrap runs + err_id <- attr(fit@boot$coef, "error.idx") + if (length(err_id) > 0L) { + fit@boot$coef <- fit@boot$coef[-err_id, , drop = FALSE] + } + + # we actually need the SEs from the bootstrap not the SEs from ML or something + N <- nrow(fit@boot$coef) + + # we multiply the var by (n-1)/n because lavaan actually uses n for the variance instead of n-1 + if (!standard) { + # for unstandardized + fit@ParTable$se[fit@ParTable$free != 0] <- apply(fit@boot$coef, 2, sd) * sqrt((N-1)/N) + } else { + fit@ParTable$se <- apply(fit@boot$coef, 2, sd) * sqrt((N-1)/N) + # the standardized solution gives all estimates not only the unconstrained, so we need to change + # the free prameters in the partable and also change the estimate + fit@ParTable$free <- seq_len(ncol(fit@boot$coef)) + std <- lavaan::standardizedSolution(fit, type = typeStd) + fit@ParTable$est <- std$est.std + } + + return(fit) } diff --git a/R/sem.R b/R/sem.R index 7bd6848f..fd875a97 100644 --- a/R/sem.R +++ b/R/sem.R @@ -155,12 +155,12 @@ SEMInternal <- function(jaspResults, dataset, options, ...) { } } - # check if standardized estimates are required but standard errors are bootstrapped, that does not work yet, - # because the CIs are not really available in lavaan - if (options[["standardizedEstimate"]] && options[["errorCalculationMethod"]] == "bootstrap") { - .quitAnalysis(gettext("Standardized estimates are not available when bootstrapping is used for standard error calculation.")) - return() - } + # # check if standardized estimates are required but standard errors are bootstrapped, that does not work yet, + # # because the CIs are not really available in lavaan + # if (options[["standardizedEstimate"]] && options[["errorCalculationMethod"]] == "bootstrap") { + # .quitAnalysis(gettext("Standardized estimates are not available when bootstrapping is used for standard error calculation.")) + # return() + # } return() } @@ -340,7 +340,15 @@ checkLavaanModel <- function(model, availableVars) { } if (options[["errorCalculationMethod"]] == "bootstrap") { - fit$value <- lavBootstrap(fit$value, options[["bootstrapSamples"]]) + type <- switch(options[["standardizedEstimateType"]], + "all" = "std.all", + "latents" = "std.lv", + "nox" = "std.nox") + fit$value <- lavBootstrap(fit$value, samples = options[["bootstrapSamples"]], + standard = options[["standardizedEstimate"]], + typeStd = type) + modelContainer$dependOn(optionsFromObject = modelContainer, + options = c("bootstrapSamples", "standardizedEstimate", "standardizedEstimateType")) } results[[i]] <- fit$value @@ -1137,10 +1145,15 @@ checkLavaanModel <- function(model, availableVars) { ifelse(options[["bootstrapCiType"]] == "percentile", "perc", "norm")) - if (!options[["standardizedEstimate"]]) { + #' we need the second option in the if statement because when we require standardized estimates and bootstrapped CIs + #' the standardization happens in each bootstrap run and the standardized estimates replace the regular raw estimates + #' in the fit object, and we only need to call parameterEstimates + if (!options[["standardizedEstimate"]] || + (options[["standardizedEstimate"]] && options[["errorCalculationMethod"]] == "bootstrap")) { pe <- lavaan::parameterestimates(fit, level = options[["ciLevel"]], boot.ci.type = bootstrapCiType) pe <- lavaan::lavMatrixRepresentation(lavaan::lav_partable_complete(pe)) + } else { type <- switch(options[["standardizedEstimateType"]], "all" = "std.all", @@ -1149,6 +1162,7 @@ checkLavaanModel <- function(model, availableVars) { pe <- lavaan::standardizedSolution(fit, level = options[["ciLevel"]], type = type) pe <- lavaan::lavMatrixRepresentation(lavaan::lav_partable_complete(pe)) colnames(pe)[colnames(pe) == "est.std"] <- "est" + } # TODO: bootstrap standardized CI: diff --git a/inst/qml/SEM.qml b/inst/qml/SEM.qml index 1a6cb254..87382069 100644 --- a/inst/qml/SEM.qml +++ b/inst/qml/SEM.qml @@ -135,23 +135,23 @@ Form label: qsTr("Estimator") id: estimator values: [ - { value: "default", label: qsTr("Default")}, - { value: "ml", label: qsTr("ML") }, - { value: "gls", label: qsTr("GLS") }, - { value: "wls", label: qsTr("WLS") }, - { value: "uls", label: qsTr("ULS") }, - { value: "dwls", label: qsTr("DWLS") }, - { value: "dls", label: qsTr("DLS") }, - { value: "pml", label: qsTr("PML") }, - { value: "mlm", label: qsTr("MLM") }, - { value: "mlmv", label: qsTr("MLMV") }, - { value: "mlmvs", label: qsTr("MLMVS") }, - { value: "mlf", label: qsTr("MLF") }, - { value: "mlr", label: qsTr("MLR") }, - { value: "wlsm", label: qsTr("WLSM") }, - { value: "wlsmv", label: qsTr("WLSMV") }, - { value: "ulsm", label: qsTr("ULSM") }, - { value: "ulsmv", label: qsTr("ULSMV") } + { label: qsTr("Default"), value: "default" }, + { label: qsTr("ML"), value: "ml" }, + { label: qsTr("GLS"), value: "gls" }, + { label: qsTr("WLS"), value: "wls" }, + { label: qsTr("ULS"), value: "uls" }, + { label: qsTr("DWLS"), value: "dwls" }, + { label: qsTr("DLS"), value: "dls" }, + { label: qsTr("PML"), value: "pml" }, + { label: qsTr("MLM"), value: "mlm" }, + { label: qsTr("MLMV"), value: "mlmv" }, + { label: qsTr("MLMVS"), value: "mlmvs" }, + { label: qsTr("MLF"), value: "mlf" }, + { label: qsTr("MLR"), value: "mlr" }, + { label: qsTr("WLSM"), value: "wlsm" }, + { label: qsTr("WLSMV"), value: "wlsmv" }, + { label: qsTr("ULSM"), value: "ulsm" }, + { label: qsTr("ULSMV"), value: "ulsmv" } ] } HelpButton @@ -167,16 +167,16 @@ Form name: "modelTest" label: qsTr("Model test") values: [ - { value: "default", label: qsTr("Default") }, - { value: "standard", label: qsTr("Standard") }, - { value: "satorraBentler", label: qsTr("Satorra-Bentler") }, - { value: "yuanBentler", label: qsTr("Yuan-Bentler") }, - { value: "yuanBentlerMplus", label: qsTr("Yuan-Bentler Mplus") }, - { value: "meanAndVarianceAdjusted", label: qsTr("Mean and variance adjusted") }, - { value: "scaledAndShifted", label: qsTr("Scaled and shifted") }, - { value: "bollenStine", label: qsTr("Bootstrap (Bollen-Stine)") }, - { value: "browneResidualAdf", label: qsTr("Browne residual based (ADF)")}, - { value: "browneResidualNt", label: qsTr("Browne residual based (NT)") } + { label: qsTr("Default"), value: "default" }, + { label: qsTr("Standard"), value: "standard" }, + { label: qsTr("Satorra-Bentler"), value: "satorraBentler" }, + { label: qsTr("Yuan-Bentler"), value: "yuanBentler" }, + { label: qsTr("Yuan-Bentler Mplus"), value: "yuanBentlerMplus" }, + { label: qsTr("Mean and variance adjusted"), value: "meanAndVarianceAdjusted" }, + { label: qsTr("Scaled and shifted"), value: "scaledAndShifted" }, + { label: qsTr("Bootstrap (Bollen-Stine)"), value: "bollenStine" }, + { label: qsTr("Browne residual based (ADF)"), value: "browneResidualAdf" }, + { label: qsTr("Browne residual based (NT)"), value: "browneResidualNt" } ] } @@ -185,10 +185,10 @@ Form label: qsTr("Information matrix") name: "informationMatrix" values: [ - { value: "default", label: qsTr("Default") }, - { value: "expected", label: qsTr("Expected") }, - { value: "observed", label: qsTr("Observed") }, - { value: "firstOrder", label: qsTr("First order") } + { label: qsTr("Default"), label: "default" }, + { label: qsTr("Expected"), value: "expected" }, + { label: qsTr("Observed"), value: "observed" }, + { label: qsTr("First order"), value: "firstOrder" } ] } @@ -198,11 +198,11 @@ Form name: "errorCalculationMethod" id: errorCalc values: [ - { value: "default", label: qsTr("Default") }, - { value: "standard", label: qsTr("Standard") }, - { value: "robust", label: qsTr("Robust") }, - { value: "robustHuberWhite", label: qsTr("Robust Huber-White") }, - { value: "bootstrap", label: qsTr("Bootstrap") } + { label: qsTr("Default"), label: "default" }, + { label: qsTr("Standard"), value: "standard" }, + { label: qsTr("Robust"), value: "robust" }, + { label: qsTr("Robust Huber-White"), value: "robustHuberWhite" }, + { label: qsTr("Bootstrap"), value: "bootstrap" } ] } IntegerField @@ -238,21 +238,6 @@ Form id: missingG CheckBox{name: "standardizedVariable"; label: qsTr("Standardize variables before estimation"); checked: false} - // property var withFiml: [ - // { label: qsTr("Listwise deletion") , value: "listwise" }, - // { label: qsTr("FIML") , value: "fiml"}, - // { label: qsTr("Pairwise") , value: "pairwise" }, - // { label: qsTr("Two-stage") , value: "twoStage" }, - // { label: qsTr("Robust two-stage") , value: "twoStageRobust" }, - // { label: qsTr("Doubly robust") , value: "doublyRobust" } - // ] - // property var noFiml: [ - // { label: qsTr("Listwise deletion") , value: "listwise" }, - // { label: qsTr("Pairwise") , value: "pairwise" }, - // { label: qsTr("Two-stage") , value: "twoStage" }, - // { label: qsTr("Robust two-stage") , value: "twoStageRobust" }, - // { label: qsTr("Doubly robust") , value: "doublyRobust" } - // ] DropDown { name: "naAction" diff --git a/tests/testthat/test-mediationanalysis.R b/tests/testthat/test-mediationanalysis.R index 6408cd9d..3d72b8c8 100644 --- a/tests/testthat/test-mediationanalysis.R +++ b/tests/testthat/test-mediationanalysis.R @@ -170,53 +170,53 @@ test_that("Multiple mediation with missing values works", { }) -test_that("Bootstrapping works", { - options <- jaspTools::analysisOptions("MediationAnalysis") - options$predictors <- "contcor1" - options$mediators <- "contcor2" - options$outcomes <- "contNormal" - options$emulation <- "lavaan" - options$estimator <- "ml" - options$errorCalculationMethod <- "bootstrap" - options$bootstrapSamples <- 100 - options$bootstrapCiType <- "percentileBiasCorrected" - options$naAction <- "fiml" - - set.seed(1) - results <- jaspTools::runAnalysis("MediationAnalysis", "test.csv", options) - - # Direct effects table results match +# bootstrapped mediation analysis works +options <- jaspTools::analysisOptions("MediationAnalysis") +options$predictors <- "contcor1" +options$mediators <- "contcor2" +options$outcomes <- "contNormal" +options$emulation <- "lavaan" +options$estimator <- "ml" +options$errorCalculationMethod <- "bootstrap" +options$bootstrapSamples <- 100 +options$bootstrapCiType <- "percentileBiasCorrected" +options$naAction <- "fiml" + +set.seed(1) +results <- jaspTools::runAnalysis("MediationAnalysis", "test.csv", options) + +test_that("Direct effects table results match", { table <- results[["results"]][["modelContainer"]][["collection"]][["modelContainer_parest"]][["collection"]][["modelContainer_parest_dir"]][["data"]] jaspTools::expect_equal_tables(table, - list(0.0761400441341873, 0.671443275931252, 0.257757857221231, "contcor1", - "", 0.0585458466298006, "contNormal", 0.136265300259528, - 1.89158837011558), - label = "Direct effects table results match") + list(0.0761400441342019, 0.671443275931248, 0.257757843176866, "contcor1", + "", 0.0589199802205698, "contNormal", 0.136467182830779, + 1.88878994810415)) +}) - # Indirect effects table results match +test_that("Indirect effects table results match", { table <- results[["results"]][["modelContainer"]][["collection"]][["modelContainer_parest"]][["collection"]][["modelContainer_parest_ind"]][["data"]] jaspTools::expect_equal_tables(table, - list(-0.299720825954288, 0.07329367848577, -0.0893136271813779, "contcor2", - "", "", 0.32161886441617, 0.0901123252539836, - "contcor1", "contNormal", -0.991136638963043), - label = "Indirect effects table results match") + list(-0.299720825954233, 0.073293678485771, -0.0893136104463747, "contcor2", + "", "", 0.321618931630391, 0.0901123208859354, + "contcor1", "contNormal", -0.991136501294072)) +}) - # Total effects table results match - table <- results[["results"]][["modelContainer"]][["collection"]][["modelContainer_parest"]][["collection"]][["modelContainer_parest_tot"]][["data"]] - jaspTools::expect_equal_tables(table, - list(-0.00126714302931763, 0.440665080716848, 0.168444230039853, "contcor1", - "", 0.102761018524938, "contNormal", 0.103237851474511, - 1.63161309184588), - label = "Total effects table results match") - # Path coefficients table results match +test_that("Path coefficients table results match", { table <- results[["results"]][["modelContainer"]][["collection"]][["modelContainer_parest"]][["collection"]][["modelContainer_parest_path"]][["data"]] jaspTools::expect_equal_tables(table, - list(-0.475862353202891, 0.110585516049068, -0.136980438302071, "contcor2", - "", 0.318469030743472, "contNormal", 0.137308753597124, - -0.997608926696571, 0.0761400441341873, 0.671443275931252, 0.257757844059186, - "contcor1", "", 0.0585458547695155, "contNormal", 0.136265297698809, - 1.89158830907129, 0.48165177793888, 0.754639896102569, 0.652017217213314, - "contcor1", "", 0, "contcor2", 0.0748152807686884, - 8.71502733818777), - label = "Path coefficients table results match") + list(-0.475862353202892, 0.110585521838546, -0.136980447123533, "contcor2", + "", 0.306661780077382, "contNormal", 0.133998745278152, + -1.02225171466488, 0.0761400441342019, 0.671443275931248, 0.257757843176866, + "contcor1", "", 0.0589199802205698, "contNormal", 0.136467182830779, + 1.88878994810415, 0.481651777938874, 0.754639896102554, 0.652017220865318, + "contcor1", "", 0, "contcor2", 0.06513006081573, 10.0110027950081 + )) +}) + +test_that("Total effects table results match", { + table <- results[["results"]][["modelContainer"]][["collection"]][["modelContainer_parest"]][["collection"]][["modelContainer_parest_tot"]][["data"]] + jaspTools::expect_equal_tables(table, + list(-0.00126714302930483, 0.440665080716848, 0.168444232730491, "contcor1", + "", 0.102761008345243, "contNormal", 0.103237850067971, + 1.63161314013793)) }) diff --git a/tests/testthat/test-sem.R b/tests/testthat/test-sem.R index 73178dc6..cb6cb434 100644 --- a/tests/testthat/test-sem.R +++ b/tests/testthat/test-sem.R @@ -615,44 +615,47 @@ test_that("Bootstrapping model fit table works", { 75, 0, 5, 5)) }) -# Residual covariances table results match -test_that("Bootstrapping residual covariances work", { +test_that("Residual covariances table results match", { table <- results[["results"]][["modelContainer"]][["collection"]][["modelContainer_params"]][["collection"]][["modelContainer_params_cov"]][["data"]] jaspTools::expect_equal_tables(table, - list(1.782009228184, 1.782009228184, 1.782009228184, "", "x2 - x3", + list(1.78200922818399, 1.78200922818399, 1.78200922818399, "", "x2 - x3", "", 0, "", 1.2564136096, 1.2564136096, 1.2564136096, "", "x2 - y1", - "", 0, "", 0.899301179999998, 0.899301179999998, 0.899301179999998, - "", "x3 - y1", "", 0, ""), - label = "Residual covariances table results match") + "", 0, "", 0.899301179999995, 0.899301179999995, 0.899301179999995, + "", "x3 - y1", "", 0, "")) }) - -# Regression coefficients table results match -test_that("Bootstrapping regression coefficients work", { +test_that("Regression coefficients table results match", { table <- results[["results"]][["modelContainer"]][["collection"]][["modelContainer_params"]][["collection"]][["modelContainer_params_reg"]][["data"]] jaspTools::expect_equal_tables(table, - list(0.249609316197927, 0.447867152439518, 0.355178154406208, "x1", - 3.13082892944294e-14, "x2", 0.0467810932812294, 7.59234403247094, - -0.0367153099343007, 0.19997796482335, 0.0779182667327594, "x1", - 0.112759784761359, "x3", 0.0491315895805092, 1.58590974560428, - 0.00301905711375345, 0.056769954444189, 0.0306885894512192, - "x1", 0.0358943944804584, "y1", 0.014626696431825, 2.09812171834281 - ), - label = "Regression coefficients table results match") + list(0.24960931619795, 0.447867152439228, 0.355178159312481, "x1", + 4.22839541158737e-12, "x2", 0.051257217491476, 6.92932969628222, + -0.0367153099339115, 0.199977964917624, 0.0779182588187343, + "x1", 0.152733815036115, "x3", 0.0544905227320188, 1.42994148178632, + 0.00301905711375213, 0.0567699544441715, 0.0306885890333793, + "x1", 0.0306723196916876, "y1", 0.0141991468224129, 2.16129809890678 + )) }) +test_that("Total effects table results match", { + table <- results[["results"]][["modelContainer"]][["collection"]][["modelContainer_params"]][["collection"]][["modelContainer_params_toteff"]][["data"]] + jaspTools::expect_equal_tables(table, + list(0.24960931619795, 0.447867152439228, 0.355178159312481, " x2 x1", + 3.13082892944294e-14, 0.046781092952712, 7.59234419066496, -0.0367153099339115, + 0.199977964917624, 0.0779182588187343, " x3 x1", 0.112759818779606, + 0.0491315892354857, 1.5859095956631, 0.00301905711375213, 0.0567699544441715, + 0.0306885890333793, " y1 x1", 0.0358943957021149, + 0.0146266963291098, 2.09812170450981)) +}) -# Residual variances table results match -test_that("Bootstrapping residual variances work", { +test_that("Residual variances table results match", { table <- results[["results"]][["modelContainer"]][["collection"]][["modelContainer_params"]][["collection"]][["modelContainer_params_var"]][["data"]] jaspTools::expect_equal_tables(table, - list(0.0579281075433616, 0.121472727074924, "x1", 0.097380852992327, - "", "x1", 9.14129882900738e-10, 0.0159022267032121, 6.12372435695794, + list(0.0579281075433116, 0.121472727075618, "x1", 0.0973808537831735, + "", "x1", 4.33176605696417e-10, 0.0156019590066371, 6.24157862110443, 2.25167664969695, 2.25167664969695, "x2", 2.25167664969695, - "", "x2", "", 0, "", 1.94967853807201, 1.94967853807201, "x3", - 1.94967853807201, "", "x3", "", 0, "", 6.78685155555555, 6.78685155555555, - "y1", 6.78685155555555, "", "y1", "", 0, ""), - label = "Residual variances table results match") + "", "x2", "", 0, "", 1.949678538072, 1.949678538072, "x3", 1.949678538072, + "", "x3", "", 0, "", 6.78685155555556, 6.78685155555556, "y1", + 6.78685155555556, "", "y1", "", 0, "")) })