diff --git a/R/bayesianProcess.R b/R/bayesianProcess.R
index 4c1e277..056e738 100644
--- a/R/bayesianProcess.R
+++ b/R/bayesianProcess.R
@@ -103,11 +103,15 @@ BayesianProcess <- function(jaspResults, dataset = NULL, options) {
   nuPrior   <- sprintf("normal(%s,%s)", options$nuPriorMu, options$nuPriorSigma)
   betaPrior <- sprintf("normal(%s,%s)", options$betaPriorMu, options$betaPriorSigma)
   psiPrior  <- sprintf("gamma(%s,%s)[sd]", options$psiPriorAlpha, options$psiPriorBeta)
-  rhoPrior  <- sprintf("beta(%s,%s)[sd]", options$rhoPriorAlpha, options$rhoPriorBeta)
+  rhoPrior  <- sprintf("beta(%s,%s)", options$rhoPriorAlpha, options$rhoPriorBeta)
   return(blavaan::dpriors(nu = nuPrior, beta = betaPrior, psi = psiPrior, rho = rhoPrior))
 }
 
 .procBayesResultsFitModel <- function(container, dataset, options, modelOptions) {
+  # Somehow the future.apply dependency of blavaan changes this option globally to NULL
+  # which throws an error; thus we change it locally
+  rlang::local_options(future.globals.method.default = "ordered")
+
   # Check if graph has error message
   if (!.procCheckGraph(container[["graph"]]$object) && jaspBase::isTryError(container[["graph"]]$object)) {
     return(.procEstimationMsg(container[["graph"]]$object))
@@ -152,7 +156,7 @@ BayesianProcess <- function(jaspResults, dataset = NULL, options) {
   }
   
   if (jaspBase::isTryError(fittedModel)) {
-    return(.procLavaanMsg(fittedModel))
+    return(.procEstimationMsg())
   }
 
   return(fittedModel)
diff --git a/R/classicProcess.R b/R/classicProcess.R
index 248c997..8254198 100644
--- a/R/classicProcess.R
+++ b/R/classicProcess.R
@@ -958,42 +958,33 @@ ClassicProcess <- function(jaspResults, dataset = NULL, options) {
   return(modOut)
 }
 
-.procMedEffectsSyntaxModPars <- function(pathEdge, sourceNode, contrFacVars, graph, modProbes) {
+.procMedEffectsSyntaxModPars <- function(pathEdge, sourceNode, contrFacVars, graph) {
   # Get moderator parameters for two-way interactions
   modPars <- lapply(pathEdge$modVars[[1]], function(v) {
     # Get edge for two way interaction between X and M
     twoWayEdge <- igraph::E(graph)[paste(sourceNode, v, sep = "__") %--% pathEdge$target]$parName
 
-    # Concatenate two way edge parName with moderator probes
-    return(apply(expand.grid(
-      twoWayEdge,
-      format(modProbes[[v]], digits = 3)
-    ), 1, paste, collapse = "*"))
+    # Concatenate two way edge parName with moderator
+    return(paste(twoWayEdge, v, sep = "*"))
   })
   
   # Concatenate dummy variables for factor moderators
   modPars <- .procMedEffectsSyntaxMergeFacMods(pathEdge$modVars[[1]], contrFacVars, modPars, sep = " + ")
   
-  # Get all possible combinations of probes from different moderators
-  modPars <- apply(expand.grid(modPars), 1, paste, collapse = " + ")
-  
   # Get name of potential three-way interaction
   threeWayInt <- paste(c(pathEdge$source, pathEdge$modVars[[1]]), collapse = "__")
 
   if (length(pathEdge$modVars[[1]]) > 1 && threeWayInt %in% igraph::E(graph)$source) { # If three-way int
     # Get edge of three way interaction X x M1 x M2
     threeWayEdge <- igraph::E(graph)[threeWayInt %--% pathEdge$target]
-    # Combine three way int parName with moderator probes
+    # Combine three way int parName with moderator
     threeWayModPars <- paste(
       threeWayEdge$parName,
-      apply(
-        expand.grid(lapply(pathEdge$modVars[[1]], function(v) format(modProbes[[v]], digits = 3))),
-        1, paste, collapse = "*"
-      ),
+      paste(pathEdge$modVars[[1]], collapse = "*"),
       sep = "*"
     )
-    # Add to previous moderator probes
-    modPars <- paste(modPars, threeWayModPars, sep = "+")
+    # Add to previous moderator terms
+    modPars <- append(modPars, threeWayModPars)
   }
   return(modPars)
 }
@@ -1067,8 +1058,11 @@ ClassicProcess <- function(jaspResults, dataset = NULL, options) {
     sourceFacVars <- numeric(0)
   }
 
+  rhs <- c()
+  uniqueMods <- c()
+
   # Right hand side of lavaan syntax
-  rhs <- lapply(2:length(path), function(i) {
+  for (i in 2:length(path)) {
     # Get edge from path[-1] to path[i]
     if (any(sourceInContrFacVars) && names(path)[i-1] %in% sourceFacVars) {
       # Get all vars with factor levels as source nodes
@@ -1081,24 +1075,141 @@ ClassicProcess <- function(jaspResults, dataset = NULL, options) {
     # Edge from previous node to current node
     pathEdge <- igraph::E(graph)[sourceNode %--% names(path)[i]]
 
+    modVars <- unlist(pathEdge$modVars)
+
     # If no moderators on edge, return only parName
-    if(any(is.na(unlist(pathEdge$modVars)))) return(pathEdge$parName)
+    if(any(is.na(modVars))) {
+      pathEdgeRhs <- pathEdge$parName
 
-    # Get pars from moderators
-    modPars <- .procMedEffectsSyntaxModPars(pathEdge, sourceNode, contrFacVars, graph, modProbes)
-    
-    # If indirect path add parentheses
-    if (i > 1) {
-      return(paste0("(", pathEdge$parName, " + ", modPars, ")"))
+      if (length(rhs) > 0) {
+        # Multiply effects on indirect paths
+        rhs <- paste(rhs, pathEdgeRhs, sep = "*")
+      } else {
+        rhs <- pathEdgeRhs
+      }
+    } else {
+      # Get pars from moderators
+      modPars <- .procMedEffectsSyntaxModPars(pathEdge, sourceNode, contrFacVars, graph)
+
+      # Get sum of moderator terms
+      modPars <- .doCallPaste(modPars, sep = " + ")
+
+      # If indirect path add parentheses
+      if (i > 1) {
+        pathEdgeRhs <- paste0("(", pathEdge$parName, " + ", unlist(modPars), ")")
+      } else {
+        pathEdgeRhs <- paste(pathEdge$parName, unlist(modPars), sep = " + ")
+      }
+
+      if (length(rhs) > 0) {
+        # Get all possible combinations with new moderators
+        if (all(!modVars %in% uniqueMods)) {
+          rhs <- apply(expand.grid(rhs, pathEdgeRhs), 1, paste, collapse = "*")
+        } else {
+          rhs <- paste(rhs, pathEdgeRhs, sep = "*")
+        }
+      } else {
+        rhs <- pathEdgeRhs
+      }
+      uniqueMods <- unique(c(uniqueMods, modVars))
     }
-    # Concanenate path edge parName with moderator probes
-    return(paste(pathEdge$parName, modPars, sep = " + "))
-  })
-  
-  # If indirect paths, multiply their steps
-  rhs <- .doCallPaste(rhs, sep = "*")
+  }
 
-  return(rhs)
+  modProbes <- modProbes[uniqueMods]
+
+  # Substitute moderator variables with probe values
+  rhsProbes <- c()
+
+  if (any(sourceInContrFacVars) && length(modProbes) > 0) {
+    # Identify which modProbes are factor dummies
+    modIdx <- sapply(names(modProbes), function(nm) {
+      idx <- which(sapply(contrFacVars, function(fac) nm %in% fac))
+
+      if (length(idx) == 0) return(FALSE)
+
+      return(idx)
+    })
+    # Remove non-factor moderators that have index zero
+    facModIdx <- modIdx[modIdx > 0]
+
+    facProbeList <- list()
+    facGroups <- list()
+    if (any(facModIdx)) {
+      # Group factor dummies by factor (with >2 levels)
+      facGroups <- contrFacVars[unique(facModIdx)]
+      names(facGroups) <- as.character(unique(facModIdx))
+      # Concatenate levels of the same factor
+      facProbeList <- lapply(facGroups, function(grp) {
+        probes <- Reduce(cbind, modProbes[grp])
+
+        if (is.null(dim(probes))) return(probes)
+
+        return(apply(probes, 1, paste, collapse = "_"))
+      })
+    }
+    # Now, collect all non-factor moderators
+    nonFacMods <- names(modProbes)[modIdx == 0]
+    nonFacProbeList <- modProbes[nonFacMods]
+    # Combine all for expansion
+    expandList <- c(facProbeList, nonFacProbeList)
+    modGrid <- expand.grid(expandList, stringsAsFactors = FALSE)
+    # Assign column names for substitution
+    colNames <- c(names(facProbeList), nonFacMods)
+    colnames(modGrid) <- colNames
+
+    # Split concatenated factor columns back into original dummy columns
+    if (length(facProbeList) > 0) {
+      for (k in 1:length(facProbeList)) {
+        facCol <- names(facProbeList)[k]
+        # Get the original dummy variable names for this factor
+        origDummies <- facGroups[[facCol]]
+        if (length(origDummies) == 1) {
+          splitVals <- matrix(modGrid[[facCol]], ncol = 1)
+        } else {
+          # Split the concatenated string into separate columns
+          splitVals <- do.call(rbind, strsplit(modGrid[[facCol]], "_", fixed = TRUE))
+          # If only one dummy, keep as vector
+        }
+        
+        # Assign columns with original dummy names
+        for (d in 1:length(origDummies)) {
+          modGrid[[origDummies[d]]] <- splitVals[, d]
+        }
+        # Remove the concatenated column
+        modGrid[[facCol]] <- NULL
+      }
+    }
+
+    for (j in 1:nrow(modGrid)) { # Probe combinations
+      for (i in 1:length(rhs)) { # Path edges
+        rowRhs <- rhs[i]
+        # Replace factor dummies as a group
+        for (p in 1:ncol(modGrid)) {
+          rowRhs <- gsub(colnames(modGrid)[p], format(modGrid[j, p], digits = 3), rowRhs)
+        }
+        rhsProbes <- c(rhsProbes, rowRhs)
+      }
+    }
+  } else {
+    # Default: iterate over rhs first, then probes
+    for (i in 1:length(rhs)) { # Path edges
+      modInPath <- sapply(names(modProbes), function(v) grepl(v, rhs[i]))
+      if (any(modInPath)) {
+        modGrid <- expand.grid(modProbes[modInPath])
+        for (j in 1:nrow(modGrid)) { # Probe combinations
+          rowRhs <- rhs[i]
+          for (p in 1:ncol(modGrid)) { # Moderator vars
+            rowRhs <- gsub(colnames(modGrid)[p], format(modGrid[j, p], digits = 3), rowRhs)
+          }
+          rhsProbes <- c(rhsProbes, rowRhs)
+        }
+      } else {
+        rhsProbes <- c(rhsProbes, rhs[i])
+      }
+    }
+  }
+
+  return(rhsProbes)
 }
 
 .procMedEffectsSyntaxSinglePath <- function(path, graph, modProbes, contrasts) {
@@ -1490,7 +1601,7 @@ ClassicProcess <- function(jaspResults, dataset = NULL, options) {
 .procResultsFitModel <- function(container, dataset, options, modelOptions) {
   # Check if graph has error message
   if (!.procCheckGraph(container[["graph"]]$object) && jaspBase::isTryError(container[["graph"]]$object)) {
-    return(.procEstimationMsg())
+    return(jaspBase::.extractErrorMessage(container[["graph"]]$object))
   }
 
   # Should model be fitted?
diff --git a/tests/testthat/test-classic-process-integration-general.R b/tests/testthat/test-classic-process-integration-general.R
index 99dcdd0..6593e02 100644
--- a/tests/testthat/test-classic-process-integration-general.R
+++ b/tests/testthat/test-classic-process-integration-general.R
@@ -1405,7 +1405,7 @@ test_that("Directed acyclic graph error message works", {
   set.seed(1)
   results <- jaspTools::runAnalysis("ClassicProcess", "debug.csv", options)
 
-  refMsg <- jaspProcess:::.procEstimationMsg(jaspProcess:::.procDagMsg())
+  refMsg <- jaspProcess:::.procDagMsg()
 
   msg <- results[["results"]][["parEstContainer"]][["collection"]][["parEstContainer_Model 1"]][["collection"]][["parEstContainer_Model 1_pathCoefficientsTable"]][["error"]][["errorMessage"]]
   expect_equal(msg, refMsg)
diff --git a/tests/testthat/test-classic-process-unit.R b/tests/testthat/test-classic-process-unit.R
index 6a7470f..9ff13f4 100644
--- a/tests/testthat/test-classic-process-unit.R
+++ b/tests/testthat/test-classic-process-unit.R
@@ -480,8 +480,8 @@ test_that("Test that .procMedEffectsSyntaxModPars works - no contrasts", {
   pathEdge <- igraph::E(graph)["contGamma" %--% "contNormal"]
   sourceNode <- "contGamma"
 
-  modPars <- jaspProcess:::.procMedEffectsSyntaxModPars(pathEdge, sourceNode, list(), graph, modProbes)
-  expect_equal(modPars, c("c3*0.1", "c3*0.5", "c3*0.9"))
+  modPars <- jaspProcess:::.procMedEffectsSyntaxModPars(pathEdge, sourceNode, list(), graph)
+  expect_equal(modPars, list("c3*contcor2"))
 })
 
 test_that("Test that .procMedEffectsSyntaxModPars works - with contrasts", {
@@ -502,8 +502,8 @@ test_that("Test that .procMedEffectsSyntaxModPars works - with contrasts", {
   pathEdge <- igraph::E(graph)["contGammaA" %--% "contNormal"]
   sourceNode <- "contGammaA"
 
-  modPars <- jaspProcess:::.procMedEffectsSyntaxModPars(pathEdge, sourceNode, contrFacVars, graph, modProbes)
-  expect_equal(modPars, c("c4*0.1", "c4*0.5", "c4*0.9"))
+  modPars <- jaspProcess:::.procMedEffectsSyntaxModPars(pathEdge, sourceNode, contrFacVars, graph)
+  expect_equal(modPars, list("c4*contcor2"))
 })
 
 test_that("Test that .procMedEffectsSyntaxGetLhs works - no contrasts", {