diff --git a/DESCRIPTION b/DESCRIPTION index 9a996e6a..b155f220 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -12,16 +12,18 @@ Encoding: UTF-8 Imports: forcats, ggplot2, - jaspBase, - jaspGraphs, lavaan, cSEM, reshape2, + jaspBase, + jaspGraphs, semPlot, semTools, stringr, tibble, - tidyr + tidyr, + SEMsens Remotes: jasp-stats/jaspBase, jasp-stats/jaspGraphs + diff --git a/R/common.R b/R/common.R index 73a60fcd..28fa0068 100644 --- a/R/common.R +++ b/R/common.R @@ -127,3 +127,435 @@ lavBootstrap <- function(fit, samples = 1000) { return(semPlotMod) } +.sa.aco <- function (data = NULL, sample.cov, sample.nobs, model, sens.model, + opt.fun, d = NULL, paths = NULL, verbose = TRUE, max.value = Inf, + max.iter = 1000, e = 1e-10, n.of.ants = 10, k = 100, q = 1e-04, + sig.level = 0.05, rate.of.conv = 0.1, measurement = FALSE, + xi = 0.5, seed = NULL, ...) { + for.n.of.sens.pars <- lavaan::lavaanify(sens.model, fixed.x = TRUE) + n.of.sens.pars <- length(for.n.of.sens.pars[which(for.n.of.sens.pars$lhs != + "phantom" & for.n.of.sens.pars$rhs == "phantom"), ]$lhs) + if (n.of.sens.pars < 2) + stop("Sensitivity model must have at least two sensitivity parameters or phantom coefficients.") + if (is.null(data)) { + old.out = lavaan::sem(model = model, sample.cov = sample.cov, + sample.nobs = sample.nobs, ...) + } + else { + old.out = lavaan::sem(model = model, data = data, ...) + } + old.par = lavaan::standardizedSolution(old.out, type = "std.all") + old.fit <- lavaan::fitMeasures(old.out) + if (!is.null(seed)) { + set.seed(seed) + } + if (is.null(paths)) { + paths <- old.par + } + if (is.character(paths)) { + paths <- lavaan::lavaanify(paths, fixed.x = TRUE) + } + e.abs <- e + e.rel <- e + eval <- 0 + iter <- 0 + last.impr <- max.iter + nl <- matrix(NA, k, k - 1) + sens.pars <- data.frame() + outcome <- vector() + model.results <- data.frame() + max.X <- rep(NA, n.of.sens.pars) + max.y <- -Inf + p.X <- vector() + sens.fit <- vector() + p <- data.frame(v = numeric(), sd = numeric(), gr = numeric()) + if (is.null(d)) { + d <- list(rep(c(-1, 1), n.of.sens.pars)) + } + else { + if (!is.list(d)) + stop("d (domain) must be in a list format; e.g.,\n d = list(-1, 1,\n -1, 1,\n -1, 1,\n -1, 1)") + } + if (rate.of.conv <= 0 | rate.of.conv > 1) + stop("Convergence rate (rate.of.conv) must be in (0, 1]") + for (i in 1:(round(1/rate.of.conv * k, 0))) { + X <- vector() + for (j in 1:n.of.sens.pars) { + X <- c(X, stats::runif(1, d[[1]][2 * j - 1], d[[1]][2 * + j])) + } + X <- t(X) + new.model = sens.model + for (l in 1:n.of.sens.pars) { + new.model = gsub(paste("phantom", l, sep = ""), paste(X[l]), + new.model, ignore.case = FALSE, perl = FALSE, + fixed = FALSE, useBytes = FALSE) + } + iter <- iter + 1 + if((2 * k) < max.iter) + progressbarTick() + + warnings <- options(warn = 2) + if (is.null(data)) { + new.out = try(lavaan::sem(model = new.model, sample.cov = sample.cov, + sample.nobs = sample.nobs, ...), silent = TRUE) + } + else { + new.out = try(lavaan::sem(model = new.model, data = data, + ...), silent = TRUE) + } + if (isTRUE(class(new.out) == "try-error")) { + next + } + on.exit(options(warnings)) + new.par = lavaan::standardizedSolution(new.out, type = "std.all") + eval <- eval + 1 + if((2 * k) >= max.iter) + progressbarTick() + + new.par$lines <- 1:length(new.par[, 1]) + new.par$evals <- eval + model.results <- rbind(model.results, new.par) + if (eval == 1) { + sens.out <- new.out + model.1 <- model.results + model.1$path <- paste(model.1$lhs, model.1$op, model.1$rhs, + sep = "") + phan.names <- model.1[which(model.1$evals == 1 & + model.1$op == "~" & model.1$rhs == "phantom"), + ]$path + if (is.data.frame(paths)) { + if (measurement) { + paths <- which(model.1$lhs %in% paths$lhs & + model.1$rhs %in% paths$rhs) + } + else { + paths <- which(model.1$lhs %in% paths$lhs & + model.1$op == "~" & model.1$rhs %in% paths$rhs) + } + } + } + sens.par <- c(X, eval = eval) + sens.pars <- rbind(sens.pars, sens.par) + fit <- c(lavaan::fitMeasures(new.out), eval = eval) + sens.fit <- rbind(sens.fit, fit) + if (!is.numeric(opt.fun)) { + y <- eval(opt.fun) + } + else if (opt.fun == 1) { + y <- mean(abs(old.par$est[paths]), na.rm = TRUE)/mean(abs(new.par$est[paths]), + na.rm = TRUE) + } + else if (opt.fun == 2) { + y <- stats::sd(new.par$est[paths] - old.par$est[paths], + na.rm = TRUE)/mean(abs(old.par$est[paths]), na.rm = TRUE) + } + else if (opt.fun == 3) { + y <- mean(abs(new.par$pvalue[paths] - old.par$pvalue[paths]), + na.rm = TRUE) + } + else if (opt.fun == 4) { + y <- 1/mean(abs(new.par$pvalue[paths] - rep(sig.level, + length(paths))), na.rm = TRUE) + } + else if (opt.fun == 5) { + y <- abs(unname(lavaan::fitmeasures(new.out)["rmsea"]) - + unname(lavaan::fitmeasures(old.out)["rmsea"])) + } + else if (opt.fun == 6) { + y <- 1/abs(unname(lavaan::fitmeasures(new.out)["rmsea"]) - + 0.05) + } + outcome <- c(outcome, y) + p.X <- rbind(p.X, X) + p <- rbind(p, data.frame(v = y, sd = 0, gr = 0)) + if (eval == k) { + break + } + } + if (length(p.X) == 0 | length(p$v) < k) + .quitAnalysis("Sensitivity analysis models do not reach the specified convergence rate.\n Please set a lower convergence rate threshhold or reduce model complexicity") + p$gr <- rank(-p$v, ties.method = "random") + for (i in 1:k) { + nl[i, ] <- (1:k)[1:k != i] + } + while (TRUE) { + dist.mean <- p.X + if (sum(apply(dist.mean, 2, stats::sd)) == 0) { + colnames(sens.pars) <- c(phan.names, "eval") + return(list(n.eval = eval, n.iter = iter, max.y = max.y, + phantom.coef = max.X, old.model.par = old.par, + old.model.fit = old.fit, model = model, sens.model = sens.model, + sens.fit = sens.fit, outcome = outcome, sens.pars = sens.pars, + model.results = model.results, old.out = old.out, + sens.out = sens.out)) + } + dist.rank <- p$gr + dim(dist.mean) <- c(length(p$v), n.of.sens.pars) + o.X <- vector() + o.X <- SEMsens::gen.sens.pars(dist.mean, dist.rank, n.of.ants, + nl, q, k, xi) + if (length(o.X) == 0) { + colnames(sens.pars) <- c(phan.names, "eval") + return(list(n.eval = eval, n.iter = iter, max.y = max.y, + phantom.coef = max.X, old.model.par = old.par, + old.model.fit = old.fit, model = model, sens.model = sens.model, + sens.fit = sens.fit, outcome = outcome, sens.pars = sens.pars, + model.results = model.results, old.out = old.out, + sens.out = sens.out)) + } + X <- o.X + dim(X) <- c(length(X)/n.of.sens.pars, n.of.sens.pars) + for (j in 1:dim(X)[1]) { + X.sens <- X[j, ] + X.model <- as.vector(X.sens) + new.model = sens.model + for (i in 1:dim(X)[2]) { + new.model = gsub(paste("phantom", i, sep = ""), + paste(X.model[i]), new.model, ignore.case = FALSE, + perl = FALSE, fixed = FALSE, useBytes = FALSE) + } + iter <- iter + 1 + if((2 * k) < max.iter) + progressbarTick() + + warnings <- options(warn = 2) + on.exit(options(warnings)) + if (is.null(data)) { + new.out = try(lavaan::sem(model = new.model, + sample.cov = sample.cov, sample.nobs = sample.nobs, + ...), TRUE) + } + else { + new.out = try(lavaan::sem(model = new.model, + data = data, ...), TRUE) + } + if (isTRUE(class(new.out) != "try-error")) { + new.par <- lavaan::standardizedSolution(new.out, + type = "std.all") + eval <- eval + 1 + if((2 * k) >= max.iter) + progressbarTick() + + p.X <- rbind(p.X, X.sens) + new.par$lines <- 1:length(new.par[, 1]) + new.par$evals <- eval + model.results <- rbind(model.results, new.par) + fit <- c(lavaan::fitMeasures(new.out), eval = eval) + sens.fit <- rbind(sens.fit, fit) + sens.par <- c(X.sens, eval = eval) + sens.pars <- rbind(sens.pars, sens.par) + if (!is.numeric(opt.fun)) { + y <- eval(opt.fun) + } + else if (opt.fun == 1) { + y <- mean(abs(old.par$est[paths]), na.rm = TRUE)/mean(abs(new.par$est[paths]), + na.rm = TRUE) + } + else if (opt.fun == 2) { + y <- stats::sd(new.par$est[paths] - old.par$est[paths], + na.rm = TRUE)/mean(abs(old.par$est[paths]), + na.rm = TRUE) + } + else if (opt.fun == 3) { + y <- mean(abs(new.par$pvalue[paths] - old.par$pvalue[paths]), + na.rm = TRUE) + } + else if (opt.fun == 4) { + y <- 1/mean(abs(new.par$pvalue[paths] - rep(sig.level, + length(paths))), na.rm = TRUE) + } + else if (opt.fun == 5) { + y <- abs(unname(lavaan::fitmeasures(new.out)["rmsea"]) - + unname(lavaan::fitmeasures(old.out)["rmsea"])) + } + else if (opt.fun == 6) { + y <- abs(unname(lavaan::fitmeasures(new.out)["rmsea"]) - + 0.05) + } + outcome <- c(outcome, y) + p <- rbind(p, data.frame(v = y, sd = 0, gr = 0)) + p$gr <- rank(-p$v, ties.method = "random") + idx.final <- p$gr <= k + p <- p[idx.final, ] + p.X <- p.X[idx.final, ] + dim(p.X) <- c(length(p.X)/n.of.sens.pars, n.of.sens.pars) + } + } + p$gr <- rank(-p$v, ties.method = "random") + for (i in 1:k) { + nl[i, ] <- (1:k)[1:k != i] + } + if (max(outcome, na.rm = TRUE) > max.y) { + max.y <- max(outcome, na.rm = TRUE) + max.X <- sens.pars[which.max(outcome), ] + colnames(max.X) <- c(phan.names, "eval") + last.impr <- eval + } + if ((abs(max.y - max.value) < abs(e.rel * max.value + + e.abs)) | (max.y > max.value)) { + colnames(sens.pars) <- c(phan.names, "eval") + return(list(n.eval = eval, n.iter = iter, max.y = max.y, + phantom.coef = max.X, old.model.par = old.par, + old.model.fit = old.fit, model = model, sens.model = sens.model, + sens.fit = sens.fit, outcome = outcome, sens.pars = sens.pars, + model.results = model.results, old.out = old.out, + sens.out = sens.out)) + } + if (max.iter > 0 & iter >= max.iter) { + colnames(sens.pars) <- c(phan.names, "eval") + return(list(n.eval = eval, n.iter = iter, max.y = max.y, + phantom.coef = max.X, old.model.par = old.par, + old.model.fit = old.fit, model = model, sens.model = sens.model, + sens.fit = sens.fit, outcome = outcome, sens.pars = sens.pars, + model.results = model.results, old.out = old.out, + sens.out = sens.out)) + } + } +} + +.sa.tabu <- function (model, sens.model, data = NULL, sample.cov = NULL, + sample.nobs = NULL, opt.fun = 1, sig.level = 0.05, ...) { + init.model <- model + init.model.par.table <- lavaan::lavaanify(init.model, auto = T, + model.type = "sem", fixed.x = TRUE) + non.phan.path.ids <- which(init.model.par.table$op == "~") + non.phan.path.names <- character(length(non.phan.path.ids)) + for (i in seq_along(non.phan.path.ids)) { + j <- non.phan.path.ids[i] + non.phan.path.names[i] <- paste(init.model.par.table$lhs[j], + init.model.par.table$op[j], init.model.par.table$rhs[j]) + } + sens.model.par.table <- lavaan::lavaanify(sens.model, auto = T, + model.type = "sem", fixed.x = TRUE) + phan.path.ids <- which(sens.model.par.table$label != "") + phan.path.names <- character(length(phan.path.ids)) + for (i in seq_along(phan.path.ids)) { + j <- phan.path.ids[i] + phan.path.names[i] <- paste(sens.model.par.table$lhs[j], + sens.model.par.table$op[j], sens.model.par.table$rhs[j]) + } + init.model.sem <- lavaan::sem(model = init.model.par.table, + data = data, sample.cov = sample.cov, sample.nobs = sample.nobs) + init.model.params <- lavaan::standardizedSolution(init.model.sem, + type = "std.all") + sens.model.template <- sens.model + f <- function(phantom.coef) { + for (j in 1:length(phantom.coef)) { + sens.model.template <- gsub(paste0("phantom", j), + paste(phantom.coef[j]), sens.model.template) + } + sens.model.template.par.table <- lavaan::lavaanify(sens.model.template, + auto = T, model.type = "sem", fixed.x = TRUE) + sens.model.sem <- try(lavaan::sem(model = sens.model.template.par.table, + data = data, sample.cov = sample.cov, sample.nobs = sample.nobs), + silent = TRUE) + sens.model.params <- lavaan::standardizedSolution(sens.model.sem, + type = "std.all") + if (opt.fun == 1) { + y <- mean(abs(sens.model.params$est[non.phan.path.ids] - + init.model.params$est[non.phan.path.ids]), na.rm = TRUE)/mean(abs(init.model.params$est[non.phan.path.ids]), + na.rm = TRUE) + } + else if (opt.fun == 2) { + y <- stats::sd(sens.model.params$est[non.phan.path.ids] - + init.model.params$est[non.phan.path.ids], na.rm = TRUE)/mean(abs(init.model.params$est[non.phan.path.ids]), + na.rm = TRUE) + } + else if (opt.fun == 3) { + y <- mean(abs(sens.model.params$pvalue[non.phan.path.ids] - + init.model.params$pvalue[non.phan.path.ids]), + na.rm = TRUE) + } + else if (opt.fun == 4) { + y <- mean(abs(sens.model.params$pvalue[non.phan.path.ids] - + rep(sig.level, length(non.phan.path.ids))), na.rm = TRUE) + } + else if (opt.fun == 5) { + y <- abs(unname(lavaan::fitmeasures(sens.model.sem)["rmsea"]) - + unname(lavaan::fitmeasures(init.model.sem)["rmsea"])) + } + else if (opt.fun == 6) { + y <- abs(unname(lavaan::fitmeasures(sens.model.sem)["rmsea"]) - + 0.05) + } + return(list(y = y, model = sens.model.params)) + } + res <- .sa.tabu.helper(length(phan.path.ids), f, maximum = TRUE, + ...) + colnames(res$best.param) <- phan.path.names + out <- list(model = model, old.model.par = init.model.params, + model.results = res$model.history, best.param = res$best.param[1, + ], best.obj = res$best.obj, sens.par = NULL, outcome = NULL) + return(out) +} + +.sa.tabu.helper <- function (n.var, f, maximum = FALSE, max.len = 1, max.tabu.size = 5, + neigh.size = NULL, max.iter = NULL, max.iter.obj = NULL, + range = c(-1, 1), r = 1e-05, verbose = TRUE, seed = NULL) +{ + if (is.null(neigh.size)) { + neigh.size <- min(n.var * 2, 10) + } + if (is.null(max.iter)) { + max.iter <- n.var * 50 + } + if (is.null(max.iter.obj)) { + max.iter.obj <- n.var * 5 + } + options(warn = 2) + tabu.list <- list() + n.iter <- 1 + n.iter.obj <- 1 + model.history <- list() + max.attempts <- 50 + if (!is.null(seed)) { + set.seed(seed) + } + for (i in 1:max.attempts) { + best.param <- current.param <- t(stats::runif(n.var, + -1, 1)) + best.obj <- try(f(best.param), silent = TRUE) + if (class(best.obj)[1] != "try-error") { + break + } + } + if (class(best.obj)[1] == "try-error") { + .quitAnalysis("Can't find a valid set of initial parameters for the sensitivity analysis! Maybe try a different seed?") + } + best.obj <- best.obj$y + tabu.list[[1]] <- current.param + if (verbose) { + cat(" n curr_obj best_obj\n") + } + while ((n.iter <= max.iter) & (n.iter.obj <= max.iter.obj)) { + best.neighbor <- SEMsens::gen.neighbors.tabu(current.param, maximum, + neigh.size, tabu.list, max.len, range, r, f) + current.param <- best.neighbor$best.param + current.obj <- best.neighbor$best.obj + best.neighbor$best.model$evals <- n.iter + model.history[[n.iter]] <- best.neighbor$best.model + if ((maximum & current.obj > best.obj) || (!maximum & + current.obj < best.obj)) { + best.obj <- current.obj + best.param <- current.param + n.iter.obj <- 1 + } + else { + n.iter.obj <- n.iter.obj + 1 + } + tabu.list <- append(tabu.list, list(current.param)) + if (length(tabu.list) > max.tabu.size) { + tabu.list <- tabu.list[-1] + } + if (verbose) { + cat(sprintf("%3d %10f %10f\n", n.iter, current.obj, + best.obj)) + } + n.iter <- n.iter + 1 + progressbarTick() + } + return(list(best.param = best.param, best.obj = best.obj, + model.history = do.call(rbind, model.history))) +} + diff --git a/R/sem.R b/R/sem.R index 492b1f7e..02deb6b3 100644 --- a/R/sem.R +++ b/R/sem.R @@ -36,6 +36,7 @@ SEMInternal <- function(jaspResults, dataset, options, ...) { .semParameters(modelContainer, dataset, options, ready) .semAdditionalFits(modelContainer, dataset, options, ready) .semRsquared(modelContainer, dataset, options, ready) + .semSensitivity(modelContainer, dataset, options, ready) .semMardiasCoefficient(modelContainer, dataset, options, ready) .semCov(modelContainer, dataset, options, ready) .semMI(modelContainer, datset, options, ready) @@ -207,7 +208,7 @@ checkLavaanModel <- function(model, availableVars) { "dependentCorrelation", "threshold", "scalingParameter", "efaConstrained", "standardizedVariable", "naAction", "estimator", "test", "errorCalculationMethod", "informationMatrix", "emulation", "group", "equalLoading", "equalIntercept", "equalResidual", "equalResidualCovariance", "equalMean", "equalThreshold", "equalRegression", - "equalVariance", "equalLatentCovariance", "dataType", "sampleSize", "freeParameters", "manifestMeanFixedToZero")) + "equalVariance", "equalLatentCovariance", "dataType", "sampleSize", "freeParameters", "sensitivityAnalysis", "manifestMeanFixedToZero")) jaspResults[["modelContainer"]] <- modelContainer } @@ -1325,6 +1326,290 @@ checkLavaanModel <- function(model, availableVars) { } } +.semSensitivity <- function(modelContainer, dataset, options, ready) { + if (!options[["sensitivityAnalysis"]] || !is.null(modelContainer[["sensitivity"]])) return() + + sensitivity <- createJaspContainer(gettext("Sensitivity analysis")) + sensitivity$position <- 0.8 + sensitivity$dependOn(c("sensitivityAnalysis", "searchAlgorithm", "optimizerFunction", "sizeOfSolutionArchive", "numberOfAnts", "alpha", "maxIterations", "setSeed", "seed", "models")) + + modelContainer[["sensitivity"]] <- sensitivity + + if (length(options[["models"]]) < 2) { + .semSensitivityTables(modelContainer[["results"]][["object"]][[1]], NULL, sensitivity, dataset, options, ready) + } else { + + for (i in seq_along(options[["models"]])) { + fit <- modelContainer[["results"]][["object"]][[i]] + model <- options[["models"]][[i]] + .semSensitivityTables(fit, model, sensitivity, dataset, options, ready) + } + } +} + +.semSensitivityTables <- function(fit, model, parentContainer, dataset, options, ready) { + if (is.null(model)) { + sencont <- parentContainer + } else { + sencont <- createJaspContainer(model[["name"]], initCollapsed = TRUE) + } + + + # Summary of sensitivity analysis + sensumtab <- createJaspTable(title = gettext("Summary of sensitivity analysis")) + + if (options[["group"]] != "") + sensumtab$addColumnInfo(name = "group", title = gettext("Group"), type = "string", combine = TRUE) + + sensumtab$addColumnInfo(name = "path", title = gettext("Path"), type = "string") + sensumtab$addColumnInfo(name = "est", title = gettext("Standardized estimate"), overtitle = gettext("Original model"), type = "number") + sensumtab$addColumnInfo(name = "pvalue", title = gettext("p"), overtitle = gettext("Original model"), type = "pvalue") + sensumtab$addColumnInfo(name = "pvaluesens", title = gettext("p\u002A"), overtitle = gettext("Sensitivity model"), type = "pvalue") + sensumtab$addColumnInfo(name = "mean", title = gettext("Mean"), overtitle = gettext("Sensitivity model"), type = "number") + sensumtab$addColumnInfo(name = "min", title = gettext("Min"), overtitle = gettext("Sensitivity model"), type = "number") + sensumtab$addColumnInfo(name = "max", title = gettext("Max"), overtitle = gettext("Sensitivity model"), type = "number") + + sencont[["sensum"]] <- sensumtab + + if (!ready || !inherits(fit, "lavaan")) return() + + # create SEMsens model + if (is.null(model)) { + analyticModel <- .semTranslateModel(options[["models"]][[1]][["syntax"]], dataset) + } else { + analyticModel <- .semTranslateModel(model[["syntax"]], dataset) + } + #enrich model + modelTable <- lavaan::lavInspect(fit, what = "list") + pathVars <- unique(c(modelTable[modelTable$op == "~", ]$rhs, modelTable[modelTable$op == "~", ]$lhs)) + + if (length(pathVars) < 2) { + .quitAnalysis(gettext("Please include at least one path in the model to perform a sensitivity analysis.")) + } + + sensModel <- analyticModel + for (i in seq_along(pathVars)) { + sensParameter <- paste0("\n", pathVars[i], " ~", " phantom", i, "*phantom\n") + sensModel <- paste0(sensModel, sensParameter) + } + sensModel <- paste0(sensModel, "\nphantom =~ 0\nphantom ~~ 1*phantom\n") + optimizerFunction <- switch(options[["optimizerFunction"]], + "percentChangeMeanEstimate" = 1, + "sdOfDeviance" = 2, + "changeOfPvalue" = 3, + "distanceOfPvalue" = 4, + "changeOfRmsea" = 5, + "distanceOfRmsea" = 6 + ) + + if (options[["group"]] != "") { + saTables <- lapply(1:5, function(x) data.frame()) + for (group in unique(dataset[, options[["group"]]])) { + data <- dataset[dataset[[options[["group"]]]] == group,] + if(options[["searchAlgorithm"]] == "antColonyOptimization") { + iter <- ifelse((2 * options[["sizeOfSolutionArchive"]]) >= options[["maxIterations"]], (options[["sizeOfSolutionArchive"]] + options[["numberOfAnts"]]), (options[["maxIterations"]] + options[["numberOfAnts"]])) + startProgressbar(iter, + sprintf("Performing sensitivity analysis (model: %1$s, group: %2$s)", + ifelse(is.null(model), + options[["models"]][[1]][["name"]], + model[["name"]]), + group)) + sa <- .sa.aco(data = data, model = analyticModel, sens.model = sensModel, n.of.ants = options[["numberOfAnts"]], k = options[["sizeOfSolutionArchive"]], rate.of.conv = options[["convergenceRateThreshold"]], opt.fun = optimizerFunction, max.iter = options[["maxIterations"]], sig.level = options[["alpha"]], seed = if (options[["setSeed"]]) options[["seed"]] else NULL) + } + if(options[["searchAlgorithm"]] == "tabuSearch") { + startProgressbar(options[["maxIterations"]], + sprintf("Performing sensitivity analysis (model: %1$s, group: %2$s)", + ifelse(is.null(model), + options[["models"]][[1]][["name"]], + model[["name"]]), + group)) + sa <- .sa.tabu(data = data, model = analyticModel, sens.model = sensModel, opt.fun = optimizerFunction, max.iter = options[["maxIterations"]], sig.level = options[["alpha"]], seed = if (options[["setSeed"]]) options[["seed"]] else NULL) + } + saTablesRows <- SEMsens::sens.tables(sa) + saTablesRows <- sapply(saTablesRows, function(x) { + cbind(x, "group" = rep(group, length(x[, 1])), "rowname" = row.names(x)) + }) + for (table in seq_along(saTablesRows)) { + saTables[[table]] <- rbind(saTables[[table]], saTablesRows[[table]]) + } + } + saTables <- sapply(saTables, as.data.frame) + saTables <- sapply(saTables, function(x) {x[order(x[["group"]], x[["rowname"]]),]}) + } else { + if (options[["dataType"]] == "raw") { + if(options[["searchAlgorithm"]] == "antColonyOptimization") { + iter <- ifelse((2 * options[["sizeOfSolutionArchive"]]) >= options[["maxIterations"]], (options[["sizeOfSolutionArchive"]] + options[["numberOfAnts"]]), (options[["maxIterations"]] + options[["numberOfAnts"]])) + startProgressbar(iter, + sprintf("Performing sensitivity analysis (model: %1$s)", + ifelse(is.null(model), + options[["models"]][[1]][["name"]], + model[["name"]]) + )) + sa <- .sa.aco(data = dataset, model = analyticModel, sens.model = sensModel, n.of.ants = options[["numberOfAnts"]], k = options[["sizeOfSolutionArchive"]], rate.of.conv = options[["convergenceRateThreshold"]], opt.fun = optimizerFunction, max.iter = options[["maxIterations"]], sig.level = options[["alpha"]], seed = if (options[["setSeed"]]) options[["seed"]] else NULL) + } + if(options[["searchAlgorithm"]] == "tabuSearch") { + startProgressbar(options[["maxIterations"]], + sprintf("Performing sensitivity analysis (model: %1$s)", + ifelse(is.null(model), + options[["models"]][[1]][["name"]], + model[["name"]]) + )) + sa <- .sa.tabu(data = dataset, model = analyticModel, sens.model = sensModel, opt.fun = optimizerFunction, max.iter = options[["maxIterations"]], sig.level = options[["alpha"]], seed = if (options[["setSeed"]]) options[["seed"]] else NULL) + } + } else { + if (is.null(model)) { + syntax <- options[["models"]][[1]][["syntax"]] + } else { + syntax <- model[["syntax"]] + } + dataset <- .semDataCovariance(dataset, syntax) + if(options[["searchAlgorithm"]] == "antColonyOptimization") { + iter <- ifelse((2 * options[["sizeOfSolutionArchive"]]) >= options[["maxIterations"]], (options[["sizeOfSolutionArchive"]] + options[["numberOfAnts"]]), (options[["maxIterations"]] + options[["numberOfAnts"]])) + startProgressbar(iter, + sprintf("Performing sensitivity analysis (model: %1$s)", + ifelse(is.null(model), + options[["models"]][[1]][["name"]], + model[["name"]]) + )) + sa <- .sa.aco(sample.cov = dataset, sample.nobs = options[["sampleSize"]], model = analyticModel, sens.model = sensModel, n.of.ants = options[["numberOfAnts"]], k = options[["sizeOfSolutionArchive"]], rate.of.conv = options[["convergenceRateThreshold"]], opt.fun = optimizerFunction, max.iter = options[["maxIterations"]], sig.level = options[["alpha"]], seed = if (options[["setSeed"]]) options[["seed"]] else NULL) + } + if(options[["searchAlgorithm"]] == "tabuSearch") { + startProgressbar(options[["maxIterations"]], + sprintf("Performing sensitivity analysis (model: %1$s)", + ifelse(is.null(model), + options[["models"]][[1]][["name"]], + model[["name"]]) + )) + sa <- .sa.tabu(sample.cov = dataset, sample.nobs = options[["sampleSize"]], model = analyticModel, sens.model = sensModel, opt.fun = optimizerFunction, max.iter = options[["maxIterations"]], sig.level = options[["alpha"]], seed = if (options[["setSeed"]]) options[["seed"]] else NULL) + } + } + saTables <- SEMsens::sens.tables(sa) + saTables <- sapply(saTables, function(x) { + cbind(x, "rowname" = row.names(x)) + }) + saTables <- sapply(saTables, as.data.frame) + saTables <- sapply(saTables, function(x) { + x[order(x[["rowname"]]),] + }) + } + + + # Fill table + + if (options[["group"]] != "") + sensumtab[["group"]] <- saTables[[1]][["group"]] + + sensumtab[["path"]] <- saTables[[1]][["rowname"]] + sensumtab[["est"]] <- saTables[[1]][["model.est"]] + sensumtab[["pvalue"]] <- saTables[[1]][["model.pvalue"]] + sensumtab[["pvaluesens"]] <- saTables[[5]][["p.changed"]] + sensumtab[["mean"]] <- saTables[[1]][["mean.est.sens"]] + sensumtab[["min"]] <- saTables[[1]][["min.est.sens"]] + sensumtab[["max"]] <- saTables[[1]][["max.est.sens"]] + + + # Sensitivity parameters that led to a change in significance + senpartab <- createJaspTable(title = gettext("Sensitivity parameters that led to a change in significance")) + + if (options[["group"]] != "") + senpartab$addColumnInfo(name = "group", title = gettext("Group"), type = "string", combine = TRUE) + + senpartab$addColumnInfo(name = "path", title = gettext("Path"), type = "string") + + sensitivityParameters <- grep("~", colnames(saTables[[5]]), value = TRUE) + for (par in sensitivityParameters) { + senpartab$addColumnInfo(name = par, title = gettext(par), overtitle = gettext("Sensitivity parameters"), type = "number") + } + + sencont[["senpar"]] <- senpartab + + # Fill table + saTable_clean <- saTables[[5]][!is.na(saTables[[5]][["p.changed"]]), ] + if (nrow(saTable_clean) == 0) sencont[["senpar"]] <- NULL + + + if (options[["group"]] != "") + senpartab[["group"]] <- saTable_clean[["group"]] + + senpartab[["path"]] <- saTable_clean[["rowname"]] + for (par in sensitivityParameters) { + senpartab[[par]] <- saTable_clean[[par]] + } + + # # Sensitivity parameters that led to min est + # senparmintab <- createJaspTable(title = gettext("Sensitivity parameters that led to the minimum estimates in the sensitivity model")) + # + # if (options[["group"]] != "") + # senparmintab$addColumnInfo(name = "group", title = gettext("Group"), type = "string", combine = TRUE) + # + # senparmintab$addColumnInfo(name = "path", title = gettext("Path"), type = "string") + # + # sensitivityParameters <- grep("~", colnames(saTables[[3]]), value = TRUE) + # for (par in sensitivityParameters) { + # senparmintab$addColumnInfo(name = par, title = gettext(par), overtitle = gettext("Sensitivity parameters"), type = "number") + # } + # + # sencont[["senparmin"]] <- senparmintab + # + # # Fill table + # if (options[["group"]] != "") + # senparmintab[["group"]] <- saTables[[3]][["group"]] + # + # senparmintab[["path"]] <- saTables[[3]][["rowname"]] + # for (par in sensitivityParameters) { + # senparmintab[[par]] <- saTables[[3]][[par]] + # } + + # # Sensitivity parameters that led to max est + # senparmaxtab <- createJaspTable(title = gettext("Sensitivity parameters that led to the maximum estimates in the sensitivity model")) + # + # if (options[["group"]] != "") + # senparmaxtab$addColumnInfo(name = "group", title = gettext("Group"), type = "string", combine = TRUE) + # + # senparmaxtab$addColumnInfo(name = "path", title = gettext("Path"), type = "string") + # + # sensitivityParameters <- grep("~", colnames(saTables[[4]]), value = TRUE) + # for (par in sensitivityParameters) { + # senparmaxtab$addColumnInfo(name = par, title = gettext(par), overtitle = gettext("Sensitivity parameters"), type = "number") + # } + # + # sencont[["senparmax"]] <- senparmaxtab + # + # # Fill table + # if (options[["group"]] != "") + # senparmaxtab[["group"]] <- saTables[[4]][["group"]] + # + # senparmaxtab[["path"]] <- saTables[[4]][["rowname"]] + # for (par in sensitivityParameters) { + # senparmaxtab[[par]] <- saTables[[4]][[par]] + # } + + # Summary of sensitivity parameters + sensumpartab <- createJaspTable(title = gettext("Summary of sensitivity parameters")) + + if (options[["group"]] != "") + sensumpartab$addColumnInfo(name = "group", title = gettext("Group"), type = "string", combine = TRUE) + + sensumpartab$addColumnInfo(name = "par", title = gettext("Sensitivity parameter"), type = "string") + sensumpartab$addColumnInfo(name = "mean", title = gettext("Mean"), type = "number") + sensumpartab$addColumnInfo(name = "min", title = gettext("Min"), type = "number") + sensumpartab$addColumnInfo(name = "max", title = gettext("Max"), type = "number") + + sencont[["sensumpar"]] <- sensumpartab + + #Fill table + + if (options[["group"]] != "") + sensumpartab[["group"]] <- saTables[[2]][["group"]] + + sensumpartab[["par"]] <- saTables[[2]][["rowname"]] + sensumpartab[["mean"]] <- saTables[[2]][["mean.phan"]] + sensumpartab[["min"]] <- saTables[[2]][["min.phan"]] + sensumpartab[["max"]] <- saTables[[2]][["max.phan"]] + + if (!is.null(model)) parentContainer[[model[["name"]]]] <- sencont +} + .semMardiasCoefficient <- function(modelContainer, dataset, options, ready) { if (!options[["mardiasCoefficient"]] || !is.null(modelContainer[["semMardiasTable"]])) return() diff --git a/inst/qml/SEM.qml b/inst/qml/SEM.qml index 1b688253..d059b6bb 100644 --- a/inst/qml/SEM.qml +++ b/inst/qml/SEM.qml @@ -94,10 +94,10 @@ Form CheckBox { name: "residualCovariance"; label: qsTr("Residual covariances") } CheckBox { name: "standardizedResidual"; label: qsTr("Standardized residuals") } CheckBox { name: "mardiasCoefficient"; label: qsTr("Mardia's coefficient") } + CheckBox {name: "standardizedEstimate"; label: qsTr("Standardized estimates") } } Group { - CheckBox{name: "standardizedEstimate"; label: qsTr("Standardized estimates"); checked: false} CheckBox { name: "pathPlot"; @@ -132,8 +132,52 @@ Form } } } + CheckBox + { + name: "sensitivityAnalysis" + label: qsTr("Sensitivity analysis") + RadioButtonGroup + { + title: qsTr("Search algorithm") + name: "searchAlgorithm" + id: search + RadioButton + { + value: "antColonyOptimization" + label: qsTr("Ant colony optimization") + checked: true + IntegerField { name: "numberOfAnts"; label: qsTr("Number of ants"); defaultValue: 10 } + IntegerField { name: "sizeOfSolutionArchive"; label: qsTr("Size of the solution archive"); defaultValue: 100 } + DoubleField { name: "convergenceRateThreshold"; label: qsTr("Convergence rate threshold"); defaultValue: 0.1; negativeValues: false } + } + RadioButton { value: "tabuSearch"; label: qsTr("Tabu search") } + } + DropDown + { + name: "optimizerFunction" + label: qsTr("Optimizer function") + values: + [ + { label: qsTr("% change mean estimate") , value: "percentChangeMeanEstimate" }, + { label: qsTr("Sd of deviance / old estimate") , value: "sdOfDeviance" }, + { label: qsTr("Change of p-value") , value: "changeOfPvalue" }, + { label: qsTr("Distance of p-value from alpha") , value: "distanceOfPvalue" }, + { label: qsTr("Change of RMSEA") , value: "changeOfRmsea" }, + { label: qsTr("Distance of RMSEA from 0.05") , value: "distanceOfRmsea" } + ] + } + DoubleField + { + name: "alpha" + label: qsTr("Significance level") + negativeValues: false + decimals: 4 + defaultValue: 0.05 + } + IntegerField { name: "maxIterations"; label: qsTr("Maximum number of iterations"); defaultValue: 1000 } + SetSeed{} + } } - } Section diff --git a/tests/testthat/test-sem.R b/tests/testthat/test-sem.R index f3f44c22..3aa95158 100644 --- a/tests/testthat/test-sem.R +++ b/tests/testthat/test-sem.R @@ -41,6 +41,48 @@ test_that("Basic SEM works", { "(Residual) variances parameter table") }) +test_that("Sensitivity analysis works", { + options <- analysisOptions("SEM") + options$samplingWeights <- "" + options$sensitivityAnalysis <- TRUE + options$sizeOfSolutionArchive <- 10 + options$maxIterations <- 10 + options$informationMatrix <- "expected" + options$estimator <- "default" + options$modelTest <- "standard" + options$emulation <- "lavaan" + options$group <- "" + options$setSeed <- TRUE + options$seed <- 1 + options$models <- list(list(name = "Model1", syntax = list(model = "\n # latent variable definitions\n ind60 =~ x1 + x2 + x3\n dem60 =~ y1 + y2 + y3 + y4\n dem65 =~ y5 + y6 + y7 + y8\n # regressions\n dem60 ~ ind60\n dem65 ~ ind60 + dem60\n # residual (co)variances\n y1 ~~ y5\n y2 ~~ y4 + y6\n y3 ~~ y7\n y4 ~~ y8\n y6 ~~ y8\n ", + columns = c("x1", "x2", "x3", "y1", "y2", "y3", "y4", "y5", "y6", "y7", "y8")))) + set.seed(1) + results <- runAnalysis("SEM", "C:/JASP/Development/Modules/jaspSem/tests/testthat/poldem_grouped.csv", options) + + sencont <- results[["results"]][["modelContainer"]][["collection"]][["modelContainer_sensitivity"]][["collection"]] + sensumtab <- sencont[["modelContainer_sensitivity_sensum"]][["data"]] + senpartab <- sencont[["modelContainer_sensitivity_senpar"]][["data"]] + sensumpartab <- sencont[["modelContainer_sensitivity_sensumpar"]][["data"]] + + expect_equal_tables(sensumtab, list(0.446712905437343, 0.817304450826043, 0.454361526415421, 0.109312503661981, + "dem60~ind60", 1.54185554406272e-05, 0.397240001459216, 0.885228943292037, + 0.940534819295127, 0.889796135759017, 0.838319341744439, "dem65~dem60", + 0, "", 0.1822594099762, 0.335620470155875, 0.172277112909256, + -0.13098545021451, "dem65~ind60", 0.00967601956304076, 0.125746802148124 + )) + + + expect_equal_tables(senpartab, list(0.401998712582474, -0.133707492109712, 0.839306883347651, "dem60~ind60", + 0.136451103548822, -0.356128009622868, -0.859719801919832, "dem65~ind60" + )) + + + expect_equal_tables(sensumpartab, list(0.401998712582474, 0.113806982564325, -0.428828062076269, "dem60~phantom", + 0.21276746153651, -0.00984860889287904, -0.356128009622868, + "dem65~phantom", 0.87878147725911, 0.249441092202601, -0.937148975256221, + "ind60~phantom")) +}) + test_that("Multigroup, multimodel SEM works", { options <- jaspTools::analysisOptions("SEM") options$emulation = "lavaan"