Handle optional parameters across predictors and simulators

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

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

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=""))

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

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) {

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)) {

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=""))

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

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

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
+}

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