Add covariates to Hayes configuration

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.
 
 ---
 

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)

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,15 +698,19 @@ 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
         ifelse(igraph::V(graph)[igraph::E(graph)$source]$isExo & igraph::V(graph)[igraph::E(graph)$target]$isMed, "a",
           # Else -> c
           "c")
       )
+    )
   )
   # Get occurence of each symbol
   parSymbolTable <- table(igraph::E(graph)$parSymbol)
@@ -703,19 +773,31 @@ 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) {
   # 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"]]) {
@@ -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."
 )