diff --git a/NEWS.md b/NEWS.md
index 6e707bfb..406b059f 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -16,8 +16,16 @@
 
 # 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)).
+* 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 a4633ab3..3b2dfb76 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
 
@@ -253,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
@@ -268,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 3d85954b..18b62555 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
 
@@ -493,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"]]
 
@@ -559,6 +615,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 +635,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 +667,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 +698,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 +710,7 @@ ClassicProcess <- function(jaspResults, dataset = NULL, options) {
           # Else -> c
           "c")
       )
+    )
   )
   # Get occurence of each symbol
   parSymbolTable <- table(igraph::E(graph)$parSymbol)
@@ -703,6 +773,22 @@ ClassicProcess <- function(jaspResults, dataset = NULL, options) {
       }
     }
   }
+
+  .procCheckFactorLevels(dataset, factors)
+}
+
+.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) {
@@ -710,12 +796,8 @@ ClassicProcess <- function(jaspResults, dataset = NULL, options) {
   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"]]) {
@@ -776,6 +858,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)
@@ -788,12 +873,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]
 
@@ -1717,6 +1796,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 +1894,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
@@ -1863,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()
 
@@ -1947,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())
@@ -2953,6 +3061,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 +3078,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 +3130,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)
+        )
       )
     )
   )
@@ -3022,8 +3143,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"
@@ -3253,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 1905c98c..e2cbcec4 100644
--- a/inst/help/BayesianProcess.md
+++ b/inst/help/BayesianProcess.md
@@ -58,8 +58,9 @@ 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.
+- 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 df24c125..67ed6fc3 100644
--- a/inst/help/ClassicProcess.md
+++ b/inst/help/ClassicProcess.md
@@ -58,8 +58,9 @@ 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.
+- 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/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..017c1e10 100644
--- a/tests/testthat/test-bayesian-process-general.R
+++ b/tests/testthat/test-bayesian-process-general.R
@@ -105,3 +105,93 @@ 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)
+})
+
+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)
+})
+
+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-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)
+})
diff --git a/tests/testthat/test-classic-process-unit.R b/tests/testthat/test-classic-process-unit.R
index ed5b5be3..a8420cf0 100644
--- a/tests/testthat/test-classic-process-unit.R
+++ b/tests/testthat/test-classic-process-unit.R
@@ -835,3 +835,118 @@ 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"))
 })
+
+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() {
+  # 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..a7e7675a 100644
--- a/tests/testthat/test-classic-process-verification.R
+++ b/tests/testthat/test-classic-process-verification.R
@@ -82,3 +82,180 @@ 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, 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"]])))
+}
+
+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)
+})
+
+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)
+  }
+})
+
+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)))
+})