diff --git a/R/predictRTConf.R b/R/predictRTConf.R
index 0cc6cb2..cc8a4d6 100644
--- a/R/predictRTConf.R
+++ b/R/predictRTConf.R
@@ -128,7 +128,6 @@ predictConf <- function(paramDf, model=NULL,
                            maxrt=Inf, subdivisions = 100L, simult_conf = FALSE, stop.on.error=FALSE,
                            .progress=TRUE){
   paramDf <- as.data.frame(paramDf)
-  paramDf <- paramDf[,!is.na(paramDf)]
   if (is.null(model) && ("model" %in% names(paramDf))) model <- paramDf$model
   if (model =="DDMConf") {
     warning("DDMConf was renamed DDConf in version 1.0.0! DDConf will be predicted!")
@@ -157,7 +156,6 @@ predictRT <- function(paramDf, model=NULL,
                          scaled = FALSE, DistConf=NULL,
                          .progress = TRUE) {
   paramDf <- as.data.frame(paramDf)
-  paramDf <- paramDf[,!is.na(paramDf)]
   if (is.null(model) && ("model" %in% names(paramDf))) model <- paramDf$model
 
   #### Check model argument
diff --git a/R/predictratingdist_DDConf.R b/R/predictratingdist_DDConf.R
index 1e5dec4..ddf8436 100644
--- a/R/predictratingdist_DDConf.R
+++ b/R/predictratingdist_DDConf.R
@@ -109,8 +109,9 @@
 #' @rdname predictDDConf
 #' @export
 predictDDConf_Conf <- function(paramDf,
-                                maxrt=Inf, subdivisions = 100L, stop.on.error=FALSE,
-                                .progress=TRUE){
+                                  maxrt=Inf, subdivisions = 100L, stop.on.error=FALSE,
+                                  .progress=TRUE){
+  paramDf <- fill_optional_params(paramDf, c(st0 = 0, sz = 0, sv = 0))
 
   nConds <- length(grep(pattern = "^v[0-9]", names(paramDf), value = T))
   symmetric_confidence_thresholds <- length(grep(pattern = "thetaUpper", names(paramDf), value = T))<1
@@ -144,11 +145,32 @@ predictDDConf_Conf <- function(paramDf,
   }
   ## Recover confidence thresholds
   if (symmetric_confidence_thresholds) {
-    thetas_upper <- c(0, t(paramDf[,paste("theta",(nRatings-1):1, sep = "")]), 1e+64)
-    thetas_lower <- c(0, t(paramDf[,paste("theta",(nRatings-1):1, sep = "")]), 1e+64)
+    theta_cols <- grep("^theta[0-9]", names(paramDf), value = TRUE)
+    theta_cols <- theta_cols[order(as.integer(sub("^theta", "", theta_cols)))]
+    if (length(theta_cols) == 0L) {
+      theta_vals <- numeric(0)
+    } else {
+      theta_vals <- t(paramDf[, rev(theta_cols), drop = FALSE])
+    }
+    thetas_upper <- c(0, theta_vals, 1e+64)
+    thetas_lower <- c(0, theta_vals, 1e+64)
   } else {
-    thetas_upper <- c(0, t(paramDf[,paste("thetaUpper",(nRatings-1):1, sep = "")]), 1e+64)
-    thetas_lower <- c(0, t(paramDf[,paste("thetaLower",(nRatings-1):1, sep="")]), 1e+64)
+    theta_upper_cols <- grep("^thetaUpper[0-9]", names(paramDf), value = TRUE)
+    theta_upper_cols <- theta_upper_cols[order(as.integer(sub("^thetaUpper", "", theta_upper_cols)))]
+    theta_lower_cols <- grep("^thetaLower[0-9]", names(paramDf), value = TRUE)
+    theta_lower_cols <- theta_lower_cols[order(as.integer(sub("^thetaLower", "", theta_lower_cols)))]
+    if (length(theta_upper_cols) == 0L) {
+      theta_upper_vals <- numeric(0)
+    } else {
+      theta_upper_vals <- t(paramDf[, rev(theta_upper_cols), drop = FALSE])
+    }
+    if (length(theta_lower_cols) == 0L) {
+      theta_lower_vals <- numeric(0)
+    } else {
+      theta_lower_vals <- t(paramDf[, rev(theta_lower_cols), drop = FALSE])
+    }
+    thetas_upper <- c(0, theta_upper_vals, 1e+64)
+    thetas_lower <- c(0, theta_lower_vals, 1e+64)
   }
   # Because we integrate over the response time, st0 does not matter
   # So, to speed up computations for high values of st0, we set it to 0
@@ -200,6 +222,7 @@ predictDDConf_RT <- function(paramDf,
                               maxrt=9, subdivisions = 100L,  minrt=NULL,
                               scaled = FALSE, DistConf=NULL,
                               .progress = TRUE) {
+  paramDf <- fill_optional_params(paramDf, c(st0 = 0, sz = 0, sv = 0))
   if (scaled && is.null(DistConf)) {
     message(paste("scaled is TRUE and DistConf is NULL. The rating distribution will",
                   " be computed, which will take additional time.", sep=""))
@@ -236,11 +259,32 @@ predictDDConf_RT <- function(paramDf,
   }
   ## Recover confidence thresholds
   if (symmetric_confidence_thresholds) {
-    thetas_upper <- c(0, t(paramDf[,paste("theta",(nRatings-1):1, sep = "")]), 1e+64)
-    thetas_lower <- c(0, t(paramDf[,paste("theta",(nRatings-1):1, sep = "")]), 1e+64)
+    theta_cols <- grep("^theta[0-9]", names(paramDf), value = TRUE)
+    theta_cols <- theta_cols[order(as.integer(sub("^theta", "", theta_cols)))]
+    if (length(theta_cols) == 0L) {
+      theta_vals <- numeric(0)
+    } else {
+      theta_vals <- t(paramDf[, rev(theta_cols), drop = FALSE])
+    }
+    thetas_upper <- c(0, theta_vals, 1e+64)
+    thetas_lower <- c(0, theta_vals, 1e+64)
   } else {
-    thetas_upper <- c(0, t(paramDf[,paste("thetaUpper",(nRatings-1):1, sep = "")]), 1e+64)
-    thetas_lower <- c(0, t(paramDf[,paste("thetaLower",(nRatings-1):1, sep="")]), 1e+64)
+    theta_upper_cols <- grep("^thetaUpper[0-9]", names(paramDf), value = TRUE)
+    theta_upper_cols <- theta_upper_cols[order(as.integer(sub("^thetaUpper", "", theta_upper_cols)))]
+    theta_lower_cols <- grep("^thetaLower[0-9]", names(paramDf), value = TRUE)
+    theta_lower_cols <- theta_lower_cols[order(as.integer(sub("^thetaLower", "", theta_lower_cols)))]
+    if (length(theta_upper_cols) == 0L) {
+      theta_upper_vals <- numeric(0)
+    } else {
+      theta_upper_vals <- t(paramDf[, rev(theta_upper_cols), drop = FALSE])
+    }
+    if (length(theta_lower_cols) == 0L) {
+      theta_lower_vals <- numeric(0)
+    } else {
+      theta_lower_vals <- t(paramDf[, rev(theta_lower_cols), drop = FALSE])
+    }
+    thetas_upper <- c(0, theta_upper_vals, 1e+64)
+    thetas_lower <- c(0, theta_lower_vals, 1e+64)
   }
 
   if (is.null(minrt)) minrt <- paramDf$t0
diff --git a/R/predictratingdist_MTLNR.R b/R/predictratingdist_MTLNR.R
index 2655d13..cd9c836 100644
--- a/R/predictratingdist_MTLNR.R
+++ b/R/predictratingdist_MTLNR.R
@@ -154,6 +154,7 @@ predictMTLNR_Conf <- function(paramDf,
                                 maxrt=Inf, subdivisions = 100L, stop.on.error=FALSE,
                                 .progress=TRUE){
 
+  paramDf <- fill_optional_params(paramDf, c(st0 = 0))
   nConds <- length(grep(pattern = "^v[0-9]", names(paramDf), value = T))
   symmetric_confidence_thresholds <- length(grep(pattern = "thetaUpper", names(paramDf), value = T))<1
   if (symmetric_confidence_thresholds) {
@@ -233,6 +234,7 @@ predictMTLNR_RT <- function(paramDf,
                                  scaled = FALSE, DistConf=NULL,
                                 .progress = TRUE) {
 
+  paramDf <- fill_optional_params(paramDf, c(st0 = 0))
   if (scaled && is.null(DistConf)) {
     message(paste("scaled is TRUE and DistConf is NULL. The rating distribution will",
     " be computed, which will take additional time.", sep=""))
diff --git a/R/predictratingdist_RM.R b/R/predictratingdist_RM.R
index a1ddd9d..b7b03d6 100644
--- a/R/predictratingdist_RM.R
+++ b/R/predictratingdist_RM.R
@@ -133,6 +133,7 @@
 predictRM_Conf <- function(paramDf, model="IRM", time_scaled = FALSE,
                                 maxrt=Inf, subdivisions = 100L, stop.on.error=FALSE,
                                 .progress=TRUE){
+  paramDf <- fill_optional_params(paramDf, c(st0 = 0, wx = 1, wrt = 0, wint = 0))
   #### Check model argument
   if (grepl("IRMt",model)) {
     time_scaled=TRUE
@@ -255,6 +256,7 @@ predictRM_RT <- function(paramDf, model="IRM", time_scaled = FALSE,
                                  maxrt=9, subdivisions = 100L,  minrt=NULL,
                                  scaled = FALSE, DistConf=NULL,
                                 .progress = TRUE) {
+  paramDf <- fill_optional_params(paramDf, c(st0 = 0, wx = 1, wrt = 0, wint = 0))
   #### Check model argument
   if (grepl("IRMt",model)) {
     time_scaled=TRUE
diff --git a/R/predictratingdist_WEV.R b/R/predictratingdist_WEV.R
index 40a20e2..3a6687f 100644
--- a/R/predictratingdist_WEV.R
+++ b/R/predictratingdist_WEV.R
@@ -150,6 +150,8 @@ predictWEV_Conf <- function(paramDf, model="dynaViTE",
   if (grepl("2DSD", model)) model <- "2DSD"
   if ("model" %in% names(paramDf)) paramDf$model <- NULL
 
+  paramDf <- fill_optional_params(paramDf, c(st0 = 0, sz = 0, sv = 0))
+
   nConds <- length(grep(pattern = "^v[0-9]", names(paramDf), value = T))
   symmetric_confidence_thresholds <- length(grep(pattern = "thetaUpper", names(paramDf), value = T))<1
   if (symmetric_confidence_thresholds) {
@@ -182,13 +184,34 @@ predictWEV_Conf <- function(paramDf, model="dynaViTE",
   }
   ## Recover confidence thresholds
   if (symmetric_confidence_thresholds) {
-    thetas_upper <- c(-1e+32, t(paramDf[,paste("theta",1:(nRatings-1), sep = "")]), 1e+32)
-    thetas_lower <- c(-1e+32, t(paramDf[,paste("theta",1:(nRatings-1), sep = "")]), 1e+32)
+    theta_cols <- grep("^theta[0-9]", names(paramDf), value = TRUE)
+    theta_cols <- theta_cols[order(as.integer(sub("^theta", "", theta_cols)))]
+    if (length(theta_cols) == 0L) {
+      theta_vals <- numeric(0)
+    } else {
+      theta_vals <- t(paramDf[, theta_cols, drop = FALSE])
+    }
+    thetas_upper <- c(-1e+32, theta_vals, 1e+32)
+    thetas_lower <- c(-1e+32, theta_vals, 1e+32)
   } else {
-    thetas_upper <- c(-1e+32, t(paramDf[,paste("thetaUpper",1:(nRatings-1), sep = "")]), 1e+32)
-    thetas_lower <- c(-1e+32, t(paramDf[,paste("thetaLower",1:(nRatings-1), sep="")]), 1e+32)
+    theta_upper_cols <- grep("^thetaUpper[0-9]", names(paramDf), value = TRUE)
+    theta_upper_cols <- theta_upper_cols[order(as.integer(sub("^thetaUpper", "", theta_upper_cols)))]
+    theta_lower_cols <- grep("^thetaLower[0-9]", names(paramDf), value = TRUE)
+    theta_lower_cols <- theta_lower_cols[order(as.integer(sub("^thetaLower", "", theta_lower_cols)))]
+    if (length(theta_upper_cols) == 0L) {
+      theta_upper_vals <- numeric(0)
+    } else {
+      theta_upper_vals <- t(paramDf[, theta_upper_cols, drop = FALSE])
+    }
+    if (length(theta_lower_cols) == 0L) {
+      theta_lower_vals <- numeric(0)
+    } else {
+      theta_lower_vals <- t(paramDf[, theta_lower_cols, drop = FALSE])
+    }
+    thetas_upper <- c(-1e+32, theta_upper_vals, 1e+32)
+    thetas_lower <- c(-1e+32, theta_lower_vals, 1e+32)
   }
-  if (thetas_lower[2]>thetas_lower[3]) {
+  if (length(thetas_lower) >= 3 && thetas_lower[2]>thetas_lower[3]) {
     # For 2DSD the parametrization for lower thetas is different (different confidence scale)
     thetas_lower <- c(-1e+32, rev(thetas_lower[2:(nRatings)]), 1e+32)
     if (symmetric_confidence_thresholds) {
@@ -277,6 +300,7 @@ predictWEV_RT <- function(paramDf, model=NULL,
   if (grepl("2DSD", model)) model <- "2DSD"
   if ("model" %in% names(paramDf)) paramDf$model <- NULL
 
+  paramDf <- fill_optional_params(paramDf, c(st0 = 0, sz = 0, sv = 0))
 
   nConds <- length(grep(pattern = "^v[0-9]", names(paramDf), value = T))
   symmetric_confidence_thresholds <- length(grep(pattern = "thetaUpper", names(paramDf), value = T))<1
@@ -311,13 +335,34 @@ predictWEV_RT <- function(paramDf, model=NULL,
 
   ## Recover confidence thresholds
   if (symmetric_confidence_thresholds) {
-    thetas_upper <- c(-1e+32, t(paramDf[,paste("theta",1:(nRatings-1), sep = "")]), 1e+32)
-    thetas_lower <- c(-1e+32, t(paramDf[,paste("theta",1:(nRatings-1), sep = "")]), 1e+32)
+    theta_cols <- grep("^theta[0-9]", names(paramDf), value = TRUE)
+    theta_cols <- theta_cols[order(as.integer(sub("^theta", "", theta_cols)))]
+    if (length(theta_cols) == 0L) {
+      theta_vals <- numeric(0)
+    } else {
+      theta_vals <- t(paramDf[, theta_cols, drop = FALSE])
+    }
+    thetas_upper <- c(-1e+32, theta_vals, 1e+32)
+    thetas_lower <- c(-1e+32, theta_vals, 1e+32)
   } else {
-    thetas_upper <- c(-1e+32, t(paramDf[,paste("thetaUpper",1:(nRatings-1), sep = "")]), 1e+32)
-    thetas_lower <- c(-1e+32, t(paramDf[,paste("thetaLower",1:(nRatings-1), sep="")]), 1e+32)
+    theta_upper_cols <- grep("^thetaUpper[0-9]", names(paramDf), value = TRUE)
+    theta_upper_cols <- theta_upper_cols[order(as.integer(sub("^thetaUpper", "", theta_upper_cols)))]
+    theta_lower_cols <- grep("^thetaLower[0-9]", names(paramDf), value = TRUE)
+    theta_lower_cols <- theta_lower_cols[order(as.integer(sub("^thetaLower", "", theta_lower_cols)))]
+    if (length(theta_upper_cols) == 0L) {
+      theta_upper_vals <- numeric(0)
+    } else {
+      theta_upper_vals <- t(paramDf[, theta_upper_cols, drop = FALSE])
+    }
+    if (length(theta_lower_cols) == 0L) {
+      theta_lower_vals <- numeric(0)
+    } else {
+      theta_lower_vals <- t(paramDf[, theta_lower_cols, drop = FALSE])
+    }
+    thetas_upper <- c(-1e+32, theta_upper_vals, 1e+32)
+    thetas_lower <- c(-1e+32, theta_lower_vals, 1e+32)
   }
-  if (thetas_lower[2]>thetas_lower[3]) {
+  if (length(thetas_lower) >= 3 && thetas_lower[2]>thetas_lower[3]) {
     # For 2DSD the parametrization for lower thetas is different (different confidence scale)
     thetas_lower <- c(-1e+32, rev(thetas_lower[2:(nRatings)]), 1e+32)
     if (symmetric_confidence_thresholds) {
diff --git a/R/simulateConfModel.R b/R/simulateConfModel.R
index 7db8e0b..c766249 100644
--- a/R/simulateConfModel.R
+++ b/R/simulateConfModel.R
@@ -101,7 +101,6 @@ simulateRTConf <- function (paramDf, n=1e+4,  model = NULL,
 {
   gc(verbose = FALSE, full=FALSE)
   if (nrow(paramDf)>1) stop("paramDf must have one row.")
-  paramDf <- paramDf[,c(!is.na(paramDf))]
   if (is.null(model) && ("model" %in% names(paramDf))) model <- paramDf$model
   if ((model %in% c("dynaViTE", "dynWEV", "WEVmu", "2DSD", "2DSDT")) && identical(stimulus, c(1,2))) stimulus <- c(-1,1)
   if (grepl("RM", model)) {
diff --git a/R/simulateMTLNR.R b/R/simulateMTLNR.R
index 621a6a8..9a8a000 100644
--- a/R/simulateMTLNR.R
+++ b/R/simulateMTLNR.R
@@ -165,6 +165,7 @@ simulateMTLNR <- function (paramDf, n=1e+4, gamma = FALSE, agg_simus=FALSE,
   if (!is.null(seed)) {
     set.seed(seed)
   }
+  paramDf <- fill_optional_params(paramDf, c(st0 = 0))
 
   if (!(all(stimulus %in% c(1, 2)))) {
     stop(paste("Not accepted value for stimulus: ", paste(stimulus, collapse=", "),". Must be either 1, 2, or c(1, 2).", sep=""))
diff --git a/R/simulateRM.R b/R/simulateRM.R
index 599bce0..51b0ded 100644
--- a/R/simulateRM.R
+++ b/R/simulateRM.R
@@ -155,6 +155,7 @@ simulateRM <- function (paramDf, n=1e+4,  model = "IRM", time_scaled=FALSE,
     model = "PCRM"
     time_scaled=TRUE
   }
+  paramDf <- fill_optional_params(paramDf, c(st0 = 0, wx = 1, wrt = 0, wint = 0))
   if (!model %in% c("IRM", "PCRM")) stop("model must be 'IRM', 'PCRM', 'IRMt' or 'PCRMt'")
 
   if (!(all(stimulus %in% c(1,2)))) {
@@ -430,4 +431,3 @@ rRM_Kiani <- function (paramDf, n=1e+4, time_scaled=FALSE,
 
 
 
-
diff --git a/R/simulateWEV.R b/R/simulateWEV.R
index 28aa8c4..e13e118 100644
--- a/R/simulateWEV.R
+++ b/R/simulateWEV.R
@@ -154,6 +154,7 @@ simulateWEV <- function (paramDf, n=1e+4,  model = "dynWEV", simult_conf = FALSE
   if (!is.null(seed)) {
     set.seed(seed)
   }
+  paramDf <- fill_optional_params(paramDf, c(st0 = 0, sz = 0, sv = 0))
   if (!("lambda" %in% names(paramDf))) {
     if (model %in% c("dynaViTE", "2DSDT")) warning("No lambda specified in paramDf. lambda=0 used")
     lambda <- 0
@@ -257,11 +258,32 @@ simulateWEV <- function (paramDf, n=1e+4,  model = "dynWEV", simult_conf = FALSE
   if (!process_results) simus <- select(simus, -c("dec", "vis", "mu"))
 
   if (symmetric_confidence_thresholds) {
-    thetas_upper <- c(-Inf, t(paramDf[,paste("theta",1:(nRatings-1), sep = "")]), Inf)
-    thetas_lower <- c(-Inf, t(paramDf[,paste("theta",1:(nRatings-1), sep = "")]), Inf)
+    theta_cols <- grep("^theta[0-9]", names(paramDf), value = TRUE)
+    theta_cols <- theta_cols[order(as.integer(sub("^theta", "", theta_cols)))]
+    if (length(theta_cols) == 0L) {
+      theta_vals <- numeric(0)
+    } else {
+      theta_vals <- t(paramDf[, theta_cols, drop = FALSE])
+    }
+    thetas_upper <- c(-Inf, theta_vals, Inf)
+    thetas_lower <- c(-Inf, theta_vals, Inf)
   } else {
-    thetas_upper <- c(-Inf, t(paramDf[,paste("thetaUpper",1:(nRatings-1), sep = "")]), Inf)
-    thetas_lower <- c(-Inf, t(paramDf[,paste("thetaLower",1:(nRatings-1), sep="")]), Inf)
+    theta_upper_cols <- grep("^thetaUpper[0-9]", names(paramDf), value = TRUE)
+    theta_upper_cols <- theta_upper_cols[order(as.integer(sub("^thetaUpper", "", theta_upper_cols)))]
+    theta_lower_cols <- grep("^thetaLower[0-9]", names(paramDf), value = TRUE)
+    theta_lower_cols <- theta_lower_cols[order(as.integer(sub("^thetaLower", "", theta_lower_cols)))]
+    if (length(theta_upper_cols) == 0L) {
+      theta_upper_vals <- numeric(0)
+    } else {
+      theta_upper_vals <- t(paramDf[, theta_upper_cols, drop = FALSE])
+    }
+    if (length(theta_lower_cols) == 0L) {
+      theta_lower_vals <- numeric(0)
+    } else {
+      theta_lower_vals <- t(paramDf[, theta_lower_cols, drop = FALSE])
+    }
+    thetas_upper <- c(-Inf, theta_upper_vals, Inf)
+    thetas_lower <- c(-Inf, theta_lower_vals, Inf)
   }
 
   levels_lower <- cumsum(as.numeric(table(thetas_lower)))
diff --git a/R/utils_param_defaults.R b/R/utils_param_defaults.R
new file mode 100644
index 0000000..8a9e394
--- /dev/null
+++ b/R/utils_param_defaults.R
@@ -0,0 +1,33 @@
+#' Fill optional parameter columns with default values.
+#'
+#' Ensures that optional scalar parameters exist and replaces `NA` entries with
+#' the provided defaults. This keeps the original column order intact (new
+#' columns are appended) and works for one-row data frames and list inputs that
+#' can be coerced to a data frame.
+#'
+#' @param paramDf Data frame or list specifying model parameters.
+#' @param defaults Named list or vector of default values.
+#'
+#' @return A data frame with missing/`NA` optional columns replaced.
+#' @keywords internal
+fill_optional_params <- function(paramDf, defaults) {
+  paramDf <- as.data.frame(paramDf, stringsAsFactors = FALSE)
+  if (length(defaults) == 0L || nrow(paramDf) == 0L) {
+    return(paramDf)
+  }
+
+  for (nm in names(defaults)) {
+    default_val <- defaults[[nm]]
+    if (!nm %in% names(paramDf)) {
+      paramDf[[nm]] <- rep_len(default_val, nrow(paramDf))
+      next
+    }
+
+    missing_idx <- is.na(paramDf[[nm]])
+    if (any(missing_idx)) {
+      paramDf[[nm]][missing_idx] <- rep_len(default_val, sum(missing_idx))
+    }
+  }
+
+  paramDf
+}
diff --git a/tests/testthat/test-RNGs.R b/tests/testthat/test-RNGs.R
index 459ec5e..e61af30 100644
--- a/tests/testthat/test-RNGs.R
+++ b/tests/testthat/test-RNGs.R
@@ -14,6 +14,39 @@ test_that("r2DSD works", {
   expect_true(all(r2DSD(20, tau=1, a=2, v=0.1, t0=1, z =1.5, sv=0.2, maxrt = 5, stop_on_error = FALSE)== 0))
 })
 
+
+test_that("simulation wrappers fill missing optional parameters", {
+  wev_params <- data.frame(
+    a = 1.5, v = 0.4, t0 = 0.2, z = 0.5, tau = 0.1,
+    w = 0.5, svis = 0.1, sigvis = 0.05, theta1 = 0.6
+  )
+  expect_silent(simulateWEV(wev_params, n = 10, model = "dynWEV", simult_conf = FALSE, gamma = FALSE, agg_simus = FALSE, stimulus = 1))
+
+  rm_params <- data.frame(a = 1.2, b = 1.3, v = 0.3, t0 = 0.2, theta1 = 0.5)
+  expect_silent(simulateRM(rm_params, n = 10, model = "IRM", gamma = FALSE, agg_simus = FALSE, stimulus = 1))
+
+  conf_params <- data.frame(model = "dynWEV", a = 1.5, v = 0.4, t0 = 0.2, z = 0.5, tau = 0.1, w = 0.5, svis = 0.1, sigvis = 0.05, theta1 = 0.6)
+  expect_silent(simulateRTConf(conf_params, n = 10, simult_conf = FALSE, stimulus = 1))
+})
+
+
+test_that("rLCA simult_conf only shifts RTs by tau", {
+  set.seed(123)
+  base <- rLCA(
+    n = 20, mu1 = 1, mu2 = 0.5, th1 = 1.2, th2 = 1.1,
+    tau = 0.2, t0 = 0.1, st0 = 0.05, simult_conf = FALSE
+  )
+  set.seed(123)
+  with_conf <- rLCA(
+    n = 20, mu1 = 1, mu2 = 0.5, th1 = 1.2, th2 = 1.1,
+    tau = 0.2, t0 = 0.1, st0 = 0.05, simult_conf = TRUE
+  )
+
+  expect_equal(with_conf$response, base$response)
+  expect_equal(with_conf$conf, base$conf)
+  expect_equal(with_conf$rt, base$rt + 0.2)
+})
+
 test_that("rdynaViTE works", {
   expect_true(all(rdynaViTE(200, tau=1, a=2, v=c(0.5,-0.5), t0=0, z =0.5, sv=0.2, w= 0.8)$response %in% c(-1,0,1)))
   expect_true(all(rdynaViTE(200, tau=1, a=2, v=0.5, t0=0, z =c(0.4,0.5, 0.6), sv=0.2, w= 0.8)$rt >= 0))
diff --git a/tests/testthat/test-prediction.R b/tests/testthat/test-prediction.R
index 9297816..34899bf 100644
--- a/tests/testthat/test-prediction.R
+++ b/tests/testthat/test-prediction.R
@@ -43,5 +43,101 @@ test_that("Prediction sums to 1", {
 
 
 
+# Regression coverage for optional parameters ---------------------------------
+
+test_that("model-specific predictors tolerate missing optional parameters", {
+  wev_params <- data.frame(
+    a = 1, v = 0.3, t0 = 0.2, z = 0.5, tau = 0.1, w = 0.5,
+    svis = 0.1, sigvis = 0.05
+  )
+  expect_silent(
+    predictWEV_Conf(
+      wev_params, model = "dynWEV", maxrt = 1, subdivisions = 10,
+      simult_conf = FALSE, stop.on.error = FALSE, precision = 2,
+      .progress = FALSE
+    )
+  )
+  expect_silent(
+    predictWEV_RT(
+      wev_params, model = "dynWEV", maxrt = 1, subdivisions = 10,
+      minrt = wev_params$t0, simult_conf = FALSE, scaled = FALSE,
+      DistConf = NULL, precision = 2, .progress = FALSE
+    )
+  )
+
+  rm_params <- data.frame(a = 1, b = 1.2, v = 0.4, t0 = 0.2, theta1 = 0.5)
+  expect_silent(
+    predictRM_Conf(
+      rm_params, model = "IRM", maxrt = 1, subdivisions = 10,
+      stop.on.error = FALSE, .progress = FALSE
+    )
+  )
+  expect_silent(
+    predictRM_RT(
+      rm_params, model = "IRM", maxrt = 1, subdivisions = 10,
+      minrt = rm_params$t0, scaled = FALSE, DistConf = NULL,
+      .progress = FALSE
+    )
+  )
+
+  mtl_params <- data.frame(
+    v = 0.4, t0 = 0.2,
+    mu_d1 = 0, mu_d2 = 0,
+    s_v1 = 1, s_v2 = 1,
+    s_d1 = 1, s_d2 = 1,
+    rho_v = 0, rho_d = 0
+  )
+  expect_silent(
+    predictMTLNR_Conf(
+      mtl_params, maxrt = 1, subdivisions = 10,
+      stop.on.error = FALSE, .progress = FALSE
+    )
+  )
+  expect_silent(
+    predictMTLNR_RT(
+      mtl_params, maxrt = 1, subdivisions = 10,
+      minrt = mtl_params$t0, scaled = FALSE, DistConf = NULL,
+      .progress = FALSE
+    )
+  )
+
+  dd_params <- data.frame(a = 1, v = 0.2, t0 = 0.2, z = 0.5, sz = NA_real_, sv = NA_real_)
+  expect_silent(
+    predictDDConf_Conf(
+      dd_params, maxrt = 1, subdivisions = 10,
+      stop.on.error = FALSE, .progress = FALSE
+    )
+  )
+  expect_silent(
+    predictDDConf_RT(
+      dd_params, maxrt = 1, subdivisions = 10,
+      minrt = dd_params$t0, scaled = FALSE, DistConf = NULL,
+      .progress = FALSE
+    )
+  )
+})
+
+
+test_that("predictConfModels and predictRTModels keep identifiers leading", {
+  params <- data.frame(
+    model = "dynWEV", sbj = 1,
+    a = 1, v = 0.4, t0 = 0.2, z = 0.5, tau = 0.1,
+    w = 0.5, svis = 0.1, sigvis = 0.05, theta1 = 0.6
+  )
+  conf_dist <- predictConfModels(
+    params, maxrt = 1, subdivisions = 5,
+    simult_conf = FALSE, stop.on.error = FALSE,
+    .progress = FALSE, parallel = FALSE
+  )
+  expect_identical(names(conf_dist)[1:2], c("sbj", "model"))
+
+  rt_dist <- predictRTModels(
+    params, maxrt = 1, subdivisions = 5, minrt = params$t0,
+    simult_conf = FALSE, scaled = FALSE, DistConf = NULL,
+    .progress = FALSE, parallel = FALSE
+  )
+  expect_identical(names(rt_dist)[1:2], c("sbj", "model"))
+})
+
 
 #as.matrix(jobs)[2,]