Skip to content

Commit

Permalink
add boostrapped CIs for standardized estimates
Browse files Browse the repository at this point in the history
  • Loading branch information
@juliuspfadt
juliuspfadt committed Aug 19, 2024
1 parent 0c64851 commit a640959
Show file tree
Hide file tree
Showing 5 changed files with 166 additions and 127 deletions.
41 changes: 39 additions & 2 deletions R/common.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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)
}

Expand Down
30 changes: 22 additions & 8 deletions R/sem.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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()
}
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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",
Expand 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:
Expand Down
87 changes: 36 additions & 51 deletions inst/qml/SEM.qml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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" }
]
}

Expand All @@ -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" }
]
}

Expand All @@ -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
Expand Down Expand Up @@ -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"
Expand Down
82 changes: 41 additions & 41 deletions tests/testthat/test-mediationanalysis.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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",
"<unicode>", 0.0585458466298006, "contNormal", 0.136265300259528,
1.89158837011558),
label = "Direct effects table results match")
list(0.0761400441342019, 0.671443275931248, 0.257757843176866, "contcor1",
"<unicode>", 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",
"<unicode>", "<unicode>", 0.32161886441617, 0.0901123252539836,
"contcor1", "contNormal", -0.991136638963043),
label = "Indirect effects table results match")
list(-0.299720825954233, 0.073293678485771, -0.0893136104463747, "contcor2",
"<unicode>", "<unicode>", 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",
"<unicode>", 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",
"<unicode>", 0.318469030743472, "contNormal", 0.137308753597124,
-0.997608926696571, 0.0761400441341873, 0.671443275931252, 0.257757844059186,
"contcor1", "<unicode>", 0.0585458547695155, "contNormal", 0.136265297698809,
1.89158830907129, 0.48165177793888, 0.754639896102569, 0.652017217213314,
"contcor1", "<unicode>", 0, "contcor2", 0.0748152807686884,
8.71502733818777),
label = "Path coefficients table results match")
list(-0.475862353202892, 0.110585521838546, -0.136980447123533, "contcor2",
"<unicode>", 0.306661780077382, "contNormal", 0.133998745278152,
-1.02225171466488, 0.0761400441342019, 0.671443275931248, 0.257757843176866,
"contcor1", "<unicode>", 0.0589199802205698, "contNormal", 0.136467182830779,
1.88878994810415, 0.481651777938874, 0.754639896102554, 0.652017220865318,
"contcor1", "<unicode>", 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",
"<unicode>", 0.102761008345243, "contNormal", 0.103237850067971,
1.63161314013793))
})
Loading

0 comments on commit a640959

Please sign in to comment.