From 940caa7db9c43d8f6cb770754c6da965b8283f3e Mon Sep 17 00:00:00 2001
From: Johnny <j.b.vandoorn@uva.nl>
Date: Fri, 6 Nov 2026 13:06:09 +0100
Subject: [PATCH 1/7] feat: add covariates to the Hayes configuration input

Covariates enter every mediator and outcome equation (the PROCESS
default) and can be continuous or categorical. They are excluded from
the mediation/total effects tables, the conceptual path plot and Hayes
number recognition, are labelled f (into mediators) and g (into the
outcome) as in the processR package, get their own plot colour, and are
drawn below X in the statistical path plot.

Closes #63

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
---
 NEWS.md                                       |   3 +
 R/classicProcess.R                            |  82 ++++++++++++--
 inst/qml/common/InputModelNumber.qml          |  13 ++-
 .../testthat/test-bayesian-process-general.R  |  30 ++++++
 tests/testthat/test-classic-process-unit.R    | 100 ++++++++++++++++++
 .../test-classic-process-verification.R       |  93 ++++++++++++++++
 6 files changed, 306 insertions(+), 15 deletions(-)

diff --git a/NEWS.md b/NEWS.md
index 6e707bfb..72fc6340 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -16,6 +16,9 @@
 
 # jaspProcess (development version)
 
+## Added
+* Covariates box for the "Hayes configuration" input type. As in PROCESS, covariates are added to every mediator and outcome equation, so the path coefficients, effects, and R² are adjusted for them. Covariates can be continuous or categorical, and they are left out of the mediation and total effects tables and the conceptual path plot ([Issue #63](https://github.com/jasp-stats/jaspProcess/issues/63), [Issue #4388](https://github.com/jasp-stats/jasp-issues/issues/4388), [Issue #3481](https://github.com/jasp-stats/jasp-issues/issues/3481)).
+
 ## Fixed
 * 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)).
 
diff --git a/R/classicProcess.R b/R/classicProcess.R
index 3d85954b..8f1dae6b 100644
--- a/R/classicProcess.R
+++ b/R/classicProcess.R
@@ -173,14 +173,44 @@ ClassicProcess <- function(jaspResults, dataset = NULL, options) {
 
   graph <- .procProcessRelationshipsToGraph(processRelationships, modelOptions[["inputType"]] == "inputVariables")
 
-  if (modelOptions[["inputType"]] == "inputModelNumber")
+  if (modelOptions[["inputType"]] == "inputModelNumber") {
     graph <- .procModelGraphInputModelNumber(graph, modelOptions, globalDependent)
-  
+    covariates <- unlist(modelOptions[["modelNumberCovariates"]])
+    if (length(covariates) > 0)
+      graph <- .procGraphAddCovariates(graph, encodeColNames(covariates))
+  }
+
   if (!igraph::is_dag(graph)) stop(.procDagMsg())
 
   return(graph)
 }
 
+.procGraphAddCovariates <- function(graph, covariates) {
+  # Covariates (statistical controls) are regressed into every mediator and
+  # outcome equation, but not into the equation of X, as in the default of the
+  # PROCESS macro. They are stored as a graph attribute, so that later steps can
+  # keep them out of the focal effects (see `isCov` in .procGraphAddAttributes)
+  if (length(covariates) == 0) return(graph)
+
+  # Endogenous variables that are not interaction terms: mediators and outcome
+  targets <- igraph::V(graph)[!isExo & !isInt]$name
+
+  covNotInGraph <- covariates[!covariates %in% igraph::V(graph)$name]
+  graph <- igraph::add_vertices(graph, length(covNotInGraph), name = covNotInGraph)
+
+  edgeSources <- rep(covariates, each = length(targets))
+  edgeTargets <- rep(targets, times = length(covariates))
+  graph <- igraph::add_edges(graph,
+    edges = as.vector(rbind(edgeSources, edgeTargets)),
+    source = edgeSources,
+    target = edgeTargets
+  )
+
+  igraph::graph_attr(graph, "covariates") <- covariates
+
+  return(.procGraphAddAttributes(graph))
+}
+
 .procProcessRelationshipsToGraph <- function(processRelationships, inputVariables = TRUE) {
   graph <- igraph::make_empty_graph()
   
@@ -319,8 +349,10 @@ ClassicProcess <- function(jaspResults, dataset = NULL, options) {
   igraph::V(graph)$isPartOfThreewayInt <- sapply(igraph::V(graph)$name, function(v) {
     any(sapply(igraph::V(graph)$intVars, function(vars) length(vars) > 2))
   })
-  # Is treatment variable (i.e., X): When exo, not int, and not part of int
-  igraph::V(graph)$isTreat <- igraph::V(graph)$isExo & !igraph::V(graph)$isInt & !igraph::V(graph)$isPartOfInt
+  # Is covariate (statistical control, see .procGraphAddCovariates)
+  igraph::V(graph)$isCov <- igraph::V(graph)$name %in% igraph::graph_attr(graph, "covariates")
+  # Is treatment variable (i.e., X): When exo, not int, not part of int, and not a covariate
+  igraph::V(graph)$isTreat <- igraph::V(graph)$isExo & !igraph::V(graph)$isInt & !igraph::V(graph)$isPartOfInt & !igraph::V(graph)$isCov
   # Is edge moderated
   igraph::E(graph)$isMod <- FALSE
 
@@ -559,6 +591,10 @@ ClassicProcess <- function(jaspResults, dataset = NULL, options) {
     # Split terms of predictor vars
     sourcVarsSplit <- strsplit(unique(igraph::E(graph)$source), "__")
 
+    # Covariates are carried over to the new graph under their dummy-coded names
+    covariates <- igraph::graph_attr(graph, "covariates")
+    facCovariates <- character(0)
+
     # Replace dummy-coded variables in graph
     # Goes through all source variables and replaces them with dummy coded variable names
     # if necessary
@@ -575,6 +611,9 @@ ClassicProcess <- function(jaspResults, dataset = NULL, options) {
       # Concatenate source var and interaction terms in dummy coded variable
       contr <- .doCallPaste(contr, sep = "__")
       sourceName <- paste(vars, collapse = "__")
+
+      if (sourceName %in% covariates)
+        facCovariates <- c(facCovariates, contr)
       
 
       # Add nodes for source variables not in graph
@@ -604,6 +643,9 @@ ClassicProcess <- function(jaspResults, dataset = NULL, options) {
       }
     }
 
+    if (length(facCovariates) > 0)
+      igraph::graph_attr(facGraph, "covariates") <- facCovariates
+
     # Add attributes to graph
     facGraph <- .procGraphAddAttributes(facGraph)
 
@@ -632,8 +674,11 @@ ClassicProcess <- function(jaspResults, dataset = NULL, options) {
 
 .procGraphAddParNamesSingleModel <- function(graph) {
   igraph::E(graph)$parSymbol <- ifelse(
-    # If source is mediator and target is global dependent -> b
-    igraph::V(graph)[igraph::E(graph)$source]$isMed & igraph::V(graph)[igraph::E(graph)$target]$isDep, "b",
+    # If source is covariate -> f into a mediator, g into the global dependent
+    # (as in the processR package), so that covariate paths are not read as a or c paths
+    igraph::V(graph)[igraph::E(graph)$source]$isCov, ifelse(igraph::V(graph)[igraph::E(graph)$target]$isDep, "g", "f"),
+    # Else If source is mediator and target is global dependent -> b
+    ifelse(igraph::V(graph)[igraph::E(graph)$source]$isMed & igraph::V(graph)[igraph::E(graph)$target]$isDep, "b",
       # Else If source is mediator and target is not global dependent (i.e., other mediator) -> d
       ifelse(igraph::V(graph)[igraph::E(graph)$source]$isMed & !igraph::V(graph)[igraph::E(graph)$target]$isDep, "d",
         # Else If source is exogenous and target is mediator -> a
@@ -641,6 +686,7 @@ ClassicProcess <- function(jaspResults, dataset = NULL, options) {
           # Else -> c
           "c")
       )
+    )
   )
   # Get occurence of each symbol
   parSymbolTable <- table(igraph::E(graph)$parSymbol)
@@ -1717,6 +1763,7 @@ ClassicProcess <- function(jaspResults, dataset = NULL, options) {
     c("processModels", modelIdx, "modelNumberMediators"),
     c("processModels", modelIdx, "modelNumberModeratorW"),
     c("processModels", modelIdx, "modelNumberModeratorZ"),
+    c("processModels", modelIdx, "modelNumberCovariates"),
     c("processModels", modelIdx, "independentCovariances"),
     c("processModels", modelIdx, "mediatorCovariances"),
     c("processModels", modelIdx, "dependentCovariances"),
@@ -1814,6 +1861,12 @@ ClassicProcess <- function(jaspResults, dataset = NULL, options) {
 }
 
 .procRecognizeModelNumber <- function(graph) {
+  # Covariates do not define the Hayes configuration, so leave them out
+  if (any(igraph::V(graph)$isCov)) {
+    graph <- igraph::delete_vertices(graph, igraph::V(graph)[isCov])
+    graph <- igraph::delete_graph_attr(graph, "covariates")
+    graph <- .procGraphAddAttributes(graph)
+  }
   # Get global dependent variable
   globalDependent <- igraph::V(graph)[isDep]$name
   # Get mediator terms
@@ -2953,6 +3006,14 @@ ClassicProcess <- function(jaspResults, dataset = NULL, options) {
       }
     }
   }
+
+  # Place covariates in the column of X, below all other nodes
+  nCov <- sum(igraph::V(graph)$isCov)
+  if (nCov > 0) {
+    igraph::V(graph)[isCov]$posX <- 0
+    igraph::V(graph)[isCov]$posY <- min(igraph::V(graph)$posY, na.rm = TRUE) - seq_len(nCov)
+  }
+
   return(graph)
 }
 
@@ -2962,6 +3023,9 @@ ClassicProcess <- function(jaspResults, dataset = NULL, options) {
   if (!.procCheckGraph(graph)) return()
   
   if (type == "conceptual") {
+    # Covariates are statistical controls, not part of the conceptual model, so they are not shown
+    if (any(igraph::V(graph)$isCov))
+      graph <- igraph::delete_vertices(graph, igraph::V(graph)[isCov])
     graph <- .procGraphLayoutConceptual(graph)
   } else {
     graph <- .procGraphLayoutStatistical(graph)
@@ -3011,8 +3075,10 @@ ClassicProcess <- function(jaspResults, dataset = NULL, options) {
   # Create node type variable for coloring
   nodeType <- ifelse(igraph::V(graph)$isMed, "Mediator",
     ifelse(igraph::V(graph)$isInt | igraph::V(graph)$isPartOfInt, "Moderator",
-      ifelse(igraph::V(graph)$isExo, "Independent",
-        ifelse(igraph::V(graph)$isDep, "Dependent", NA)
+      ifelse(igraph::V(graph)$isCov, "Covariate",
+        ifelse(igraph::V(graph)$isExo, "Independent",
+          ifelse(igraph::V(graph)$isDep, "Dependent", NA)
+        )
       )
     )
   )
diff --git a/inst/qml/common/InputModelNumber.qml b/inst/qml/common/InputModelNumber.qml
index ee71bbd0..e64cddf6 100644
--- a/inst/qml/common/InputModelNumber.qml
+++ b/inst/qml/common/InputModelNumber.qml
@@ -54,13 +54,6 @@ Group
 			title: 				qsTr("Mediators M")
 			allowedColumns: 	["scale"]
 		}
-		// TODO
-		// AssignedVariablesList
-		// {
-		//     name: 				"modelNumberCovariates"
-		//     title: 				qsTr("Covariates")
-		//     allowedColumns: 	["scale", "ordinal"]
-		// }
 		AssignedVariablesList
 		{
 			id: modelNumberModeratorW
@@ -75,5 +68,11 @@ Group
 			title: 				qsTr("Moderator Z")
 			singleVariable: 	true
 		}
+		// Covariates enter every mediator and outcome equation (PROCESS default)
+		AssignedVariablesList
+		{
+			name: 				"modelNumberCovariates"
+			title: 				qsTr("Covariates")
+		}
 	}
 }
diff --git a/tests/testthat/test-bayesian-process-general.R b/tests/testthat/test-bayesian-process-general.R
index 5c929935..157a3e58 100644
--- a/tests/testthat/test-bayesian-process-general.R
+++ b/tests/testthat/test-bayesian-process-general.R
@@ -105,3 +105,33 @@ test_that("Incomplete Hayes configuration works", {
 	testPlot <- results[["state"]][["figures"]][[plotName]][["obj"]]
 	jaspTools::expect_equal_plots(testPlot, "statistical-path-plot-incomplete")
 })
+
+test_that("Covariates in Hayes configuration work", {
+  skip_on_ci()
+  options <- getOptionsBayesian()
+  options$hayesNumber <- TRUE
+  options$processModels <- list(getProcessModel(list()))
+  options$processModels[[1]]$inputType <- "inputModelNumber"
+  options$processModels[[1]]$modelNumber <- 4
+  options$processModels[[1]]$modelNumberIndependent <- "contGamma"
+  options$processModels[[1]]$modelNumberMediators <- list("debCollin1")
+  options$processModels[[1]]$modelNumberCovariates <- list("contcor2", "facGender")
+  set.seed(1)
+  results <- jaspTools::runAnalysis("BayesianProcess", "debug.csv", options)
+
+  # Covariates enter the mediator and outcome equations
+  coefTable <- results[["results"]][["parEstContainer"]][["collection"]][["parEstContainer_Model 1"]][["collection"]][["parEstContainer_Model 1_pathCoefficientsTable"]][["data"]]
+  paths <- sapply(coefTable, function(row) paste(row[["lhs"]], row[["rhs"]], sep = "->"))
+  expect_setequal(paths, c(
+    "contGamma->contNormal", "debCollin1->contNormal", "contcor2->contNormal", "facGenderm->contNormal",
+    "contGamma->debCollin1", "contcor2->debCollin1", "facGenderm->debCollin1"
+  ))
+
+  # Covariates are not focal: only the direct and indirect effect of X are reported
+  medTable <- results[["results"]][["parEstContainer"]][["collection"]][["parEstContainer_Model 1"]][["collection"]][["parEstContainer_Model 1_mediationEffectsTable"]][["data"]]
+  expect_length(medTable, 2)
+  expect_true(all(sapply(medTable, function(row) row[["lhs_1"]]) == "contGamma"))
+
+  summaryTable <- results[["results"]][["modelSummaryTable"]][["data"]]
+  expect_equal(summaryTable[[1]][["modelNumber"]], 4)
+})
diff --git a/tests/testthat/test-classic-process-unit.R b/tests/testthat/test-classic-process-unit.R
index ed5b5be3..68956ebe 100644
--- a/tests/testthat/test-classic-process-unit.R
+++ b/tests/testthat/test-classic-process-unit.R
@@ -835,3 +835,103 @@ test_that("Test that .procSafeFactorLevels only relabels declared factor columns
   expect_equal(levels(result[["X"]]), c("a.b", "c.d"))
   expect_identical(result[["Z"]], c("e f", "g h"))
 })
+
+# Covariates in the Hayes configuration input ------------------------------
+
+createDummyGraphModelMediationCovariates <- function() {
+  # Simple mediation (X = contGamma, M = contcor1, Y = contNormal) with two covariates
+  edgeList <- matrix(c("contGamma", "contNormal",
+                       "contGamma", "contcor1",
+                       "contcor1", "contNormal"
+  ), ncol = 2, byrow = TRUE)
+  graph <- create_graph_from_edgeList(edgeList)
+  return(jaspProcess:::.procGraphAddCovariates(graph, c("contcor2", "debCollin1")))
+}
+
+test_that("Test that .procGraphAddCovariates works", {
+  graph <- createDummyGraphModelMediationCovariates()
+  edges <- paste(igraph::E(graph)$source, igraph::E(graph)$target, sep = "->")
+
+  # Covariates enter the mediator and outcome equations, but not the equation of X
+  expect_setequal(edges, c(
+    "contGamma->contNormal", "contGamma->contcor1", "contcor1->contNormal",
+    "contcor2->contNormal", "contcor2->contcor1", "debCollin1->contNormal", "debCollin1->contcor1"
+  ))
+  expect_equal(igraph::V(graph)[isCov]$name, c("contcor2", "debCollin1"))
+  # Covariates are exogenous but they are not treatment variables
+  expect_true(all(igraph::V(graph)[isCov]$isExo))
+  expect_equal(igraph::V(graph)[isTreat]$name, "contGamma")
+  expect_equal(igraph::V(graph)[isMed]$name, "contcor1")
+})
+
+test_that("Test that .procGraphAddCovariates works - no covariates", {
+  graph <- createDummyGraphModelModeratedMediation()
+  expect_identical(jaspProcess:::.procGraphAddCovariates(graph, character(0)), graph)
+  expect_false(any(igraph::V(graph)$isCov))
+})
+
+test_that("Test that .procGraphAddCovariates works - moderated mediation", {
+  graph <- jaspProcess:::.procGraphAddCovariates(createDummyGraphModelModeratedMediation(), "debCollin1")
+
+  # No edges into moderators or interaction terms
+  expect_setequal(igraph::E(graph)[.from("debCollin1")]$target, c("contNormal", "contcor1"))
+  # Moderation structure is unchanged
+  expect_true(igraph::E(graph)[source == "contGamma" & target == "contNormal"]$isMod)
+  expect_equal(igraph::V(graph)[isTreat]$name, "contGamma")
+})
+
+test_that("Test that .procModelGraphSingleModel adds covariates for Hayes configurations", {
+  modelOptions <- list(
+    inputType = "inputModelNumber",
+    modelNumber = 4,
+    modelNumberIndependent = "contGamma",
+    modelNumberMediators = list("contcor1"),
+    modelNumberModeratorW = "",
+    modelNumberModeratorZ = "",
+    modelNumberCovariates = list("contcor2")
+  )
+  graph <- jaspProcess:::.procModelGraphSingleModel(modelOptions, globalDependent = "contNormal")
+
+  expect_setequal(igraph::E(graph)[.from("contcor2")]$target, c("contNormal", "contcor1"))
+  expect_equal(igraph::V(graph)[isCov]$name, "contcor2")
+
+  # Analyses saved before the option existed have no covariates
+  modelOptions[["modelNumberCovariates"]] <- NULL
+  graph <- jaspProcess:::.procModelGraphSingleModel(modelOptions, globalDependent = "contNormal")
+
+  expect_false(any(igraph::V(graph)$isCov))
+  expect_equal(igraph::ecount(graph), 3)
+})
+
+test_that("Test that .procGraphAddParNamesSingleModel works - covariates", {
+  graph <- jaspProcess:::.procGraphAddParNamesSingleModel(createDummyGraphModelMediationCovariates())
+  # Covariate paths: f into the mediator, g into the outcome
+  expect_equal(igraph::E(graph)$parSymbol, c("c", "a", "b", "g", "f", "g", "f"))
+  expect_equal(igraph::E(graph)$parName, c("c1", "a1", "b1", "g1", "f1", "g2", "f2"))
+})
+
+test_that("Test that .procRecognizeModelNumber ignores covariates", {
+  edgeList <- matrix(c("contGamma", "contNormal",
+                       "contGamma", "contcor1",
+                       "contcor1", "contNormal"
+  ), ncol = 2, byrow = TRUE)
+  graphMediation <- create_graph_from_edgeList(edgeList)
+
+  expect_equal(
+    jaspProcess:::.procRecognizeModelNumber(createDummyGraphModelMediationCovariates()),
+    jaspProcess:::.procRecognizeModelNumber(graphMediation)
+  )
+  expect_true(4 %in% jaspProcess:::.procRecognizeModelNumber(createDummyGraphModelMediationCovariates()))
+
+  graph_5_14 <- jaspProcess:::.procGraphAddCovariates(createDummyGraphModelModeratedMediation(), "debCollin1")
+  expect_equal(jaspProcess:::.procRecognizeModelNumber(graph_5_14), c(5, 14))
+})
+
+test_that("Test that .procGraphLayoutStatistical places covariates below X", {
+  graph <- jaspProcess:::.procGraphLayoutStatistical(createDummyGraphModelMediationCovariates())
+  covNodes <- igraph::V(graph)[isCov]
+
+  expect_equal(covNodes$posX, c(0, 0))
+  expect_true(all(covNodes$posY < min(igraph::V(graph)[!isCov]$posY)))
+  expect_false(anyNA(igraph::V(graph)$posY))
+})
diff --git a/tests/testthat/test-classic-process-verification.R b/tests/testthat/test-classic-process-verification.R
index ce6616b6..1570f5f5 100644
--- a/tests/testthat/test-classic-process-verification.R
+++ b/tests/testthat/test-classic-process-verification.R
@@ -82,3 +82,96 @@ test_that("Test that moderated mediation effect matches validated results", {
 	expect_equal(medEstimates, rev(refDf$ind), tolerance = 1e-7)
 	expect_equal(medEstimatesStd, rev(refDf$std), tolerance = 1e-7)
 })
+
+# Covariates in the Hayes configuration input enter every mediator and outcome
+# equation (the PROCESS default), so the path coefficients must equal the OLS
+# estimates of these equations (jasp-stats/jaspProcess#63)
+
+getHayesModelCovariates <- function(modelNumber, independent, mediators = list(), modW = "", covariates = list()) {
+  model <- getProcessModel(list())
+  model$inputType <- "inputModelNumber"
+  model$modelNumber <- modelNumber
+  model$modelNumberIndependent <- independent
+  model$modelNumberMediators <- mediators
+  model$modelNumberModeratorW <- modW
+  model$modelNumberModeratorZ <- ""
+  model$modelNumberCovariates <- covariates
+  return(model)
+}
+
+getPathCoefficients <- function(results, dependent) {
+  coefTable <- results[["results"]][["parEstContainer"]][["collection"]][["parEstContainer_Model 1"]][["collection"]][["parEstContainer_Model 1_pathCoefficientsTable"]][["data"]]
+  coefTable <- Filter(function(row) row[["rhs"]] == dependent, coefTable)
+  return(setNames(sapply(coefTable, function(row) row[["est"]]), sapply(coefTable, function(row) row[["lhs"]])))
+}
+
+test_that("Test that moderation with covariates in Hayes configuration matches linear regression", {
+  options <- getOptionsClassical()
+  options$dependent <- "GovernmentAction"
+  options$covariates <- list("NegativeEmotion", "Age", "PositiveEmotion", "Ideology", "Sex")
+  options$factors <- list()
+  options$naAction <- "listwise"
+  options$meanCenteredModeration <- FALSE
+  options$standardizedModelEstimates <- FALSE
+  options$processModels <- list(getHayesModelCovariates(1, "NegativeEmotion", modW = "Age",
+                                                        covariates = list("PositiveEmotion", "Ideology", "Sex")))
+  set.seed(1)
+  results <- jaspTools::runAnalysis("ClassicProcess", test_path("global_warming.csv"), options)
+
+  data <- read.csv(test_path("global_warming.csv"))
+  fit <- lm(GovernmentAction ~ NegativeEmotion * Age + PositiveEmotion + Ideology + Sex, data = data)
+  refCoefs <- coef(fit)[-1]
+
+  coefs <- getPathCoefficients(results, "GovernmentAction")
+  expect_setequal(names(coefs), names(refCoefs))
+  expect_equal(coefs[names(refCoefs)], refCoefs, tolerance = 1e-6)
+
+  # Covariates are not focal: only conditional effects of X are reported
+  medTable <- results[["results"]][["parEstContainer"]][["collection"]][["parEstContainer_Model 1"]][["collection"]][["parEstContainer_Model 1_mediationEffectsTable"]][["data"]]
+  expect_length(medTable, 3)
+  expect_true(all(sapply(medTable, function(row) row[["lhs_1"]]) == "NegativeEmotion"))
+  probes <- quantile(data$Age, c(0.16, 0.5, 0.84))
+  expect_equal(sapply(medTable, function(row) row[["est"]]),
+               unname(refCoefs["NegativeEmotion"] + refCoefs["NegativeEmotion:Age"] * probes), tolerance = 1e-6)
+
+  # Covariates do not change the recognized Hayes configuration
+  summaryTable <- results[["results"]][["modelSummaryTable"]][["data"]]
+  expect_equal(summaryTable[[1]][["modelNumber"]], 1)
+
+  rSquaredTable <- results[["results"]][["rSquaredTable"]][["data"]]
+  expect_equal(rSquaredTable[[1]][["rsq_1"]], summary(fit)$r.squared, tolerance = 1e-6)
+})
+
+test_that("Test that mediation with continuous and categorical covariates in Hayes configuration matches linear regression", {
+  options <- getOptionsClassical()
+  options$naAction <- "listwise"
+  options$meanCenteredModeration <- FALSE
+  options$standardizedModelEstimates <- FALSE
+  options$processModels <- list(getHayesModelCovariates(4, "contGamma", mediators = list("debCollin1"),
+                                                        covariates = list("contcor2", "facGender")))
+  set.seed(1)
+  results <- jaspTools::runAnalysis("ClassicProcess", "debug.csv", options)
+
+  data <- read.csv(system.file("jaspData", "debug.csv", package = "jaspTools"))
+  data <- data[complete.cases(data[, c("contNormal", "contGamma", "contcor1", "contcor2", "debCollin1", "facGender", "facExperim")]), ]
+  data$facGender <- factor(data$facGender)
+  refCoefsM <- coef(lm(debCollin1 ~ contGamma + contcor2 + facGender, data = data))[-1]
+  refCoefsY <- coef(lm(contNormal ~ contGamma + debCollin1 + contcor2 + facGender, data = data))[-1]
+
+  coefsM <- getPathCoefficients(results, "debCollin1")
+  coefsY <- getPathCoefficients(results, "contNormal")
+  expect_setequal(names(coefsM), names(refCoefsM))
+  expect_setequal(names(coefsY), names(refCoefsY))
+  expect_equal(coefsM[names(refCoefsM)], refCoefsM, tolerance = 1e-6)
+  expect_equal(coefsY[names(refCoefsY)], refCoefsY, tolerance = 1e-6)
+
+  # Covariates are not focal: only the direct and indirect effect of X are reported
+  medTable <- results[["results"]][["parEstContainer"]][["collection"]][["parEstContainer_Model 1"]][["collection"]][["parEstContainer_Model 1_mediationEffectsTable"]][["data"]]
+  expect_length(medTable, 2)
+  expect_true(all(sapply(medTable, function(row) row[["lhs_1"]]) == "contGamma"))
+  expect_equal(sapply(medTable, function(row) row[["est"]]),
+               unname(c(refCoefsY["contGamma"], refCoefsM["contGamma"] * refCoefsY["debCollin1"])), tolerance = 1e-6)
+
+  summaryTable <- results[["results"]][["modelSummaryTable"]][["data"]]
+  expect_equal(summaryTable[[1]][["modelNumber"]], 4)
+})

From a77464eb803ff329d406e3d04a72a2d000924858 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Franti=C5=A1ek=20Barto=C5=A1?= <f.bartos96@gmail.com>
Date: Sat, 7 Nov 2026 13:27:28 +0100
Subject: [PATCH 2/7] fix: center and probe moderators on each model's own
 complete cases

With listwise deletion, .procMeanCenter and .procModProbesSingleModel took
complete.cases() over the whole dataset, i.e. over every variable in the
variable boxes and every dummy/interaction column of every model. lavaan,
however, deletes listwise only on the variables of the fitted model. When any
other listed variable had missing values (for example a covariate used only
in another model), a model with complete data of its own was centered on the
smaller sample, while the excluded rows entered the fit uncentered. The
moderator probes used the same restricted sample.

Example (debug.csv, Hayes model 1, X = contGamma, W = contcor1, next to a
second model with covariate debMiss30, 30 NAs): the X coefficient of the first
model was -0.0718 instead of -0.0285 (lm on the centered n = 100 sample).

The new helper .procModelCompleteCases restricts the complete cases to the
model's own variables, and is used for both centering and probing. The new
verification test compares both models with lm() on their own centered samples.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
---
 R/classicProcess.R                            | 31 +++++++-----
 .../test-classic-process-verification.R       | 47 ++++++++++++++++++-
 2 files changed, 64 insertions(+), 14 deletions(-)

diff --git a/R/classicProcess.R b/R/classicProcess.R
index 8f1dae6b..e67d64e2 100644
--- a/R/classicProcess.R
+++ b/R/classicProcess.R
@@ -751,17 +751,27 @@ ClassicProcess <- function(jaspResults, dataset = NULL, options) {
   }
 }
 
+.procModelCompleteCases <- function(graph, dataset, options) {
+  # Rows that lavaan uses for this model. With listwise deletion, these are the
+  # cases that are complete on the model's own (dummy-coded) variables only.
+  # Missing values in variables of other models, or in variables that are only
+  # listed in the variable boxes, must not restrict this model's sample:
+  # otherwise its moderators are centered and probed on a smaller sample than
+  # the one lavaan fits, and the excluded rows enter the fit uncentered.
+  if (is.null(options$naAction) || options$naAction != "listwise") {
+    return(!logical(nrow(dataset)))
+  }
+  modelVars <- intersect(igraph::V(graph)$name, names(dataset))
+  return(complete.cases(dataset[, modelVars, drop = FALSE]))
+}
+
 .procMeanCenter <- function(graph, dataset, options) {
   # Center interacting variables
   if (options$meanCenteredModeration) {
     # Only use complete cases for standardization with listwise deletion
     # Otherwise data entered into lavaan is not actually standardized
-    if (!is.null(options$naAction) && options$naAction == "listwise") {
-      isComplete <- complete.cases(dataset)
-    } else {
-      isComplete <- !logical(nrow(dataset))
-    }
-    
+    isComplete <- .procModelCompleteCases(graph, dataset, options)
+
     for (v in unique(unlist(igraph::V(graph)[intLength > 1]$intVars))) {
       # Only center continuous variables
       if (v %in% options[["covariates"]]) {
@@ -822,6 +832,9 @@ ClassicProcess <- function(jaspResults, dataset = NULL, options) {
   # Mean-center data set before quantile computation
   dataset <- .procMeanCenter(container[["graph"]]$object, dataset, options)
 
+  # If listwise missing value deletion, only use the model's complete cases for quantiles
+  isComplete <- .procModelCompleteCases(graph, dataset, options)
+
   modProbes <- lapply(modVars, function(nms) {
     # Is moderator factor
     matchFac <- sapply(options[["factors"]], grepl, x = nms)
@@ -834,12 +847,6 @@ ClassicProcess <- function(jaspResults, dataset = NULL, options) {
       # Return matrix with dummy coding for each factor
       return(conMat[, colIdx])
     }
-    # If listwise missing value deletion, only use complete cases for quantiles
-    if (options$naAction == "listwise") {
-      isComplete <- complete.cases(dataset)
-    } else {
-      isComplete <- !logical(nrow(dataset))
-    }
 
     vComplete <- dataset[[nms]][isComplete]
 
diff --git a/tests/testthat/test-classic-process-verification.R b/tests/testthat/test-classic-process-verification.R
index 1570f5f5..55e4dd75 100644
--- a/tests/testthat/test-classic-process-verification.R
+++ b/tests/testthat/test-classic-process-verification.R
@@ -99,8 +99,8 @@ getHayesModelCovariates <- function(modelNumber, independent, mediators = list()
   return(model)
 }
 
-getPathCoefficients <- function(results, dependent) {
-  coefTable <- results[["results"]][["parEstContainer"]][["collection"]][["parEstContainer_Model 1"]][["collection"]][["parEstContainer_Model 1_pathCoefficientsTable"]][["data"]]
+getPathCoefficients <- function(results, dependent, model = "Model 1") {
+  coefTable <- results[["results"]][["parEstContainer"]][["collection"]][[paste0("parEstContainer_", model)]][["collection"]][[paste0("parEstContainer_", model, "_pathCoefficientsTable")]][["data"]]
   coefTable <- Filter(function(row) row[["rhs"]] == dependent, coefTable)
   return(setNames(sapply(coefTable, function(row) row[["est"]]), sapply(coefTable, function(row) row[["lhs"]])))
 }
@@ -175,3 +175,46 @@ test_that("Test that mediation with continuous and categorical covariates in Hay
   summaryTable <- results[["results"]][["modelSummaryTable"]][["data"]]
   expect_equal(summaryTable[[1]][["modelNumber"]], 4)
 })
+
+test_that("Test that moderators are centered and probed on each model's own complete cases", {
+  # Model 2 adds a covariate with 30 missing values. With listwise deletion, these
+  # must not change the sample on which Model 1 (complete data) is centered and probed.
+  options <- getOptionsClassical()
+  options$covariates <- c(options$covariates, list("debMiss30"))
+  options$naAction <- "listwise"
+  options$meanCenteredModeration <- TRUE
+  options$standardizedModelEstimates <- FALSE
+  model2 <- getHayesModelCovariates(1, "contGamma", modW = "contcor1", covariates = list("debMiss30"))
+  model2$name <- "Model 2"
+  options$processModels <- list(getHayesModelCovariates(1, "contGamma", modW = "contcor1"), model2)
+  set.seed(1)
+  results <- jaspTools::runAnalysis("ClassicProcess", "debug.csv", options)
+
+  data <- read.csv(system.file("jaspData", "debug.csv", package = "jaspTools"))
+  centerXW <- function(d) transform(d, contGamma = contGamma - mean(contGamma), contcor1 = contcor1 - mean(contcor1))
+  data1 <- centerXW(data)
+  data2 <- centerXW(data[!is.na(data$debMiss30), ])
+  refCoefs1 <- coef(lm(contNormal ~ contGamma * contcor1, data = data1))[-1]
+  refCoefs2 <- coef(lm(contNormal ~ contGamma * contcor1 + debMiss30, data = data2))[-1]
+
+  coefs1 <- getPathCoefficients(results, "contNormal", "Model 1")
+  coefs2 <- getPathCoefficients(results, "contNormal", "Model 2")
+  expect_setequal(names(coefs1), names(refCoefs1))
+  expect_setequal(names(coefs2), names(refCoefs2))
+  expect_equal(coefs1[names(refCoefs1)], refCoefs1, tolerance = 1e-6)
+  expect_equal(coefs2[names(refCoefs2)], refCoefs2, tolerance = 1e-6)
+
+  summaryTable <- results[["results"]][["modelSummaryTable"]][["data"]]
+  expect_equal(sapply(summaryTable, function(row) row[["N"]]), c(100, 70))
+
+  # Conditional effects of X at the percentiles of the centered moderator in each model's sample
+  for (i in 1:2) {
+    model <- paste("Model", i)
+    refCoefs <- list(refCoefs1, refCoefs2)[[i]]
+    # Probe values enter the lavaan syntax with 3 significant digits
+    probes <- signif(quantile(list(data1, data2)[[i]]$contcor1, c(0.16, 0.5, 0.84)), 3)
+    medTable <- results[["results"]][["parEstContainer"]][["collection"]][[paste0("parEstContainer_", model)]][["collection"]][[paste0("parEstContainer_", model, "_mediationEffectsTable")]][["data"]]
+    expect_equal(sapply(medTable, function(row) row[["est"]]),
+                 unname(refCoefs["contGamma"] + refCoefs["contGamma:contcor1"] * probes), tolerance = 1e-6)
+  }
+})

From a41883497918ca09a608bb41793a55006af58d49 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Franti=C5=A1ek=20Barto=C5=A1?= <f.bartos96@gmail.com>
Date: Sat, 7 Nov 2026 13:39:20 +0100
Subject: [PATCH 3/7] fix: stop early for factor levels that clash with the
 interaction separator

"__" separates the terms of interaction variables (<X>__<W>), and factor
dummy variables are named <factor><level>. A factor level that contains "__"
(e.g. a covariate level "treat__A") or ends with "_" (e.g. "treat_" for a
factor X in an interaction, giving "<X>treat___<W>") is therefore split at
the wrong place when the interaction and centering variables are built. The
analysis then crashed with an unrelated error ("object ... not found" or
"subscript out of bounds").

.procErrorHandling now checks the factor levels, after the same name-safe
relabelling that .procSafeFactorLevels applies (shared helper
.procSafeLevelLabels), and stops with a message that names the levels to
rename.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
---
 R/classicProcess.R                         | 28 +++++++++++++++++++++-
 tests/testthat/test-classic-process-unit.R | 15 ++++++++++++
 2 files changed, 42 insertions(+), 1 deletion(-)

diff --git a/R/classicProcess.R b/R/classicProcess.R
index e67d64e2..cc5afa62 100644
--- a/R/classicProcess.R
+++ b/R/classicProcess.R
@@ -525,12 +525,36 @@ ClassicProcess <- function(jaspResults, dataset = NULL, options) {
   for (f in unlist(options[["factors"]])) {
     if (is.null(dataset[[f]])) next
     dataset[[f]] <- as.factor(dataset[[f]])
-    safeLevels <- gsub("[^A-Za-z0-9._]+", ".", levels(dataset[[f]]))
+    safeLevels <- .procSafeLevelLabels(levels(dataset[[f]]))
     levels(dataset[[f]]) <- make.unique(safeLevels, sep = "_")
   }
   return(dataset)
 }
 
+.procSafeLevelLabels <- function(levels) {
+  # Replace runs of characters that are invalid within a name (see .procSafeFactorLevels)
+  return(gsub("[^A-Za-z0-9._]+", ".", levels))
+}
+
+.procCheckFactorLevels <- function(dataset, factors) {
+  # "__" separates the terms of interaction variables (e.g. <X>__<W>), and
+  # factor dummy variables are named <factor><level>. A level label that
+  # contains "__" or ends with "_" (also after .procSafeLevelLabels) is
+  # therefore split at the wrong place when the interaction and centering
+  # variables are built, which crashes the analysis. Stop early with a clear
+  # message instead.
+  for (f in factors) {
+    levs <- levels(as.factor(dataset[[f]]))
+    badLevels <- levs[grepl("__|_$", .procSafeLevelLabels(levs))]
+    if (length(badLevels) > 0) {
+      .quitAnalysis(gettextf(
+        "Levels of the categorical variable %1$s must not contain two consecutive underscores (__) or end with an underscore (_), because __ separates the variables in interaction terms. Please rename these levels: %2$s.",
+        f, paste(badLevels, collapse = ", ")
+      ))
+    }
+  }
+}
+
 .procAddFactorDummyIntVars <- function(jaspResults, dataset, options) {
   modelsContainer <- jaspResults[["modelsContainer"]]
 
@@ -749,6 +773,8 @@ ClassicProcess <- function(jaspResults, dataset = NULL, options) {
       }
     }
   }
+
+  .procCheckFactorLevels(dataset, factors)
 }
 
 .procModelCompleteCases <- function(graph, dataset, options) {
diff --git a/tests/testthat/test-classic-process-unit.R b/tests/testthat/test-classic-process-unit.R
index 68956ebe..a8420cf0 100644
--- a/tests/testthat/test-classic-process-unit.R
+++ b/tests/testthat/test-classic-process-unit.R
@@ -836,6 +836,21 @@ test_that("Test that .procSafeFactorLevels only relabels declared factor columns
   expect_identical(result[["Z"]], c("e f", "g h"))
 })
 
+test_that("Test that .procCheckFactorLevels rejects levels that clash with the interaction separator", {
+  dataset <- data.frame(
+    ok    = c("ctrl", "treat_A", "_treat", "a b"),
+    dbl   = c("ctrl", "treat__A", "treat__B", "ctrl"),
+    trail = c("ctrl", "treat_", "ctrl", "ctrl"),
+    space = c("ctrl", "treat _", "ctrl", "ctrl"),
+    stringsAsFactors = FALSE
+  )
+  expect_silent(jaspProcess:::.procCheckFactorLevels(dataset, "ok"))
+  expect_error(jaspProcess:::.procCheckFactorLevels(dataset, c("ok", "dbl")), "treat__A, treat__B", class = "validationError")
+  expect_error(jaspProcess:::.procCheckFactorLevels(dataset, "trail"), "treat_\\.", class = "validationError")
+  # Also after the name-safe relabelling ("treat _" -> "treat._" ends with "_")
+  expect_error(jaspProcess:::.procCheckFactorLevels(dataset, "space"), "treat _", class = "validationError")
+})
+
 # Covariates in the Hayes configuration input ------------------------------
 
 createDummyGraphModelMediationCovariates <- function() {

From ea9bd0a78c5f584702e01ae7787318a2dfbd416e Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Franti=C5=A1ek=20Barto=C5=A1?= <f.bartos96@gmail.com>
Date: Sat, 7 Nov 2026 13:39:27 +0100
Subject: [PATCH 4/7] fix: keep node colours in path plots independent of
 covariates

The node fill colours come from JASPcolors(asFunction = TRUE), which
interpolates the palette, so every colour depends on the number of node types
in a plot. The statistical plot has a "Covariate" node type and the conceptual
plot (which leaves covariates out) does not, so the same variable got
different colours in the two plots of one model: in Hayes model 7, the
moderator was grey in the conceptual plot while grey meant "Covariate" in the
statistical plot. Adding a covariate also recoloured the mediator, moderator
and dependent variable in the statistical plot.

Covariates now get a fixed light grey outside the palette, and the palette is
interpolated over the other node types only. Plots without covariates are
unchanged, so the existing plot snapshots still match. The new test checks
that adding a covariate leaves the conceptual plot and the other nodes of the
statistical plot unchanged.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
---
 R/classicProcess.R                            | 11 +++++-
 ...test-classic-process-integration-general.R | 39 +++++++++++++++++++
 2 files changed, 48 insertions(+), 2 deletions(-)

diff --git a/R/classicProcess.R b/R/classicProcess.R
index cc5afa62..4a286524 100644
--- a/R/classicProcess.R
+++ b/R/classicProcess.R
@@ -3121,8 +3121,15 @@ ClassicProcess <- function(jaspResults, dataset = NULL, options) {
 
   nodeTypeLevels <- unique(na.omit(nodeType))
   if (options[["useColorPalette"]]) {
-    typeCols <- colorFun(length(nodeTypeLevels))
-    names(typeCols) <- nodeTypeLevels
+    # colorFun interpolates the palette, so each colour depends on the number of
+    # node types. Covariates therefore get a fixed light grey outside the palette:
+    # adding covariates, or dropping them in the conceptual plot, must not
+    # recolour the other node types
+    paletteLevels <- nodeTypeLevels[nodeTypeLevels != "Covariate"]
+    typeCols <- c(colorFun(length(paletteLevels)), "#D9D9D9")
+    names(typeCols) <- c(paletteLevels, "Covariate")
+    # Same order as nodeTypeLevels, which labels the legend
+    typeCols <- typeCols[nodeTypeLevels]
     igraph::V(graph)$nodeFill <- ifelse(is.na(nodeType), "transparent", typeCols[nodeType])
   } else {
     igraph::V(graph)$nodeFill <- "transparent"
diff --git a/tests/testthat/test-classic-process-integration-general.R b/tests/testthat/test-classic-process-integration-general.R
index 407c6a1a..23d2580f 100644
--- a/tests/testthat/test-classic-process-integration-general.R
+++ b/tests/testthat/test-classic-process-integration-general.R
@@ -1593,3 +1593,42 @@ test_that("Moderated mediation index works - three-level factor independent", {
 			 "contcor1", "contNormal", "<unicode>", "<unicode>", 0.862779808586041,
 			 0.0556237127253498, -0.172836674711031))
 })
+
+test_that("Covariates do not change the node colours of the other variables in path plots", {
+  # The palette is interpolated, so the colour of each node type depends on the
+  # number of node types; covariates get a fixed colour outside the palette
+  nodeFills <- function(results, plot) {
+    plotName <- results[["results"]][["pathPlotContainer"]][["collection"]][["pathPlotContainer_Model 1"]][["collection"]][[paste0("pathPlotContainer_Model 1_", plot)]][["data"]]
+    built <- ggplot2::ggplot_build(results[["state"]][["figures"]][[plotName]][["obj"]])
+    # Layer 1 draws the node tiles and layer 3 the node labels, both in node order
+    fills <- setNames(built$data[[1]]$fill, built$data[[3]]$label)
+    return(fills[names(fills) != ""]) # drop invisible helper nodes
+  }
+  runModel7 <- function(covariates) {
+    options <- getOptionsClassical()
+    options$useColorPalette <- TRUE
+    options$pathPlotsLabelLength <- 20 # full labels, so that nodes can be matched across plots
+    options$processModels <- list(getProcessModel(list()))
+    options$processModels[[1]]$inputType <- "inputModelNumber"
+    options$processModels[[1]]$modelNumber <- 7
+    options$processModels[[1]]$modelNumberIndependent <- "contGamma"
+    options$processModels[[1]]$modelNumberMediators <- list("debCollin1")
+    options$processModels[[1]]$modelNumberModeratorW <- "contcor1"
+    options$processModels[[1]]$modelNumberCovariates <- covariates
+    set.seed(1)
+    return(jaspTools::runAnalysis("ClassicProcess", "debug.csv", options))
+  }
+  resultsNoCov <- runModel7(list())
+  resultsCov   <- runModel7(list("contcor2"))
+
+  # The conceptual plot leaves covariates out
+  conceptFills <- nodeFills(resultsCov, "conceptPathPlot")
+  expect_identical(conceptFills, nodeFills(resultsNoCov, "conceptPathPlot"))
+
+  # In the statistical plot, only the covariate node is added, in its own colour
+  statFillsNoCov <- nodeFills(resultsNoCov, "statPathPlot")
+  statFillsCov   <- nodeFills(resultsCov, "statPathPlot")
+  expect_identical(statFillsCov[names(statFillsNoCov)], statFillsNoCov)
+  expect_identical(unname(statFillsCov["contcor2"]), "#D9D9D9")
+  expect_false("#D9D9D9" %in% conceptFills)
+})

From 217f7cb808824a6f78f5763afb574ab123fa4917 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Franti=C5=A1ek=20Barto=C5=A1?= <f.bartos96@gmail.com>
Date: Sat, 7 Nov 2026 13:46:25 +0100
Subject: [PATCH 5/7] fix: fit the Bayesian analysis on mean-centered data

The Bayesian analysis has the "Mean-centered moderation" option (on by
default), and the moderator probes are computed by the shared
.procModProbes on mean-centered data, but .procBayesResultsFitModel passed
the raw data to blavaan::bsem(). The conditional effects were therefore
evaluated at moderator values shifted by the moderator mean. In the new test
(global_warming.csv, Hayes model 1, W = Age), all three conditional effects
were off by b_XW * mean(Age) = 0.35 (posterior SD about 0.04).

The data are now centered with .procMeanCenter before fitting, as in the
classic analysis. The new test compares the posterior means with the OLS
conditional effects at the raw percentiles of Age (skip_on_ci, like the other
Bayesian tests).

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
---
 R/bayesianProcess.R                           |  6 ++++
 .../testthat/test-bayesian-process-general.R  | 31 +++++++++++++++++++
 2 files changed, 37 insertions(+)

diff --git a/R/bayesianProcess.R b/R/bayesianProcess.R
index a4633ab3..6582c986 100644
--- a/R/bayesianProcess.R
+++ b/R/bayesianProcess.R
@@ -135,6 +135,12 @@ BayesianProcess <- function(jaspResults, dataset = NULL, options) {
       do.fit          = doFit
     ))
   } else {
+    # Mean-center moderators as in the classic analysis: the moderator probes
+    # (.procModProbes) are computed on the centered data, so the conditional
+    # effects are only evaluated at the requested moderator values if the
+    # model is fitted on the centered data too
+    dataset <- .procMeanCenter(container[["graph"]]$object, dataset, options)
+
     # Necessary for JASP to find function blavaan
     blavaan <- blavaan::blavaan
 
diff --git a/tests/testthat/test-bayesian-process-general.R b/tests/testthat/test-bayesian-process-general.R
index 157a3e58..ebc21678 100644
--- a/tests/testthat/test-bayesian-process-general.R
+++ b/tests/testthat/test-bayesian-process-general.R
@@ -135,3 +135,34 @@ test_that("Covariates in Hayes configuration work", {
   summaryTable <- results[["results"]][["modelSummaryTable"]][["data"]]
   expect_equal(summaryTable[[1]][["modelNumber"]], 4)
 })
+
+test_that("Mean-centered moderation evaluates conditional effects at the requested moderator values", {
+  skip_on_ci()
+  # The moderator probes are computed on the mean-centered moderator, so the model
+  # must be fitted on centered data too. Otherwise the conditional effects are
+  # evaluated at moderator values shifted by the moderator mean (here by
+  # b_XW * mean(Age), about 0.35). With weak priors and n = 815, the posterior
+  # means are close to the OLS conditional effects at the raw percentiles.
+  options <- getOptionsBayesian()
+  options$dependent <- "GovernmentAction"
+  options$covariates <- list("NegativeEmotion", "Age")
+  options$factors <- list()
+  options$meanCenteredModeration <- TRUE
+  options$mcmcBurnin <- 500
+  options$mcmcSamples <- 1000
+  options$processModels <- list(getProcessModel(list(list(processDependent = "GovernmentAction",
+          processIndependent = "NegativeEmotion", processType = "moderators",
+          processVariable = "Age"))))
+  set.seed(1)
+  results <- jaspTools::runAnalysis("BayesianProcess", test_path("global_warming.csv"), options)
+
+  data <- read.csv(test_path("global_warming.csv"))
+  refCoefs <- coef(lm(GovernmentAction ~ NegativeEmotion * Age, data = data))
+  probes <- quantile(data$Age, c(0.16, 0.5, 0.84))
+  refEffects <- unname(refCoefs["NegativeEmotion"] + refCoefs["NegativeEmotion:Age"] * probes)
+
+  medTable <- results[["results"]][["parEstContainer"]][["collection"]][["parEstContainer_Model 1"]][["collection"]][["parEstContainer_Model 1_mediationEffectsTable"]][["data"]]
+  effects <- sapply(medTable, function(row) row[["mean"]])
+  expect_length(effects, 3)
+  expect_lt(max(abs(effects - refEffects)), 0.03)
+})

From ccfbc04792b28a575baade252341230099e7155a Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Franti=C5=A1ek=20Barto=C5=A1?= <f.bartos96@gmail.com>
Date: Sat, 7 Nov 2026 13:46:32 +0100
Subject: [PATCH 6/7] docs: describe the Covariates box in the help pages, add
 NEWS entries

The help pages of both analyses listed the Hayes configuration boxes X, M, W
and Z but not the new Covariates box. Describe what covariates do (as in
PROCESS: added to every mediator and outcome equation, paths f and g), what
they are left out of, and that a covariate cannot also have another role.
Add NEWS entries for the centering, Bayesian centering and factor-level fixes
in this branch.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
---
 NEWS.md                      | 3 +++
 inst/help/BayesianProcess.md | 1 +
 inst/help/ClassicProcess.md  | 1 +
 3 files changed, 5 insertions(+)

diff --git a/NEWS.md b/NEWS.md
index 72fc6340..13743910 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -21,6 +21,9 @@
 
 ## Fixed
 * 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)).
+* With listwise deletion, mean-centering of moderators and the moderator probe values now use each model's own complete cases. Previously, missing values in any other variable in the variable boxes (for example a covariate used only in another model) restricted the centering sample and changed the estimates of models with complete data.
+* The Bayesian analysis now fits the mean-centered data when "Mean-centered moderation" is selected (the default). Previously only the moderator probe values were centered, so conditional effects were evaluated at moderator values shifted by the moderator mean.
+* Factor levels that contain two consecutive underscores or end with an underscore now stop the analysis with a message asking to rename them, instead of causing an unrelated error later.
 
 ---
 
diff --git a/inst/help/BayesianProcess.md b/inst/help/BayesianProcess.md
index 1905c98c..ab22a267 100644
--- a/inst/help/BayesianProcess.md
+++ b/inst/help/BayesianProcess.md
@@ -58,6 +58,7 @@ This section allows users to specify multiple process models through one of two
 		- Mediators M: Mediator variables.
 		- Moderator W: First moderator variable.
 		- Moderator Z: Second moderator variable.
+		- Covariates: Control variables (continuous or categorical). As in the PROCESS macro, every covariate is added as a predictor to each mediator and outcome equation (paths labelled f into mediators and g into the outcome), so all path coefficients and effects are adjusted for the covariates. Covariates do not appear in the mediation effects tables or the conceptual path plot, and they do not change the Hayes configuration number. A covariate cannot also be the independent variable, a mediator, a moderator, or the dependent variable.
 
 - Residual Covariances: Adds covariances between the residuals of variables indicating that they have a common cause not specified in the model. Also adds variances for independent (exogenous) variables. This adds degrees of freedom corresponding to the number of covariances and variances added.
 - Parameter Estimates: Determines which parameter estimates are shown in the output. Requesting intercepts adds the mean structure to the model, which adds parameters and will therefore affect some fit indices (such as the information criteria AIC and BIC).
diff --git a/inst/help/ClassicProcess.md b/inst/help/ClassicProcess.md
index df24c125..ea835121 100644
--- a/inst/help/ClassicProcess.md
+++ b/inst/help/ClassicProcess.md
@@ -58,6 +58,7 @@ This section allows users to specify multiple process models through one of two
 		- Mediators M: Mediator variables.
 		- Moderator W: First moderator variable.
 		- Moderator Z: Second moderator variable.
+		- Covariates: Control variables (continuous or categorical). As in the PROCESS macro, every covariate is added as a predictor to each mediator and outcome equation (paths labelled f into mediators and g into the outcome), so all path coefficients and effects are adjusted for the covariates. Covariates do not appear in the mediation effects tables or the conceptual path plot, and they do not change the Hayes configuration number. A covariate cannot also be the independent variable, a mediator, a moderator, or the dependent variable.
 
 - Residual Covariances: Adds covariances between the residuals of variables indicating that they have a common cause not specified in the model. Also adds variances for independent (exogenous) variables. This adds degrees of freedom corresponding to the number of covariances and variances added.
 	- Independent variables: Adds residual (co-) variances for independent (exogenous) variables. If the model includes an interaction term between a moderator and a mediator (e.g., second stage moderated mediation), residual covariances between the (error term of the) mediator, the moderator and the interaction term are added (see Kwan & Chan, 2018).

From ea1997ed744f07bc764c8f91a0d5fa96a7cdf399 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Franti=C5=A1ek=20Barto=C5=A1?= <f.bartos96@gmail.com>
Date: Sat, 7 Nov 2026 22:06:09 +0100
Subject: [PATCH 7/7] fix: flag model summaries whose information criteria are
 not comparable

With "Independent variables" selected under "Residual Covariances" (the
default), lavaan and blavaan fit with fixed.x = FALSE, so the likelihood
also describes the independent variables, moderators, and covariates.
Models with different covariates then describe different data, and their
AIC, BIC, and log-likelihoods cannot be compared. In a simulation where a
covariate truly predicts M and Y, adding it raised AIC by 484 (lm(): a
decrease of 93), and rescaling the covariate changed that difference.
Different cases have the same effect. Without independent covariances the
likelihood is conditional on the predictors, as in lm(), so models that
differ only in their predictors or covariates remain comparable.

The model summary tables now add a footnote when the compared models
differ in the variables whose distribution they model or in their cases.

The Bayesian analysis also stopped with an error from loo_compare() when
the models had different numbers of cases: without independent
covariances, blavaan removes cases with missing predictors, for example
when only one model has a covariate with missing values. The LOO
difference is now left out in that case.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
---
 NEWS.md                                       |  2 +
 R/bayesianProcess.R                           | 11 ++++-
 R/classicProcess.R                            | 24 +++++++++++
 inst/help/BayesianProcess.md                  |  2 +-
 inst/help/ClassicProcess.md                   |  2 +-
 .../testthat/test-bayesian-process-general.R  | 29 +++++++++++++
 .../test-classic-process-verification.R       | 41 +++++++++++++++++++
 7 files changed, 108 insertions(+), 3 deletions(-)

diff --git a/NEWS.md b/NEWS.md
index 13743910..406b059f 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -18,12 +18,14 @@
 
 ## Added
 * Covariates box for the "Hayes configuration" input type. As in PROCESS, covariates are added to every mediator and outcome equation, so the path coefficients, effects, and R² are adjusted for them. Covariates can be continuous or categorical, and they are left out of the mediation and total effects tables and the conceptual path plot ([Issue #63](https://github.com/jasp-stats/jaspProcess/issues/63), [Issue #4388](https://github.com/jasp-stats/jasp-issues/issues/4388), [Issue #3481](https://github.com/jasp-stats/jasp-issues/issues/3481)).
+* The model summary table adds a footnote when the compared models differ in the variables they describe or in the cases used, because their information criteria are then not comparable. With `Independent variables` selected under `Residual Covariances` (the default), this includes models with different independent variables, moderators, or covariates, because the likelihood then also describes these variables.
 
 ## Fixed
 * 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)).
 * With listwise deletion, mean-centering of moderators and the moderator probe values now use each model's own complete cases. Previously, missing values in any other variable in the variable boxes (for example a covariate used only in another model) restricted the centering sample and changed the estimates of models with complete data.
 * The Bayesian analysis now fits the mean-centered data when "Mean-centered moderation" is selected (the default). Previously only the moderator probe values were centered, so conditional effects were evaluated at moderator values shifted by the moderator mean.
 * Factor levels that contain two consecutive underscores or end with an underscore now stop the analysis with a message asking to rename them, instead of causing an unrelated error later.
+* The Bayesian analysis no longer stops with an error when the models have different numbers of cases, for example without independent covariances and with a covariate with missing values in only one model. The LOO difference is then left out of the model summary table.
 
 ---
 
diff --git a/R/bayesianProcess.R b/R/bayesianProcess.R
index 6582c986..3b2dfb76 100644
--- a/R/bayesianProcess.R
+++ b/R/bayesianProcess.R
@@ -259,6 +259,10 @@ BayesianProcess <- function(jaspResults, dataset = NULL, options) {
     summaryTable[["wBIC"]] <- .computeWeights(bic)
   }
 
+  if (length(procResults) > 1 && !.procModelsShareData(procResults)) {
+    summaryTable$addFootnote(message = .procIncomparableModelsMsg())
+  }
+
   logLikResults <- lapply(procResults, .procBayesCalcModelCaseLogLik)
 
   # Try to calc loo, will fail if LL contains NAs
@@ -274,7 +278,12 @@ BayesianProcess <- function(jaspResults, dataset = NULL, options) {
 
   waicResults <- lapply(logLikResults, function(mod) try(loo::waic(mod)))
 
-  if (length(looResults[looIsValid]) > 1) {
+  # loo_compare() stops when the models have different numbers of cases. Without
+  # independent covariances (fixed.x = TRUE), blavaan removes cases with missing
+  # predictors, so a covariate with missing values in only one model leads to this
+  looNobs <- sapply(looResults[looIsValid], function(res) nrow(res$pointwise))
+
+  if (length(looResults[looIsValid]) > 1 && all(looNobs == looNobs[1])) {
     ovtDiff <- gettext("LOO Difference")
 
     summaryTable$addColumnInfo(name = "elpdDiff", title = gettext("ELPD"), type = "number", overtitle = ovtDiff)
diff --git a/R/classicProcess.R b/R/classicProcess.R
index 4a286524..18b62555 100644
--- a/R/classicProcess.R
+++ b/R/classicProcess.R
@@ -1949,6 +1949,24 @@ ClassicProcess <- function(jaspResults, dataset = NULL, options) {
   }
 }
 
+.procModeledVariables <- function(fit) {
+  # Variables whose distribution enters the likelihood. With independent covariances
+  # (fixed.x = FALSE, the default) the likelihood also includes the exogenous
+  # variables (X, moderators, covariates, interaction terms); otherwise it is
+  # conditional on them, as in a regression.
+  what <- if (fit@Options[["fixed.x"]]) "ov.nox" else "ov"
+  return(sort(lavaan::lavNames(fit, what)))
+}
+
+.procModelsShareData <- function(fits) {
+  # Information criteria and log-likelihoods are only comparable between models
+  # of the same data: the same modeled variables and the same cases. Models with
+  # different predictors or covariates qualify only without independent covariances.
+  modeledVars <- lapply(fits, .procModeledVariables)
+  caseIdx     <- lapply(fits, lavaan::lavInspect, what = "case.idx")
+  return(all(sapply(modeledVars, identical, modeledVars[[1]])) && all(sapply(caseIdx, identical, caseIdx[[1]])))
+}
+
 .procModelSummaryTable <- function(jaspResults, options, modelsContainer, errors) {
   if (!is.null(jaspResults[["modelSummaryTable"]])) return()
 
@@ -2033,6 +2051,10 @@ ClassicProcess <- function(jaspResults, dataset = NULL, options) {
     if (options[["bicWeights"]]) {
       summaryTable[["wBIC"]][converged] <- .computeWeights(bic)
     }
+
+    if (sum(converged) > 1 && !.procModelsShareData(procResults[converged])) {
+      summaryTable$addFootnote(message = .procIncomparableModelsMsg())
+    }
   }
   if (options$estimator %in% c("dwls", "gls", "wls", "uls")) {
     summaryTable$addFootnote(message = .procInformationCriteriaMsg())
@@ -3359,6 +3381,8 @@ ClassicProcess <- function(jaspResults, dataset = NULL, options) {
 
 .procInformationCriteriaMsg <- function() gettext("The AIC, BIC and additional information criteria are only available with ML-type estimators.")
 
+.procIncomparableModelsMsg <- function() gettext("The models differ in the variables they describe or in the cases used, so their information criteria and log-likelihoods cannot be compared. With 'Independent variables' selected under 'Residual Covariances', a model also describes its independent variables, moderators, and covariates; otherwise only its mediators and dependent variables.")
+
 .procCausalAssumptionsMsg <- function() gettext(
   "<b>Important</b>: Parameter estimates can only be interpreted as causal effects if all confounding effects are accounted for and if the causal effect directions are correctly specified."
 )
diff --git a/inst/help/BayesianProcess.md b/inst/help/BayesianProcess.md
index ab22a267..e2cbcec4 100644
--- a/inst/help/BayesianProcess.md
+++ b/inst/help/BayesianProcess.md
@@ -60,7 +60,7 @@ This section allows users to specify multiple process models through one of two
 		- Moderator Z: Second moderator variable.
 		- Covariates: Control variables (continuous or categorical). As in the PROCESS macro, every covariate is added as a predictor to each mediator and outcome equation (paths labelled f into mediators and g into the outcome), so all path coefficients and effects are adjusted for the covariates. Covariates do not appear in the mediation effects tables or the conceptual path plot, and they do not change the Hayes configuration number. A covariate cannot also be the independent variable, a mediator, a moderator, or the dependent variable.
 
-- Residual Covariances: Adds covariances between the residuals of variables indicating that they have a common cause not specified in the model. Also adds variances for independent (exogenous) variables. This adds degrees of freedom corresponding to the number of covariances and variances added.
+- Residual Covariances: Adds covariances between the residuals of variables indicating that they have a common cause not specified in the model. Also adds variances for independent (exogenous) variables. This adds degrees of freedom corresponding to the number of covariances and variances added. With `Independent variables` selected, the independent variables, moderators, and covariates become part of the modeled data, so information criteria can then only be compared between models with the same variables. The model summary table adds a footnote when the compared models differ in the variables they describe or in the cases used.
 - Parameter Estimates: Determines which parameter estimates are shown in the output. Requesting intercepts adds the mean structure to the model, which adds parameters and will therefore affect some fit indices (such as the information criteria AIC and BIC).
 - 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.
diff --git a/inst/help/ClassicProcess.md b/inst/help/ClassicProcess.md
index ea835121..67ed6fc3 100644
--- a/inst/help/ClassicProcess.md
+++ b/inst/help/ClassicProcess.md
@@ -60,7 +60,7 @@ This section allows users to specify multiple process models through one of two
 		- Moderator Z: Second moderator variable.
 		- Covariates: Control variables (continuous or categorical). As in the PROCESS macro, every covariate is added as a predictor to each mediator and outcome equation (paths labelled f into mediators and g into the outcome), so all path coefficients and effects are adjusted for the covariates. Covariates do not appear in the mediation effects tables or the conceptual path plot, and they do not change the Hayes configuration number. A covariate cannot also be the independent variable, a mediator, a moderator, or the dependent variable.
 
-- Residual Covariances: Adds covariances between the residuals of variables indicating that they have a common cause not specified in the model. Also adds variances for independent (exogenous) variables. This adds degrees of freedom corresponding to the number of covariances and variances added.
+- Residual Covariances: Adds covariances between the residuals of variables indicating that they have a common cause not specified in the model. Also adds variances for independent (exogenous) variables. This adds degrees of freedom corresponding to the number of covariances and variances added. With `Independent variables` selected, the independent variables, moderators, and covariates become part of the modeled data, so information criteria can then only be compared between models with the same variables. The model summary table adds a footnote when the compared models differ in the variables they describe or in the cases used.
 	- Independent variables: Adds residual (co-) variances for independent (exogenous) variables. If the model includes an interaction term between a moderator and a mediator (e.g., second stage moderated mediation), residual covariances between the (error term of the) mediator, the moderator and the interaction term are added (see Kwan & Chan, 2018).
 	- Mediators: Adds residual (co-) variances for mediator variables.
 	- Dependent variables: Adds residual (co-) variances for dependent variables.
diff --git a/tests/testthat/test-bayesian-process-general.R b/tests/testthat/test-bayesian-process-general.R
index ebc21678..017c1e10 100644
--- a/tests/testthat/test-bayesian-process-general.R
+++ b/tests/testthat/test-bayesian-process-general.R
@@ -166,3 +166,32 @@ test_that("Mean-centered moderation evaluates conditional effects at the request
   expect_length(effects, 3)
   expect_lt(max(abs(effects - refEffects)), 0.03)
 })
+
+test_that("Model summary works for models fitted to different numbers of cases", {
+  skip_on_ci()
+  # Without independent covariances, cases with a missing predictor are removed, here
+  # only from Model 2, which has a covariate with missing values. loo_compare() stops for
+  # models with different numbers of cases, so the LOO difference must be left out, and
+  # a footnote says that the information criteria are not comparable.
+  options <- getOptionsBayesian()
+  options$covariates <- c(options$covariates, list("debMiss30"))
+  hayesModel4 <- function(name, covariates = list()) {
+    model <- getProcessModel(list(), name = name)
+    model$inputType <- "inputModelNumber"
+    model$modelNumber <- 4
+    model$modelNumberIndependent <- "contGamma"
+    model$modelNumberMediators <- list("debCollin1")
+    model$modelNumberCovariates <- covariates
+    model$independentCovariances <- FALSE
+    return(model)
+  }
+  options$processModels <- list(hayesModel4("Model 1"), hayesModel4("Model 2", list("debMiss30")))
+  set.seed(1)
+  results <- jaspTools::runAnalysis("BayesianProcess", "debug.csv", options)
+
+  summaryTable <- results[["results"]][["modelSummaryTable"]]
+  expect_equal(sapply(summaryTable[["data"]][1:2], function(row) row[["N"]]), c(100, 70))
+  expect_false(any(sapply(summaryTable[["data"]], function(row) "elpdDiff" %in% names(row))))
+  footnotes <- sapply(summaryTable[["footnotes"]], function(fn) fn[["text"]])
+  expect_true(jaspProcess:::.procIncomparableModelsMsg() %in% footnotes)
+})
diff --git a/tests/testthat/test-classic-process-verification.R b/tests/testthat/test-classic-process-verification.R
index 55e4dd75..a7e7675a 100644
--- a/tests/testthat/test-classic-process-verification.R
+++ b/tests/testthat/test-classic-process-verification.R
@@ -218,3 +218,44 @@ test_that("Test that moderators are centered and probed on each model's own comp
                  unname(refCoefs["contGamma"] + refCoefs["contGamma:contcor1"] * probes), tolerance = 1e-6)
   }
 })
+
+test_that("Test that the model summary flags information criteria of models fitted to different data", {
+  # Information criteria are only comparable between models of the same data. With
+  # independent covariances (the default), the likelihood also includes the predictors
+  # and covariates, so adding a covariate changes the data. Without them, the likelihood
+  # is conditional on the predictors, as in lm(), so only the cases must be the same.
+  runModels <- function(covariate, independentCovariances) {
+    options <- getOptionsClassical()
+    options$covariates <- c(options$covariates, list("debMiss30"))
+    options$naAction <- "listwise"
+    options$standardizedModelEstimates <- FALSE
+    model1 <- getHayesModelCovariates(4, "contGamma", mediators = list("debCollin1"))
+    model2 <- getHayesModelCovariates(4, "contGamma", mediators = list("debCollin1"), covariates = list(covariate))
+    model2$name <- "Model 2"
+    model1$independentCovariances <- model2$independentCovariances <- independentCovariances
+    options$processModels <- list(model1, model2)
+    set.seed(1)
+    return(jaspTools::runAnalysis("ClassicProcess", "debug.csv", options))
+  }
+  hasIncomparableFootnote <- function(results) {
+    footnotes <- sapply(results[["results"]][["modelSummaryTable"]][["footnotes"]], function(fn) fn[["text"]])
+    return(jaspProcess:::.procIncomparableModelsMsg() %in% footnotes)
+  }
+
+  # Without independent covariances, adding a covariate changes AIC as in lm()
+  results <- runModels("contcor2", independentCovariances = FALSE)
+  expect_false(hasIncomparableFootnote(results))
+  data <- read.csv(system.file("jaspData", "debug.csv", package = "jaspTools"))
+  refAic <- c(
+    AIC(lm(debCollin1 ~ contGamma, data = data)) + AIC(lm(contNormal ~ contGamma + debCollin1, data = data)),
+    AIC(lm(debCollin1 ~ contGamma + contcor2, data = data)) + AIC(lm(contNormal ~ contGamma + debCollin1 + contcor2, data = data))
+  )
+  aic <- sapply(results[["results"]][["modelSummaryTable"]][["data"]], function(row) row[["AIC"]])
+  expect_equal(diff(aic), diff(refAic), tolerance = 1e-6)
+
+  # With independent covariances, the models describe different variables
+  expect_true(hasIncomparableFootnote(runModels("contcor2", independentCovariances = TRUE)))
+
+  # A covariate with missing values removes cases from Model 2 only
+  expect_true(hasIncomparableFootnote(runModels("debMiss30", independentCovariances = FALSE)))
+})