Add residual plots

NEWS.md

@@ -16,7 +16,13 @@
 
 # jaspProcess (development version)
 
+## Added
+* Residual plots for each model (`Residuals vs. predicted` and `Q-Q plot standardized residuals`), shown for every dependent variable (the outcome and each mediator), to check the homoscedasticity and normality assumptions of the regression equations. Available in the classical and Bayesian analysis, per model or for all models at once. As in the other JASP analyses, the Q-Q plot can show a confidence band (classical) or the credible interval of each residual (Bayesian).
+* `Set seed` option (Options > Repeatability) in the classical analysis, so bootstrapped standard errors and confidence intervals, and bootstrapped local tests, can be reproduced.
+
 ## Fixed
+* Changing the number of bootstrap replications in the classical analysis now refits the models; before, the results of the previous number of replications were kept.
+* Requesting output for a single model (e.g., a plot or a parameter estimates table) no longer recomputes the model summary and R² tables. In the Bayesian analysis, this avoids recomputing WAIC and LOO for every model.
 * Process models with a categorical predictor no longer error when a factor level label contains a space or other special character (e.g. "emotion regulation"). Such labels are now coerced to syntactically valid names before being used in mean-centering and lavaan model syntax ([Issue #3597](https://github.com/jasp-stats/jasp-issues/issues/3597)).
 
 ---

R/bayesianProcess.R

@@ -64,6 +64,8 @@ BayesianProcess <- function(jaspResults, dataset = NULL, options) {
     .procBayesParameterEstimateTables(parEstContainer, options, modelsContainer)
 
     .procBayesMcmcPlots(jaspResults, options, modelsContainer)
+    # Create residual plots for each dependent variable of each model if requested
+    .procResidualPlots(jaspResults, options, modelsContainer)
 
     # Create html output with lavaan syntax for each model
     .procPlotSyntax(jaspResults, options, modelsContainer)
@@ -216,7 +218,10 @@ BayesianProcess <- function(jaspResults, dataset = NULL, options) {
   tableRowIsValid <- modelNumberIsValid & resultIsValid
 
   summaryTable <- createJaspTable(title = gettext("Model summary"), rowNames = modelNames[tableRowIsValid])
-  summaryTable$dependOn(c(.procGetDependencies(), "processModels", "bicWeights", "hayesNumber"))
+  summaryTable$dependOn(
+    options = c(.procGetDependencies(), "bicWeights", "hayesNumber"),
+    nestedOptions = .procGetAllModelsDependencies(options)
+  )
   summaryTable$position <- 1
 
   ovtWaic <- gettext("WAIC")
@@ -671,6 +676,58 @@ BayesianProcess <- function(jaspResults, dataset = NULL, options) {
   }
 }
 
+.procBayesStandardizedResiduals <- function(fittedModel, dep, ciLevel = NULL) {
+  parTbl <- lavaan::parameterTable(fittedModel)
+  # Only extract the draws of the model parameters: the Stan output also contains, e.g., the log-likelihood of each observation
+  draws  <- as.matrix(fittedModel@external$mcmcout, pars = unique(sub("\\[.*$", "", na.omit(parTbl$pxnames))))
+
+  return(.procBayesStandardizedResidualsFromDraws(parTbl, lavaan::lavInspect(fittedModel, "data"), draws, dep, ciLevel))
+}
+
+.procBayesStandardizedResidualsFromDraws <- function(parTbl, data, draws, dep, ciLevel = NULL, blockSize = 1000) {
+  # Posterior mean (and central credible interval) of the standardized residual of each observation in the regression
+  # equation of dep: for each draw, the observed minus the predicted value, divided by the residual standard deviation
+  parDraws <- function(rows) {
+    # Draws x parameters; fixed parameters have the same value in every draw
+    matrix(vapply(rows, function(r) {
+      if (is.na(parTbl$pxnames[r])) rep(parTbl$est[r], nrow(draws)) else draws[, parTbl$pxnames[r]]
+    }, numeric(nrow(draws))), nrow = nrow(draws))
+  }
+
+  regRows <- which(parTbl$lhs == dep & parTbl$op == "~")
+  preds   <- data[, parTbl$rhs[regRows], drop = FALSE]
+  beta    <- parDraws(regRows)
+  intRow  <- which(parTbl$lhs == dep & parTbl$op == "~1")
+
+  # Without a mean structure the intercepts are implied by the sample means (as in .procResiduals)
+  if (length(intRow) == 0)
+    intercept <- mean(data[, dep], na.rm = TRUE) - drop(beta %*% colMeans(preds, na.rm = TRUE))
+  else
+    intercept <- drop(parDraws(intRow))
+
+  residualSd <- sqrt(drop(parDraws(which(parTbl$lhs == dep & parTbl$op == "~~" & parTbl$rhs == dep))))
+  complete   <- which(stats::complete.cases(data[, dep], preds))
+  probs      <- if (!is.null(ciLevel)) c((1 - ciLevel) / 2, (1 + ciLevel) / 2)
+
+  # Observations x draws matrices, in blocks of observations to limit memory use for large data sets
+  blocks <- split(complete, ceiling(seq_along(complete) / blockSize))
+  blockSummaries <- lapply(blocks, function(obs) {
+    predicted    <- preds[obs, , drop = FALSE] %*% t(beta) + rep(intercept, each = length(obs))
+    standardized <- (data[obs, dep] - predicted) / rep(residualSd, each = length(obs))
+    blockSummary <- data.frame(mean = rowMeans(standardized))
+
+    if (!is.null(probs)) {
+      bounds             <- apply(standardized, 1, stats::quantile, probs = probs, names = FALSE)
+      blockSummary$lower <- bounds[1, ]
+      blockSummary$upper <- bounds[2, ]
+    }
+
+    return(blockSummary)
+  })
+
+  return(do.call(rbind, unname(blockSummaries)))
+}
+
 # Footnotes ----
 
 .procBayesBadWaicFootnote <- function() gettext("Warning: WAIC estimate unreliable -- at least one effective parameter estimate (p_waic) larger than 0.4. We recommend using LOO instead.")

R/classicProcess.R

@@ -70,6 +70,8 @@ ClassicProcess <- function(jaspResults, dataset = NULL, options) {
     localTestContainer <- .procContainerLocalTests(jaspResults, options)
     # Create tables with local test results
     .procLocalTestTables(localTestContainer, dataset, options, modelsContainer)
+    # Create residual plots for each dependent variable of each model if requested
+    .procResidualPlots(jaspResults, options, modelsContainer)
     # Create html output with lavaan syntax for each model
     .procPlotSyntax(jaspResults, options, modelsContainer)
   }
@@ -81,8 +83,8 @@ ClassicProcess <- function(jaspResults, dataset = NULL, options) {
 .procGetDependencies <- function() {
   return(c(
     "dependent", "covariates", "factors", "naAction", "emulation", "estimator",
-    "meanCenteredModeration", "standardizedModelEstimates", "errorCalculationMethod", "bootstrapCiType",
-    "mcmcBurnin", "mcmcSamples", "mcmcChains", "seed", "setSeed", "nuPriorMu",
+    "meanCenteredModeration", "standardizedModelEstimates", "errorCalculationMethod", "bootstrapSamples",
+    "bootstrapCiType", "mcmcBurnin", "mcmcSamples", "mcmcChains", "seed", "setSeed", "nuPriorMu",
     "nuPriorSigma", "betaPriorMu", "betaPriorSigma", "psiPriorAlpha",
     "psiPriorBeta", "rhoPriorAlpha", "rhoPriorBeta", "moderationProbes",
     "moderationProbeType", "moderationProbesMeanSD"
@@ -138,6 +140,10 @@ ClassicProcess <- function(jaspResults, dataset = NULL, options) {
         nestedOptions = .procGetSingleModelsDependencies(as.character(i))[-1] # omit model name
       )
       modelsContainer[[modelName]] <- container
+      # Tables that combine all models only depend on the options of the models they were created for (see
+      # .procGetAllModelsDependencies), so a model that was added has to be included by recreating them
+      jaspResults[["modelSummaryTable"]] <- NULL
+      jaspResults[["rSquaredTable"]]     <- NULL
     }
   }
 
@@ -1700,6 +1706,8 @@ ClassicProcess <- function(jaspResults, dataset = NULL, options) {
 
   if (doFit) {
     if (options$errorCalculationMethod == "bootstrap") {
+      # Seed per model, so a model's bootstrap does not depend on which other models were (re)fitted in this run
+      jaspBase::.setSeedJASP(options)
       fittedModel <- .procBootstrap(fittedModel, options$bootstrapSamples)
     }
   }
@@ -1724,6 +1732,12 @@ ClassicProcess <- function(jaspResults, dataset = NULL, options) {
   ))
 }
 
+# Dependencies of output that combines all models: only the options that change a fitted model,
+# so that requesting output for a single model (e.g., a plot) does not update it
+.procGetAllModelsDependencies <- function(options) {
+  return(unlist(lapply(seq_along(options[["processModels"]]), function(i) .procGetSingleModelsDependencies(as.character(i))), recursive = FALSE))
+}
+
 .procComputeResults <- function(jaspResults, dataset, options) {
   modelsContainer <- jaspResults[["modelsContainer"]]
   nModels <- length(options[["processModels"]])
@@ -1887,7 +1901,10 @@ ClassicProcess <- function(jaspResults, dataset = NULL, options) {
   tableRowIsValid <- modelNumberIsValid & resultIsValid
 
   summaryTable <- createJaspTable(title = gettext("Model summary"), rowNames = modelNames[tableRowIsValid])
-  summaryTable$dependOn(c(.procGetDependencies(), "processModels", "aicWeights", "bicWeights", "naAction", "hayesNumber"))
+  summaryTable$dependOn(
+    options = c(.procGetDependencies(), "aicWeights", "bicWeights", "naAction", "hayesNumber"),
+    nestedOptions = .procGetAllModelsDependencies(options)
+  )
   summaryTable$position <- 0
 
   if (options[["naAction"]] == "fiml" && !options[["estimator"]] %in% c("default", "ml")) {
@@ -1968,7 +1985,7 @@ ClassicProcess <- function(jaspResults, dataset = NULL, options) {
   if (any(validModel) && is.null(parEstContainer[["warning"]])) {
     warningHtml <- createJaspHtml(text = .procCausalAssumptionsMsg())
     warningHtml$position <- 0
-    warningHtml$dependOn("processModels")
+    warningHtml$dependOn(nestedOptions = .procGetAllModelsDependencies(options))
     parEstContainer[["warning"]] <- warningHtml
   }
 
@@ -2607,6 +2624,8 @@ ClassicProcess <- function(jaspResults, dataset = NULL, options) {
 
   graph <- dagitty::dagitty(paste("dag {", paste(arrows, collapse = "\n")," } "))
 
+  # Bootstrapped local tests (and cis.loess) resample the data
+  jaspBase::.setSeedJASP(options)
   localTestResult <- dagitty::localTests(
     graph,
     dataset,
@@ -2658,7 +2677,10 @@ ClassicProcess <- function(jaspResults, dataset = NULL, options) {
                         overtitle = "R\u00B2", type = "number")
   }
 
-  tabr2$dependOn(c(.procGetDependencies(), "rSquared", "processModels"))
+  tabr2$dependOn(
+    options = c(.procGetDependencies(), "rSquared"),
+    nestedOptions = .procGetAllModelsDependencies(options)
+  )
   tabr2$position <- 1
 
   jaspResults[["rSquaredTable"]] <- tabr2
@@ -2756,6 +2778,171 @@ ClassicProcess <- function(jaspResults, dataset = NULL, options) {
   procPathPlot$plotObject <- .procLavToGraph(modelContainer, type = "statistical", estimates = options[["statisticalPathPlotsParameterEstimates"]], options)
 }
 
+# Residual plots ----
+.procContainerResidualPlots <- function(jaspResults, options) {
+  if (!is.null(jaspResults[["residualPlotContainer"]]))
+    return(jaspResults[["residualPlotContainer"]])
+
+  residualPlotContainer <- createJaspContainer(gettext("Residual plots"))
+  residualPlotContainer$dependOn(.procGetDependencies())
+
+  jaspResults[["residualPlotContainer"]] <- residualPlotContainer
+
+  return(residualPlotContainer)
+}
+
+.procResidualPlots <- function(jaspResults, options, modelsContainer) {
+  showPlots <- sapply(options[["processModels"]], function(mod) mod[["residualVsFittedPlot"]] || mod[["residualQqPlot"]])
+
+  if (!any(showPlots)) return()
+
+  container <- .procContainerResidualPlots(jaspResults, options)
+
+  for (i in which(showPlots)) {
+    modelName <- options[["processModels"]][[i]][["name"]]
+
+    # Models that cannot be fitted yet (e.g., incomplete Hayes configurations) have no residuals
+    if (!.procCheckFitModel(modelsContainer[[modelName]][["graph"]]$object)) next
+
+    if (is.null(container[[modelName]])) {
+      modelContainer <- createJaspContainer(title = modelName)
+      modelContainer$dependOn(
+        nestedOptions = .procGetSingleModelsDependencies(as.character(i))
+      )
+      container[[modelName]] <- modelContainer
+    } else {
+      modelContainer <- container[[modelName]]
+    }
+
+    fittedModel <- modelsContainer[[modelName]][["fittedModel"]]$object
+    .procSetContainerError(modelContainer, fittedModel)
+
+    if (modelContainer$getError()) next
+
+    residuals <- try(.procResiduals(fittedModel))
+
+    if (jaspBase::isTryError(residuals)) {
+      modelContainer$setError(jaspBase::.extractErrorMessage(residuals))
+      next
+    }
+
+    for (j in seq_along(residuals)) {
+      .procResidualPlotsSingleVariable(modelContainer, options, residuals[[j]], names(residuals)[j], j, i, fittedModel)
+    }
+  }
+}
+
+.procResiduals <- function(fittedModel) {
+  # Residuals of the regression equation of each dependent variable (outcome and mediators):
+  # observed values minus the values predicted by the path coefficients
+  parTbl <- lavaan::parameterTable(fittedModel)
+  data   <- lavaan::lavInspect(fittedModel, "data")
+  regs   <- parTbl[parTbl$op == "~", ]
+
+  residuals <- lapply(unique(regs$lhs), function(dep) {
+    reg       <- regs[regs$lhs == dep, ]
+    preds     <- data[, reg$rhs, drop = FALSE]
+    intercept <- parTbl$est[parTbl$lhs == dep & parTbl$op == "~1"]
+
+    # Without a mean structure the intercepts are implied by the sample means
+    if (length(intercept) == 0)
+      intercept <- mean(data[, dep], na.rm = TRUE) - sum(reg$est * colMeans(preds, na.rm = TRUE))
+
+    predicted <- drop(intercept + preds %*% reg$est)
+    residual  <- data[, dep] - predicted
+    resVar    <- parTbl$est[parTbl$lhs == dep & parTbl$op == "~~" & parTbl$rhs == dep]
+    complete  <- !is.na(residual)
+
+    return(data.frame(
+      predicted    = predicted[complete],
+      residual     = residual[complete],
+      standardized = residual[complete] / sqrt(resVar)
+    ))
+  })
+  names(residuals) <- unique(regs$lhs)
+
+  return(residuals)
+}
+
+.procResidualPlotsSingleVariable <- function(container, options, residuals, varName, varIdx, modelIdx, fittedModel) {
+  varKey <- paste0("variable", varIdx)
+
+  if (is.null(container[[varKey]])) {
+    varContainer <- createJaspContainer(title = decodeColNames(varName))
+    varContainer$position <- varIdx
+    container[[varKey]] <- varContainer
+  } else {
+    varContainer <- container[[varKey]]
+  }
+
+  modelOptions <- options[["processModels"]][[modelIdx]]
+
+  if (modelOptions[["residualVsFittedPlot"]] && is.null(varContainer[["residualVsFittedPlot"]])) {
+    residualPlot <- createJaspPlot(title = gettext("Residuals vs. Predicted"), width = 530, height = 400)
+    residualPlot$dependOn(nestedOptions = list(c("processModels", as.character(modelIdx), "residualVsFittedPlot")))
+    residualPlot$position <- 1
+    varContainer[["residualVsFittedPlot"]] <- residualPlot
+
+    residualPlot$plotObject <- .procResidualVsFittedPlot(residuals)
+  }
+
+  if (modelOptions[["residualQqPlot"]] && is.null(varContainer[["residualQqPlot"]])) {
+    qqPlot <- createJaspPlot(title = gettext("Q-Q Plot Standardized Residuals"), width = 400, height = 400)
+    qqPlot$dependOn(nestedOptions = list(
+      c("processModels", as.character(modelIdx), "residualQqPlot"),
+      c("processModels", as.character(modelIdx), "qqPlotCi"),
+      c("processModels", as.character(modelIdx), "qqPlotCiLevel")
+    ))
+    qqPlot$position <- 2
+    varContainer[["residualQqPlot"]] <- qqPlot
+
+    ciLevel <- if (modelOptions[["qqPlotCi"]]) modelOptions[["qqPlotCiLevel"]] else NULL
+
+    if (inherits(fittedModel, "blavaan")) {
+      # Posterior means of the standardized residuals, with the credible interval of each as an error bar
+      stdResiduals <- try(.procBayesStandardizedResiduals(fittedModel, varName, ciLevel))
+
+      if (jaspBase::isTryError(stdResiduals)) {
+        qqPlot$setError(jaspBase::.extractErrorMessage(stdResiduals))
+      } else {
+        qqPlot$plotObject <- jaspGraphs::plotQQnorm(
+          stdResiduals[["mean"]],
+          lower         = stdResiduals[["lower"]],
+          upper         = stdResiduals[["upper"]],
+          yName         = gettext("Standardized residuals"),
+          ablineColor   = "darkred",
+          ablineOrigin  = TRUE,
+          identicalAxes = TRUE
+        )
+      }
+    } else {
+      qqPlot$plotObject <- jaspGraphs::plotQQnorm(
+        residuals[["standardized"]],
+        yName         = gettext("Standardized residuals"),
+        ablineColor   = "darkred",
+        ablineOrigin  = TRUE,
+        identicalAxes = TRUE,
+        ciLevel       = ciLevel
+      )
+    }
+  }
+}
+
+.procResidualVsFittedPlot <- function(residuals) {
+  xBreaks <- jaspGraphs::getPrettyAxisBreaks(residuals[["predicted"]])
+  yBreaks <- jaspGraphs::getPrettyAxisBreaks(c(residuals[["residual"]], 0))
+
+  plotObj <- ggplot2::ggplot(residuals, ggplot2::aes(x = predicted, y = residual)) +
+    ggplot2::geom_hline(yintercept = 0, colour = "darkred") +
+    jaspGraphs::geom_point() +
+    ggplot2::scale_x_continuous(name = gettext("Predicted Values"), breaks = xBreaks, limits = range(xBreaks)) +
+    ggplot2::scale_y_continuous(name = gettext("Residuals"), breaks = yBreaks, limits = range(yBreaks)) +
+    jaspGraphs::geom_rangeframe() +
+    jaspGraphs::themeJaspRaw()
+
+  return(plotObj)
+}
+
 .procMainGraphLayoutPosHelper <- function(nodes) {
   # This function positions nodes alternatingly above and below zero with increasing distance
   nUniqueNodes <- length(unique(nodes))

inst/help/BayesianProcess.md

@@ -64,6 +64,10 @@ This section allows users to specify multiple process models through one of two
 - Path plots: Displays path plots in the output.
 	- Conceptual: Shows a conceptual path plot. These plots contain no parameter labels or estimates and simplify moderation effects in the model structure by omitting direct effects and interactions. Instead moderators point to the path they are moderating. These plots can be shown even when the model is not complete yet (when `Input type` is `Model number`) and are easier to understand than the statistical path plots.
 	- Statistical: Shows a statistical path plot including parameter labels or estimates for all paths.
+- Residual plots: Displays residual plots for each dependent variable in the model (the outcome and each mediator), to check the assumptions of its regression equation. The predicted values are based on the posterior means of the parameters.
+	- Residuals vs. predicted: Scatterplot of the residuals against the predicted values. If the variance of the residuals is constant (homoscedasticity), the residuals scatter evenly around zero across the predicted values.
+	- Q-Q plot standardized residuals: Q-Q plot of the residuals divided by the residual standard deviation. The standardized residuals are computed for each posterior draw, and the points show their posterior means. If the residuals are normally distributed, the points lie close to the diagonal.
+		- Credible interval: Adds the central credible interval of each standardized residual as an error bar.
 
 #### Options