diff --git a/.lintr b/.lintr
index 7049743a2..e1cd584ba 100644
--- a/.lintr
+++ b/.lintr
@@ -1,6 +1,7 @@
 linters: linters_with_defaults(
   line_length_linter = line_length_linter(120),
   cyclocomp_linter = NULL,
+  indentation_linter = NULL,
   object_usage_linter = NULL,
   object_length_linter = NULL,
   object_name_linter = object_name_linter(c("CamelCase", "camelCase", "snake_case"))
diff --git a/DESCRIPTION b/DESCRIPTION
index c2c85dc83..fcca6d997 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -110,6 +110,7 @@ Collate:
     'checkmate.R'
     'crmPack-package.R'
     'helpers_covr.R'
+    'helpers_design.R'
     'helpers_model.R'
     'logger.R'
     'utils-pipe.R'
diff --git a/NAMESPACE b/NAMESPACE
index 22a2fa797..b6e7c10b3 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -29,6 +29,7 @@ export(.DefaultCohortSizeParts)
 export(.DefaultCohortSizeRange)
 export(.DefaultDALogisticLogNormal)
 export(.DefaultDataGrouped)
+export(.DefaultDesignGrouped)
 export(.DefaultDualEndpoint)
 export(.DefaultDualEndpointBeta)
 export(.DefaultDualEndpointEmax)
@@ -87,6 +88,7 @@ export(.DefaultStoppingTargetBiomarker)
 export(.DefaultStoppingTargetProb)
 export(.DefaultTITELogisticLogNormal)
 export(.Design)
+export(.DesignGrouped)
 export(.DualDesign)
 export(.DualEndpoint)
 export(.DualEndpointBeta)
@@ -191,6 +193,7 @@ export(DataMixture)
 export(DataOrdinal)
 export(DataParts)
 export(Design)
+export(DesignGrouped)
 export(DualDesign)
 export(DualEndpoint)
 export(DualEndpointBeta)
@@ -378,6 +381,7 @@ exportClasses(DataMixture)
 exportClasses(DataOrdinal)
 exportClasses(DataParts)
 exportClasses(Design)
+exportClasses(DesignGrouped)
 exportClasses(DualDesign)
 exportClasses(DualEndpoint)
 exportClasses(DualEndpointBeta)
diff --git a/NEWS.md b/NEWS.md
index 25af45d59..67e80121a 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,5 +1,5 @@
 # Version 1.0.9000.9133
-* Added new `DataGrouped` class to support simultaneous dose escalation with monotherapy and combination therapy.
+* Added new `DataGrouped` and `DesignGrouped` classes with corresponding model `LogisticLogNormalGrouped` to support simultaneous dose escalation with monotherapy and combination therapy arms.
 * Created the `CrmPackClass` class as the ultimate ancestor of all other 
 `crmPack` classes to allow identification of crmPack classes and simpler 
 definition of generic methods.
diff --git a/R/Design-class.R b/R/Design-class.R
index a0687ccbc..05db7632b 100644
--- a/R/Design-class.R
+++ b/R/Design-class.R
@@ -529,3 +529,90 @@ DADesign <- function(model, data,
     safetyWindow = safetyWindow
   )
 }
+
+# DesignGrouped ----
+
+## class ----
+
+#' `DesignGrouped`
+#'
+#' @description `r lifecycle::badge("experimental")`
+#'
+#' [`DesignGrouped`] combines two [`Design`] objects: one for the mono and one
+#' for the combo arm of a joint dose escalation design.
+#'
+#' @slot model (`LogisticLogNormalGrouped`)\cr the model to be used, currently only one
+#'   class is allowed.
+#' @slot mono (`Design`)\cr defines the dose escalation rules for the mono arm, see
+#'   details.
+#' @slot combo (`Design`)\cr defines the dose escalation rules for the combo arm, see
+#'   details.
+#' @slot first_cohort_mono_only (`flag`)\cr whether first test one mono agent cohort, and then
+#'   once its DLT data has been collected, we proceed from the second cohort onwards with
+#'   concurrent mono and combo cohorts.
+#' @slot same_dose (`flag`)\cr whether the lower dose of the separately determined mono and combo
+#'   doses should be used as the next dose for both mono and combo.
+#'
+#' @details Note that the model slots inside the `mono` and `combo` parameters
+#'   are ignored (because we don't fit separate regression models for the mono and
+#'   combo arms). Instead, the `model` parameter is used to fit a joint regression
+#'   model for the mono and combo arms together.
+#'
+#' @aliases DesignGrouped
+#' @export
+#'
+.DesignGrouped <- setClass(
+  Class = "DesignGrouped",
+  slots = c(
+    model = "LogisticLogNormalGrouped",
+    mono = "Design",
+    combo = "Design",
+    first_cohort_mono_only = "logical",
+    same_dose = "logical"
+  ),
+  prototype = prototype(
+    model = .DefaultLogisticLogNormalGrouped(),
+    mono = .Design(),
+    combo = .Design(),
+    first_cohort_mono_only = TRUE,
+    same_dose = TRUE
+  ),
+  validity = v_design_grouped,
+  contains = "CrmPackClass"
+)
+
+## constructor ----
+
+#' @rdname DesignGrouped-class
+#'
+#' @param model (`LogisticLogNormalGrouped`)\cr see slot definition.
+#' @param mono (`Design`)\cr see slot definition.
+#' @param combo (`Design`)\cr see slot definition.
+#' @param first_cohort_mono_only (`flag`)\cr see slot definition.
+#' @param same_dose (`flag`)\cr see slot definition.
+#' @param ... not used.
+#'
+#' @export
+#' @example examples/Design-class-DesignGrouped.R
+#'
+DesignGrouped <- function(model,
+                          mono,
+                          combo = mono,
+                          first_cohort_mono_only = TRUE,
+                          same_dose = TRUE,
+                          ...) {
+  .DesignGrouped(
+    model = model,
+    mono = mono,
+    combo = combo,
+    first_cohort_mono_only = first_cohort_mono_only,
+    same_dose = same_dose
+  )
+}
+
+## default constructor ----
+
+#' @rdname DesignGrouped-class
+#' @note Typically, end-users will not use the `.DefaultDesignGrouped()` function.
+#' @export
+.DefaultDesignGrouped <- .DesignGrouped
diff --git a/R/Design-methods.R b/R/Design-methods.R
index 03e9b5be1..7d07790b1 100644
--- a/R/Design-methods.R
+++ b/R/Design-methods.R
@@ -7,147 +7,6 @@
 #' @include mcmc.R
 NULL
 
-# helper functions ----
-
-## set_seed ----
-
-#' Helper function to set and save the RNG seed
-#'
-#' @description `r lifecycle::badge("stable")`
-#'
-#' This code is basically copied from `stats:::simulate.lm`.
-#'
-#' @param seed an object specifying if and how the random number generator
-#' should be initialized ("seeded"). Either `NULL` (default) or an
-#' integer that will be used in a call to [set.seed()] before
-#' simulating the response vectors. If set, the value is saved as the
-#' `seed` slot of the returned object. The default, `NULL` will
-#' not change the random generator state.
-#' @return The integer vector containing the random number generate state will
-#' be returned, in order to call this function with this input to reproduce
-#' the obtained simulation results.
-#'
-#' @export
-#' @keywords programming
-set_seed <- function(seed = NULL) {
-  assert_number(seed, null.ok = TRUE)
-
-  if (!exists(".Random.seed", envir = .GlobalEnv, inherits = FALSE)) {
-    runif(1)
-  }
-
-  if (is.null(seed)) {
-    get(".Random.seed", envir = .GlobalEnv)
-  } else {
-    seed <- as.integer(seed)
-    r_seed <- get(".Random.seed", envir = .GlobalEnv)
-    # Make sure r_seed exists in parent frame.
-    assign(".r_seed", r_seed, envir = parent.frame())
-    set.seed(seed)
-    # Here we need the r_seed in the parent.frame!
-    do.call(
-      "on.exit",
-      list(quote(assign(".Random.seed", .r_seed, envir = .GlobalEnv))),
-      envir = parent.frame()
-    )
-    structure(seed, kind = as.list(RNGkind()))
-  }
-}
-
-## get_result_list ----
-
-#' Helper function to obtain simulation results list
-#'
-#' @description `r lifecycle::badge("stable")`
-#'
-#' The function `fun` can use variables that are visible to itself.
-#' The names of these variables have to be given in the vector `vars`.
-#'
-#' @param fun (`function`)\cr the simulation function for a single iteration, which takes as
-#' single parameter the iteration index.
-#' @param nsim number of simulations to be conducted.
-#' @param vars names of the variables.
-#' @param parallel should the simulation runs be parallelized across the
-#' clusters of the computer?
-#' @param n_cores how many cores should be used for parallel computing?
-#' @return The list with all simulation results (one iteration corresponds
-#' to one list element).
-#'
-#' @importFrom parallel makeCluster
-#' @importFrom parallelly availableCores
-#' @keywords internal programming
-get_result_list <- function(
-    fun,
-    nsim,
-    vars,
-    parallel,
-    n_cores) {
-  assert_flag(parallel)
-  assert_integerish(n_cores, lower = 1)
-
-  if (!parallel) {
-    lapply(
-      X = seq_len(nsim),
-      FUN = fun
-    )
-  } else {
-    # Process all simulations.
-    cores <- min(
-      safeInteger(n_cores),
-      parallelly::availableCores()
-    )
-
-    # Start the cluster.
-    cl <- parallel::makeCluster(cores)
-
-    # Load the required R package.
-    parallel::clusterEvalQ(cl, {
-      library(crmPack)
-      NULL
-    })
-
-    # Export local variables from the caller environment.
-    # Note: parent.frame() is different from parent.env() which returns
-    # the environment where this function has been defined!
-    parallel::clusterExport(
-      cl = cl,
-      varlist = vars,
-      envir = parent.frame()
-    )
-
-    # Export all global variables.
-    parallel::clusterExport(
-      cl = cl,
-      varlist = ls(.GlobalEnv)
-    )
-
-    # Load user extensions from global options.
-    crmpack_extensions <- getOption("crmpack_extensions")
-    if (is.null(crmpack_extensions) != TRUE) {
-      tryCatch(
-        {
-          parallel::clusterCall(cl, crmpack_extensions)
-        },
-        error = function(e) {
-          stop("Failed to export crmpack_extensions: ", e$message)
-        }
-      )
-    }
-
-    # Do the computations in parallel.
-    res <- parallel::parLapply(
-      cl = cl,
-      X = seq_len(nsim),
-      fun = fun
-    )
-
-    # Stop the cluster.
-    parallel::stopCluster(cl)
-
-    res
-  }
-}
-
 # nolint start
 
 ## ============================================================
@@ -4708,3 +4567,170 @@ setMethod("simulate",
 
 ## --------------------------------------------------------------------------
 # nolint end
+
+# simulate ----
+
+## DesignGrouped ----
+
+#' Simulate Method for the [`DesignGrouped`] Class
+#'
+#' @description `r lifecycle::badge("experimental")`
+#'
+#' A simulate method for [`DesignGrouped`] designs.
+#'
+#' @param object (`DesignGrouped`)\cr the design we want to simulate trials from.
+#' @param nsim (`number`)\cr how many trials should be simulated.
+#' @param seed (`RNGstate`)\cr generated with [set_seed()].
+#' @param truth (`function`)\cr a function which takes as input a dose (vector) and
+#'   returns the true probability (vector) for toxicity for the mono arm.
+#'   Additional arguments can be supplied in `args`.
+#' @param combo_truth (`function`)\cr same as `truth` but for the combo arm.
+#' @param args (`data.frame`)\cr optional `data.frame` with arguments that work
+#'   for both the `truth` and `combo_truth` functions. The column names correspond to
+#'   the argument names, the rows to the values of the arguments. The rows are
+#'   appropriately recycled in the `nsim` simulations.
+#' @param firstSeparate (`flag`)\cr whether to enroll the first patient separately
+#'   from the rest of the cohort and close the cohort in case a DLT occurs in this
+#'   first patient.
+#' @param mcmcOptions (`McmcOptions`)\cr MCMC options for each evaluation in the trial.
+#' @param parallel (`flag`)\cr whether the simulation runs are parallelized across the
+#'   cores of the computer.
+#' @param nCores (`number`)\cr how many cores should be used for parallel computing.
+#' @param ... not used.
+#'
+#' @return A list of `mono` and `combo` simulation results as [`Simulations`] objects.
+#'
+#' @aliases simulate-DesignGrouped
+#' @export
+#' @example examples/Design-method-simulate-DesignGrouped.R
+#'
+setMethod(
+  "simulate",
+  signature =
+    signature(
+      object = "DesignGrouped",
+      nsim = "ANY",
+      seed = "ANY"
+    ),
+  def =
+    function(object,
+             nsim = 1L,
+             seed = NULL,
+             truth,
+             combo_truth,
+             args = data.frame(),
+             firstSeparate = FALSE,
+             mcmcOptions = McmcOptions(),
+             parallel = FALSE,
+             nCores = min(parallelly::availableCores(), 5),
+             ...) {
+      nsim <- safeInteger(nsim)
+      assert_function(truth)
+      assert_function(combo_truth)
+      assert_data_frame(args)
+      assert_count(nsim, positive = TRUE)
+      assert_flag(firstSeparate)
+      assert_flag(parallel)
+      assert_count(nCores, positive = TRUE)
+
+      n_args <- max(nrow(args), 1L)
+      rng_state <- set_seed(seed)
+      sim_seeds <- sample.int(n = 2147483647, size = nsim)
+
+      run_sim <- function(iter_sim) {
+        set.seed(sim_seeds[iter_sim])
+        current <- list(mono = list(), combo = list())
+        # Define true toxicity functions.
+        current$args <- args[(iter_sim - 1) %% n_args + 1, , drop = FALSE]
+        current$mono$truth <- function(dose) do.call(truth, c(dose, current$args))
+        current$combo$truth <- function(dose) do.call(combo_truth, c(dose, current$args))
+        # Start the simulated data with the provided one.
+        current$mono$data <- object@mono@data
+        current$combo$data <- object@combo@data
+        # We are in the first cohort and continue for mono and combo.
+        current$first <- TRUE
+        current$mono$stop <- current$combo$stop <- FALSE
+        # What are the next doses to be used? Initialize with starting doses.
+        if (object@same_dose) {
+          current$mono$dose <- current$combo$dose <- min(object@mono@startingDose, object@combo@startingDose)
+        } else {
+          current$mono$dose <- object@mono@startingDose
+          current$combo$dose <- object@combo@startingDose
+        }
+        # Inside this loop we simulate the whole trial, until stopping.
+        while (!(current$mono$stop && current$combo$stop)) {
+          if (!current$mono$stop) {
+            current$mono$data <- current$mono$data |>
+              h_add_dlts(current$mono$dose, current$mono$truth, object@mono@cohort_size, firstSeparate)
+          }
+          if (!current$combo$stop && (!current$first || !object@first_cohort_mono_only)) {
+            current$combo$data <- current$combo$data |>
+              h_add_dlts(current$combo$dose, current$combo$truth, object@combo@cohort_size, firstSeparate)
+          }
+          if (current$first) current$first <- FALSE
+          current$grouped <- h_group_data(current$mono$data, current$combo$data)
+          current$samples <- mcmc(current$grouped, object@model, mcmcOptions)
+          if (!current$mono$stop) {
+            current$mono$limit <- maxDose(object@mono@increments, data = current$mono$data)
+            current$mono$dose <- object@mono@nextBest |>
+              nextBest(current$mono$limit, current$samples, object@model, current$grouped, group = "mono")
+            current$mono$dose <- current$mono$dose$value
+            current$mono$stop <- object@mono@stopping |>
+              stopTrial(current$mono$dose, current$samples, object@model, current$mono$data, group = "mono")
+            current$mono$results <- h_unpack_stopit(current$mono$stop)
+          }
+          if (!current$combo$stop) {
+            current$combo$limit <- if (is.na(current$mono$dose)) {
+              0
+            } else {
+              maxDose(object@combo@increments, current$combo$data) |>
+                min(current$mono$dose, na.rm = TRUE)
+            }
+            current$combo$dose <- object@combo@nextBest |>
+              nextBest(current$combo$limit, current$samples, object@model, current$grouped, group = "combo")
+            current$combo$dose <- current$combo$dose$value
+            current$combo$stop <- object@combo@stopping |>
+              stopTrial(current$combo$dose, current$samples, object@model, current$combo$data, group = "combo")
+            current$combo$results <- h_unpack_stopit(current$combo$stop)
+          }
+          if (object@same_dose && !current$mono$stop && !current$combo$stop) {
+            current$mono$dose <- current$combo$dose <- min(current$mono$dose, current$combo$dose)
+          }
+        }
+        current$mono$fit <- fit(current$samples, object@model, current$grouped, group = "mono")
+        current$combo$fit <- fit(current$samples, object@model, current$grouped, group = "combo")
+        lapply(
+          X = current[c("mono", "combo")], FUN = with,
+          list(
+            data = data, dose = dose, fit = subset(fit, select = -dose),
+            stop = attr(stop, "message"), results = results
+          )
+        )
+      }
+      vars_needed <- c("simSeeds", "args", "nArgs", "truth", "combo_truth", "firstSeparate", "object", "mcmcOptions")
+      result_list <- get_result_list(run_sim, nsim, vars_needed, parallel, nCores)
+      # Now we have a list with each element containing mono and combo. Reorder this a bit:
+      result_list <- list(
+        mono = lapply(result_list, "[[", "mono"),
+        combo = lapply(result_list, "[[", "combo")
+      )
+      # Put everything in a list with both mono and combo Simulations:
+      lapply(result_list, function(this_list) {
+        data_list <- lapply(this_list, "[[", "data")
+        recommended_doses <- as.numeric(sapply(this_list, "[[", "dose"))
+        fit_list <- lapply(this_list, "[[", "fit")
+        stop_reasons <- lapply(this_list, "[[", "stop")
+        report_results <- lapply(this_list, "[[", "results")
+        stop_report <- as.matrix(do.call(rbind, report_results))
+
+        Simulations(
+          data = data_list,
+          doses = recommended_doses,
+          fit = fit_list,
+          stop_reasons = stop_reasons,
+          stop_report = stop_report,
+          seed = rng_state
+        )
+      })
+    }
+)
diff --git a/R/Design-validity.R b/R/Design-validity.R
index d20b0f4ca..e5e53be4c 100644
--- a/R/Design-validity.R
+++ b/R/Design-validity.R
@@ -27,3 +27,19 @@ v_rule_design <- function(object) {
   )
   v$result()
 }
+
+
+#' @describeIn v_design validates that the [`DesignGrouped`] object
+#'   contains valid flags.
+v_design_grouped <- function(object) {
+  v <- Validate()
+  v$check(
+    test_flag(object@first_cohort_mono_only),
+    "first_cohort_mono_only must be a flag"
+  )
+  v$check(
+    test_flag(object@same_dose),
+    "same_dose must be a flag"
+  )
+  v$result()
+}
diff --git a/R/Rules-methods.R b/R/Rules-methods.R
index 8d5d32127..dd3af7367 100644
--- a/R/Rules-methods.R
+++ b/R/Rules-methods.R
@@ -67,7 +67,7 @@ setMethod(
   ),
   definition = function(nextBest, doselimit = Inf, samples, model, data, ...) {
     # Generate the MTD samples and derive the next best dose.
-    dose_target_samples <- dose(x = nextBest@target, model, samples)
+    dose_target_samples <- dose(x = nextBest@target, model, samples, ...)
     dose_target <- nextBest@derive(dose_target_samples)
 
     # Round to the next possible grid point.
@@ -292,7 +292,7 @@ setMethod("nextBest",
   ),
   definition = function(nextBest, doselimit = Inf, samples, model, data, ...) {
     # Matrix with samples from the dose-tox curve at the dose grid points.
-    prob_samples <- sapply(data@doseGrid, prob, model = model, samples = samples)
+    prob_samples <- sapply(data@doseGrid, prob, model = model, samples = samples, ...)
     # Compute probabilities to be in target and overdose tox interval.
     prob_underdosing <- colMeans(prob_samples < nextBest@target[1])
     prob_target <- colMeans(h_in_range(prob_samples, nextBest@target))
@@ -553,7 +553,7 @@ setMethod(
   ),
   definition = function(nextBest, doselimit = Inf, samples, model, data, ...) {
     # Matrix with samples from the dose-tox curve at the dose grid points.
-    prob_samples <- sapply(data@doseGrid, prob, model = model, samples = samples)
+    prob_samples <- sapply(data@doseGrid, prob, model = model, samples = samples, ...)
     dlt_prob <- colMeans(prob_samples)
 
     # Determine the dose with the closest distance.
@@ -651,7 +651,7 @@ setMethod(
   ),
   definition = function(nextBest, doselimit = Inf, samples, model, data, ...) {
     # Matrix with samples from the dose-tox curve at the dose grid points.
-    prob_samples <- sapply(data@doseGrid, prob, model = model, samples = samples)
+    prob_samples <- sapply(data@doseGrid, prob, model = model, samples = samples, ...)
 
     criterion <- colMeans(h_info_theory_dist(prob_samples, nextBest@target, nextBest@asymmetry))
 
@@ -697,8 +697,8 @@ setMethod(
 
     # Target dose estimates, i.e. the dose with probability of the occurrence of
     # a DLT that equals to the prob_target_drt or prob_target_eot.
-    dose_target_drt <- dose(x = prob_target_drt, model)
-    dose_target_eot <- dose(x = prob_target_eot, model)
+    dose_target_drt <- dose(x = prob_target_drt, model, ...)
+    dose_target_eot <- dose(x = prob_target_eot, model, ...)
 
     # Find the next best doses in the doseGrid. The next best dose is the dose
     # at level closest and below the target dose estimate.
@@ -732,7 +732,7 @@ setMethod(
       prob_target_eot = prob_target_eot,
       dose_target_eot = dose_target_eot,
       data = data,
-      prob_dlt = prob(dose = data@doseGrid, model = model),
+      prob_dlt = prob(dose = data@doseGrid, model = model, ...),
       doselimit = doselimit,
       next_dose = next_dose_drt
     )
@@ -777,12 +777,12 @@ setMethod(
     model = "LogisticIndepBeta",
     data = "Data"
   ),
-  definition = function(nextBest, doselimit = Inf, samples, model, data, ...) {
+  definition = function(nextBest, doselimit = Inf, samples, model, data, in_sim, ...) {
     # Generate target dose samples, i.e. the doses with probability of the
     # occurrence of a DLT that equals to the nextBest@prob_target_drt
     # (or nextBest@prob_target_eot, respectively).
-    dose_target_drt_samples <- dose(x = nextBest@prob_target_drt, model, samples)
-    dose_target_eot_samples <- dose(x = nextBest@prob_target_eot, model, samples)
+    dose_target_drt_samples <- dose(x = nextBest@prob_target_drt, model, samples, ...)
+    dose_target_eot_samples <- dose(x = nextBest@prob_target_eot, model, samples, ...)
 
     # Derive the prior/posterior estimates based on two above samples.
     dose_target_drt <- nextBest@derive(dose_target_drt_samples)
@@ -865,15 +865,15 @@ setMethod(
 
     # Target dose estimates, i.e. the dose with probability of the occurrence of
     # a DLT that equals to the prob_target_drt or prob_target_eot.
-    dose_target_drt <- dose(x = prob_target_drt, model)
-    dose_target_eot <- dose(x = prob_target_eot, model)
+    dose_target_drt <- dose(x = prob_target_drt, model, ...)
+    dose_target_eot <- dose(x = prob_target_eot, model, ...)
 
     # Find the dose which gives the maximum gain.
     dosegrid_range <- dose_grid_range(data)
     opt <- optim(
       par = dosegrid_range[1],
       fn = function(DOSE) {
-        -gain(DOSE, model_dle = model, model_eff = model_eff)
+        -gain(DOSE, model_dle = model, model_eff = model_eff, ...)
       },
       method = "L-BFGS-B",
       lower = dosegrid_range[1],
@@ -986,15 +986,15 @@ setMethod(
 
     # Generate target dose samples, i.e. the doses with probability of the
     # occurrence of a DLT that equals to the prob_target_drt or prob_target_eot.
-    dose_target_drt_samples <- dose(x = prob_target_drt, model, samples = samples)
-    dose_target_eot_samples <- dose(x = prob_target_eot, model, samples = samples)
+    dose_target_drt_samples <- dose(x = prob_target_drt, model, samples = samples, ...)
+    dose_target_eot_samples <- dose(x = prob_target_eot, model, samples = samples, ...)
 
     # Derive the prior/posterior estimates based on two above samples.
     dose_target_drt <- nextBest@derive(dose_target_drt_samples)
     dose_target_eot <- nextBest@derive(dose_target_eot_samples)
 
     # Gain samples.
-    gain_samples <- sapply(data@doseGrid, gain, model, samples, model_eff, samples_eff)
+    gain_samples <- sapply(data@doseGrid, gain, model, samples, model_eff, samples_eff, ...)
     # For every sample, get the dose (from the dose grid) that gives the maximum gain value.
     dose_lev_mg_samples <- apply(gain_samples, 1, which.max)
     dose_mg_samples <- data@doseGrid[dose_lev_mg_samples]
@@ -1085,7 +1085,7 @@ setMethod(
   ),
   definition = function(nextBest, doselimit, samples, model, data, ...) {
     # Matrix with samples from the dose-tox curve at the dose grid points.
-    prob_samples <- sapply(data@doseGrid, prob, model = model, samples = samples)
+    prob_samples <- sapply(data@doseGrid, prob, model = model, samples = samples, ...)
 
     # Determine the maximum dose level with a toxicity probability below or
     # equal to the target and calculate how often a dose is selected as MTD
@@ -1237,7 +1237,7 @@ setMethod(
   ),
   definition = function(nextBest, doselimit, samples, model, data, ...) {
     # Matrix with samples from the dose-tox curve at the dose grid points.
-    prob_samples <- sapply(data@doseGrid, prob, model = model, samples = samples)
+    prob_samples <- sapply(data@doseGrid, prob, model = model, samples = samples, ...)
 
     # Determine which dose level has the minimum distance to target.
     dose_min_mtd_dist <- apply(
@@ -2284,7 +2284,8 @@ setMethod("stopTrial",
       mtdSamples <- dose(
         x = stopping@target,
         model,
-        samples
+        samples,
+        ...
       )
 
       ## what is the absolute threshold?
@@ -2354,7 +2355,8 @@ setMethod(
     mtd_samples <- dose(
       x = stopping@target,
       model,
-      samples
+      samples,
+      ...
     )
     # CV of MTD expressed as percentage, derived based on MTD posterior samples.
     mtd_cv <- (mad(mtd_samples) / median(mtd_samples)) * 100
@@ -2459,7 +2461,7 @@ setMethod(
   definition = function(stopping, dose, samples, model, data, ...) {
     # Compute the target biomarker prob at this dose.
     # Get the biomarker level samples at the dose grid points.
-    biom_level_samples <- biomarker(xLevel = seq_len(data@nGrid), model, samples)
+    biom_level_samples <- biomarker(xLevel = seq_len(data@nGrid), model, samples, ...)
 
     # If target is relative to maximum.
     if (stopping@is_relative) {
@@ -2910,7 +2912,8 @@ setMethod(
     dose_target_samples <- dose(
       x = stopping@prob_target,
       model = model,
-      samples = samples
+      samples = samples,
+      ...
     )
     # 95% credibility interval.
     dose_target_ci <- quantile(dose_target_samples, probs = c(0.025, 0.975))
@@ -2955,7 +2958,7 @@ setMethod("stopTrial",
       assert_probability(stopping@prob_target)
 
       prob_target <- stopping@prob_target
-      dose_target_samples <- dose(x = prob_target, model = model)
+      dose_target_samples <- dose(x = prob_target, model = model, ...)
       ## Find the variance of the log of the dose_target_samples(eta)
       M1 <- matrix(c(-1 / (model@phi2), -(log(prob_target / (1 - prob_target)) - model@phi1) / (model@phi2)^2), 1, 2)
       M2 <- model@Pcov
@@ -3024,7 +3027,8 @@ setMethod("stopTrial",
       TDtargetEndOfTrialSamples <- dose(
         x = prob_target,
         model = model,
-        samples = samples
+        samples = samples,
+        ...
       )
       ## Find the TDtarget End of trial estimate
       TDtargetEndOfTrialEstimate <- TDderive(TDtargetEndOfTrialSamples)
@@ -3047,7 +3051,8 @@ setMethod("stopTrial",
           model,
           samples,
           Effmodel,
-          Effsamples
+          Effsamples,
+          ...
         )
       }
 
@@ -3137,12 +3142,13 @@ setMethod("stopTrial",
       ## find the TDtarget End of Trial
       TDtargetEndOfTrial <- dose(
         x = prob_target,
-        model = model
+        model = model,
+        ...
       )
 
       ## Find the dose with maximum gain value
       Gainfun <- function(DOSE) {
-        -gain(DOSE, model_dle = model, model_eff = Effmodel)
+        -gain(DOSE, model_dle = model, model_eff = Effmodel, ...)
       }
 
       # if(data@placebo) {
diff --git a/R/crmPack-package.R b/R/crmPack-package.R
index 44fbf1ee0..a07ec3797 100644
--- a/R/crmPack-package.R
+++ b/R/crmPack-package.R
@@ -138,7 +138,8 @@ globalVariables(c(
   "comp",
   "X",
   "skel_probs",
-  "is_combo"
+  "is_combo",
+  "results"
 ))
 
 # nolint end
diff --git a/R/helpers_design.R b/R/helpers_design.R
new file mode 100644
index 000000000..b30026b78
--- /dev/null
+++ b/R/helpers_design.R
@@ -0,0 +1,168 @@
+#' Helper Function to Set and Save the RNG Seed
+#'
+#' @description `r lifecycle::badge("stable")`
+#'
+#' This code is basically copied from `stats:::simulate.lm`.
+#'
+#' @param seed an object specifying if and how the random number generator
+#' should be initialized ("seeded"). Either `NULL` (default) or an
+#' integer that will be used in a call to [set.seed()] before
+#' simulating the response vectors. If set, the value is saved as the
+#' `seed` slot of the returned object. The default, `NULL` will
+#' not change the random generator state.
+#' @return The integer vector containing the random number generate state will
+#' be returned, in order to call this function with this input to reproduce
+#' the obtained simulation results.
+#'
+#' @export
+set_seed <- function(seed = NULL) {
+  assert_number(seed, null.ok = TRUE)
+
+  if (!exists(".Random.seed", envir = .GlobalEnv, inherits = FALSE)) {
+    runif(1)
+  }
+
+  if (is.null(seed)) {
+    get(".Random.seed", envir = .GlobalEnv)
+  } else {
+    seed <- as.integer(seed)
+    r_seed <- get(".Random.seed", envir = .GlobalEnv)
+    # Make sure r_seed exists in parent frame.
+    assign(".r_seed", r_seed, envir = parent.frame())
+    set.seed(seed)
+    # Here we need the r_seed in the parent.frame!
+    do.call(
+      "on.exit",
+      list(quote(assign(".Random.seed", .r_seed, envir = .GlobalEnv))),
+      envir = parent.frame()
+    )
+    structure(seed, kind = as.list(RNGkind()))
+  }
+}
+
+#' Helper Function to Obtain Simulation Results List
+#'
+#' The function `fun` can use variables that are visible to itself.
+#' The names of these variables have to be given in the vector `vars`.
+#'
+#' @param fun (`function`)\cr the simulation function for a single iteration, which takes as
+#' single parameter the iteration index.
+#' @param nsim number of simulations to be conducted.
+#' @param vars names of the variables.
+#' @param parallel should the simulation runs be parallelized across the
+#' clusters of the computer?
+#' @param n_cores how many cores should be used for parallel computing?
+#' @return The list with all simulation results (one iteration corresponds
+#' to one list element).
+#'
+#' @importFrom parallel makeCluster
+#' @importFrom parallelly availableCores
+#' @keywords internal programming
+get_result_list <- function(
+    fun,
+    nsim,
+    vars,
+    parallel,
+    n_cores) {
+  assert_flag(parallel)
+  assert_integerish(n_cores, lower = 1)
+
+  if (!parallel) {
+    lapply(
+      X = seq_len(nsim),
+      FUN = fun
+    )
+  } else {
+    # Process all simulations.
+    cores <- min(
+      safeInteger(n_cores),
+      parallelly::availableCores()
+    )
+
+    # Start the cluster.
+    cl <- parallel::makeCluster(cores)
+
+    # Load the required R package.
+    parallel::clusterEvalQ(cl, {
+      library(crmPack)
+      NULL
+    })
+
+    # Export local variables from the caller environment.
+    # Note: parent.frame() is different from parent.env() which returns
+    # the environment where this function has been defined!
+    parallel::clusterExport(
+      cl = cl,
+      varlist = vars,
+      envir = parent.frame()
+    )
+
+    # Export all global variables.
+    parallel::clusterExport(
+      cl = cl,
+      varlist = ls(.GlobalEnv)
+    )
+
+    # Load user extensions from global options.
+    crmpack_extensions <- getOption("crmpack_extensions")
+    if (is.null(crmpack_extensions) != TRUE) {
+      tryCatch(
+        {
+          parallel::clusterCall(cl, crmpack_extensions)
+        },
+        error = function(e) {
+          stop("Failed to export crmpack_extensions: ", e$message)
+        }
+      )
+    }
+
+    # Do the computations in parallel.
+    res <- parallel::parLapply(
+      cl = cl,
+      X = seq_len(nsim),
+      fun = fun
+    )
+
+    # Stop the cluster.
+    parallel::stopCluster(cl)
+
+    res
+  }
+}
+
+#' Helper Function to Add Randomly Generated DLTs During Simulations
+#'
+#' @param data (`Data`)\cr what data to start from.
+#' @param dose (`number`)\cr current dose.
+#' @param truth (`function`)\cr defines the true probability for a DLT at a dose.
+#' @param cohort_size (`CohortSize`)\cr the cohort size rule to use.
+#' @param first_separate (`flag`)\cr whether the first patient is enrolled separately.
+#'
+#' @return The updated `data`.
+#'
+#' @keywords internal
+h_add_dlts <- function(data,
+                       dose,
+                       truth,
+                       cohort_size,
+                       first_separate) {
+  assert_class(data, "Data")
+  assert_number(dose)
+  assert_function(truth)
+  assert_class(cohort_size, "CohortSize")
+  assert_flag(first_separate)
+
+  prob <- truth(dose)
+  size <- size(cohort_size, dose = dose, data = data)
+  dlts <- if (first_separate && size > 1) {
+    first_dlts <- rbinom(n = 1, size = 1, prob = prob)
+    if (first_dlts == 0) {
+      c(first_dlts, rbinom(n = size - 1, size = 1, prob = prob))
+    } else {
+      first_dlts
+    }
+  } else {
+    rbinom(n = size, size = 1, prob = prob)
+  }
+  update(data, x = dose, y = dlts)
+}
diff --git a/_pkgdown.yaml b/_pkgdown.yaml
index a8f2f343b..33dc7e0f6 100644
--- a/_pkgdown.yaml
+++ b/_pkgdown.yaml
@@ -129,6 +129,7 @@ reference:
       - DualDesign
       - TDsamplesDesign
       - TDDesign
+      - DesignGrouped
   - title: Internal Helper Functions
     contents:
       - h_blind_plot_data
@@ -453,6 +454,7 @@ reference:
       - simulate,RuleDesign-method
       - simulate,TDDesign-method
       - simulate,TDsamplesDesign-method
+      - simulate-DesignGrouped
       - summary,DualSimulations-method
       - summary,GeneralSimulations-method
       - summary,PseudoDualFlexiSimulations-method
diff --git a/examples/Design-class-DesignGrouped.R b/examples/Design-class-DesignGrouped.R
new file mode 100644
index 000000000..d116306ca
--- /dev/null
+++ b/examples/Design-class-DesignGrouped.R
@@ -0,0 +1,55 @@
+empty_data <- Data(doseGrid = c(1, 3, 5, 10, 15, 20, 25, 40, 50, 80, 100))
+
+# Initialize the joint model.
+my_model <- LogisticLogNormalGrouped(
+  mean = c(-0.85, 0, 1, 0),
+  cov = diag(1, 4),
+  ref_dose = 56
+)
+
+# Choose the rule for selecting the next dose.
+my_next_best <- NextBestNCRM(
+  target = c(0.2, 0.35),
+  overdose = c(0.35, 1),
+  max_overdose_prob = 0.25
+)
+
+# Choose the rule for the cohort-size.
+my_size1 <- CohortSizeRange(
+  intervals = c(0, 30),
+  cohort_size = c(1, 3)
+)
+my_size2 <- CohortSizeDLT(
+  intervals = c(0, 1),
+  cohort_size = c(1, 3)
+)
+my_size <- maxSize(my_size1, my_size2)
+
+# Choose the rule for stopping.
+my_stopping1 <- StoppingMinCohorts(nCohorts = 3)
+my_stopping2 <- StoppingTargetProb(
+  target = c(0.2, 0.35),
+  prob = 0.5
+)
+my_stopping3 <- StoppingMinPatients(nPatients = 20)
+my_stopping <- (my_stopping1 & my_stopping2) | my_stopping3
+
+# Choose the rule for dose increments.
+my_increments <- IncrementsRelative(
+  intervals = c(0, 20),
+  increments = c(1, 0.33)
+)
+
+# Initialize the design.
+design <- DesignGrouped(
+  model = my_model,
+  mono = Design(
+    model = .DefaultModelLogNormal(), # Ignored.
+    nextBest = my_next_best,
+    stopping = my_stopping,
+    increments = my_increments,
+    cohort_size = my_size,
+    data = empty_data,
+    startingDose = 3
+  )
+)
diff --git a/examples/Design-method-simulate-DesignGrouped.R b/examples/Design-method-simulate-DesignGrouped.R
new file mode 100644
index 000000000..f37b725f3
--- /dev/null
+++ b/examples/Design-method-simulate-DesignGrouped.R
@@ -0,0 +1,74 @@
+# Assemble ingredients for our group design.
+my_stopping <- StoppingTargetProb(target = c(0.2, 0.35), prob = 0.5) |
+  StoppingMinPatients(20) |
+  StoppingMissingDose()
+my_increments <- IncrementsDoseLevels(levels = 3L)
+my_next_best <- NextBestNCRM(
+  target = c(0.2, 0.3),
+  overdose = c(0.3, 1),
+  max_overdose_prob = 0.3
+)
+my_cohort_size <- CohortSizeConst(3)
+empty_data <- Data(doseGrid = c(0.1, 0.5, 1.5, 3, 6, seq(from = 10, to = 80, by = 2)))
+my_model <- LogisticLogNormalGrouped(
+  mean = c(-4, -4, -4, -4),
+  cov = diag(rep(6, 4)),
+  ref_dose = 0.1
+)
+
+# Put together the design. Note that if we only specify the mono arm,
+# then the combo arm is having the same settings.
+my_design <- DesignGrouped(
+  model = my_model,
+  mono = Design(
+    model = my_model,
+    stopping = my_stopping,
+    increments = my_increments,
+    nextBest = my_next_best,
+    cohort_size = my_cohort_size,
+    data = empty_data,
+    startingDose = 0.1
+  ),
+  first_cohort_mono_only = TRUE,
+  same_dose = TRUE
+)
+
+# Set up a realistic simulation scenario.
+my_truth <- function(x) plogis(-4 + 0.2 * log(x / 0.1))
+my_combo_truth <- function(x) plogis(-4 + 0.5 * log(x / 0.1))
+matplot(
+  x = empty_data@doseGrid,
+  y = cbind(
+    mono = my_truth(empty_data@doseGrid),
+    combo = my_combo_truth(empty_data@doseGrid)
+  ),
+  type = "l",
+  ylab = "true DLT prob",
+  xlab = "dose"
+)
+legend("topright", c("mono", "combo"), lty = c(1, 2), col = c(1, 2))
+
+# Start the simulations.
+set.seed(123)
+my_sims <- simulate(
+  my_design,
+  nsim = 4, # This should be at least 100 in actual applications.
+  seed = 123,
+  truth = my_truth,
+  combo_truth = my_combo_truth
+)
+
+# Looking at the summary of the simulations:
+mono_sims_sum <- summary(my_sims$mono, truth = my_truth)
+combo_sims_sum <- summary(my_sims$combo, truth = my_combo_truth)
+
+mono_sims_sum
+combo_sims_sum
+
+plot(mono_sims_sum)
+plot(combo_sims_sum)
+
+# Looking at specific simulated trials:
+trial_index <- 1
+plot(my_sims$mono@data[[trial_index]])
+plot(my_sims$combo@data[[trial_index]])
diff --git a/man/DesignGrouped-class.Rd b/man/DesignGrouped-class.Rd
new file mode 100644
index 000000000..d44ef56eb
--- /dev/null
+++ b/man/DesignGrouped-class.Rd
@@ -0,0 +1,124 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/Design-class.R
+\docType{class}
+\name{DesignGrouped-class}
+\alias{DesignGrouped-class}
+\alias{.DesignGrouped}
+\alias{DesignGrouped}
+\alias{.DefaultDesignGrouped}
+\title{\code{DesignGrouped}}
+\usage{
+DesignGrouped(
+  model,
+  mono,
+  combo = mono,
+  first_cohort_mono_only = TRUE,
+  same_dose = TRUE,
+  ...
+)
+}
+\arguments{
+\item{model}{(\code{LogisticLogNormalGrouped})\cr see slot definition.}
+
+\item{mono}{(\code{Design})\cr see slot definition.}
+
+\item{combo}{(\code{Design})\cr see slot definition.}
+
+\item{first_cohort_mono_only}{(\code{flag})\cr see slot definition.}
+
+\item{same_dose}{(\code{flag})\cr see slot definition.}
+
+\item{...}{not used.}
+}
+\description{
+\ifelse{html}{\href{https://lifecycle.r-lib.org/articles/stages.html#experimental}{\figure{lifecycle-experimental.svg}{options: alt='[Experimental]'}}}{\strong{[Experimental]}}
+
+\code{\link{DesignGrouped}} combines two \code{\link{Design}} objects: one for the mono and one
+for the combo arm of a joint dose escalation design.
+}
+\details{
+Note that the model slots inside the \code{mono} and \code{combo} parameters
+are ignored (because we don't fit separate regression models for the mono and
+combo arms). Instead, the \code{model} parameter is used to fit a joint regression
+model for the mono and combo arms together.
+}
+\section{Slots}{
+
+\describe{
+\item{\code{model}}{(\code{LogisticLogNormalGrouped})\cr the model to be used, currently only one
+class is allowed.}
+
+\item{\code{mono}}{(\code{Design})\cr defines the dose escalation rules for the mono arm, see
+details.}
+
+\item{\code{combo}}{(\code{Design})\cr defines the dose escalation rules for the combo arm, see
+details.}
+
+\item{\code{first_cohort_mono_only}}{(\code{flag})\cr whether first test one mono agent cohort, and then
+once its DLT data has been collected, we proceed from the second cohort onwards with
+concurrent mono and combo cohorts.}
+
+\item{\code{same_dose}}{(\code{flag})\cr whether the lower dose of the separately determined mono and combo
+doses should be used as the next dose for both mono and combo.}
+}}
+
+\note{
+Typically, end-users will not use the \code{.DefaultDesignGrouped()} function.
+}
+\examples{
+empty_data <- Data(doseGrid = c(1, 3, 5, 10, 15, 20, 25, 40, 50, 80, 100))
+
+# Initialize the joint model.
+my_model <- LogisticLogNormalGrouped(
+  mean = c(-0.85, 0, 1, 0),
+  cov = diag(1, 4),
+  ref_dose = 56
+)
+
+# Choose the rule for selecting the next dose.
+my_next_best <- NextBestNCRM(
+  target = c(0.2, 0.35),
+  overdose = c(0.35, 1),
+  max_overdose_prob = 0.25
+)
+
+# Choose the rule for the cohort-size.
+my_size1 <- CohortSizeRange(
+  intervals = c(0, 30),
+  cohort_size = c(1, 3)
+)
+my_size2 <- CohortSizeDLT(
+  intervals = c(0, 1),
+  cohort_size = c(1, 3)
+)
+my_size <- maxSize(my_size1, my_size2)
+
+# Choose the rule for stopping.
+my_stopping1 <- StoppingMinCohorts(nCohorts = 3)
+my_stopping2 <- StoppingTargetProb(
+  target = c(0.2, 0.35),
+  prob = 0.5
+)
+my_stopping3 <- StoppingMinPatients(nPatients = 20)
+my_stopping <- (my_stopping1 & my_stopping2) | my_stopping3
+
+# Choose the rule for dose increments.
+my_increments <- IncrementsRelative(
+  intervals = c(0, 20),
+  increments = c(1, 0.33)
+)
+
+# Initialize the design.
+design <- DesignGrouped(
+  model = my_model,
+  mono = Design(
+    model = .DefaultModelLogNormal(), # Ignored.
+    nextBest = my_next_best,
+    stopping = my_stopping,
+    increments = my_increments,
+    cohort_size = my_size,
+    data = empty_data,
+    startingDose = 3
+  )
+)
+}
diff --git a/man/get_result_list.Rd b/man/get_result_list.Rd
index da89b4edc..e47267ea0 100644
--- a/man/get_result_list.Rd
+++ b/man/get_result_list.Rd
@@ -1,8 +1,8 @@
 % Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/Design-methods.R
+% Please edit documentation in R/helpers_design.R
 \name{get_result_list}
 \alias{get_result_list}
-\title{Helper function to obtain simulation results list}
+\title{Helper Function to Obtain Simulation Results List}
 \usage{
 get_result_list(fun, nsim, vars, parallel, n_cores)
 }
@@ -24,8 +24,6 @@ The list with all simulation results (one iteration corresponds
 to one list element).
 }
 \description{
-\ifelse{html}{\href{https://lifecycle.r-lib.org/articles/stages.html#stable}{\figure{lifecycle-stable.svg}{options: alt='[Stable]'}}}{\strong{[Stable]}}
-
 The function \code{fun} can use variables that are visible to itself.
 The names of these variables have to be given in the vector \code{vars}.
 }
diff --git a/man/h_add_dlts.Rd b/man/h_add_dlts.Rd
new file mode 100644
index 000000000..a7809024e
--- /dev/null
+++ b/man/h_add_dlts.Rd
@@ -0,0 +1,26 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/helpers_design.R
+\name{h_add_dlts}
+\alias{h_add_dlts}
+\title{Helper Function to Add Randomly Generated DLTs During Simulations}
+\usage{
+h_add_dlts(data, dose, truth, cohort_size, first_separate)
+}
+\arguments{
+\item{data}{(\code{Data})\cr what data to start from.}
+
+\item{dose}{(\code{number})\cr current dose.}
+
+\item{truth}{(\code{function})\cr defines the true probability for a DLT at a dose.}
+
+\item{cohort_size}{(\code{CohortSize})\cr the cohort size rule to use.}
+
+\item{first_separate}{(\code{flag})\cr whether the first patient is enrolled separately.}
+}
+\value{
+The updated \code{data}.
+}
+\description{
+Helper Function to Add Randomly Generated DLTs During Simulations
+}
+\keyword{internal}
diff --git a/man/nextBest.Rd b/man/nextBest.Rd
index fd438dbfc..b9b730b3d 100644
--- a/man/nextBest.Rd
+++ b/man/nextBest.Rd
@@ -52,7 +52,7 @@ nextBest(nextBest, doselimit, samples, model, data, ...)
 
 \S4method{nextBest}{NextBestTD,numeric,missing,LogisticIndepBeta,Data}(nextBest, doselimit = Inf, model, data, in_sim = FALSE, ...)
 
-\S4method{nextBest}{NextBestTDsamples,numeric,Samples,LogisticIndepBeta,Data}(nextBest, doselimit = Inf, samples, model, data, ...)
+\S4method{nextBest}{NextBestTDsamples,numeric,Samples,LogisticIndepBeta,Data}(nextBest, doselimit = Inf, samples, model, data, in_sim, ...)
 
 \S4method{nextBest}{NextBestMaxGain,numeric,missing,ModelTox,DataDual}(
   nextBest,
diff --git a/man/set_seed.Rd b/man/set_seed.Rd
index f63be09de..982d23bad 100644
--- a/man/set_seed.Rd
+++ b/man/set_seed.Rd
@@ -1,8 +1,8 @@
 % Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/Design-methods.R
+% Please edit documentation in R/helpers_design.R
 \name{set_seed}
 \alias{set_seed}
-\title{Helper function to set and save the RNG seed}
+\title{Helper Function to Set and Save the RNG Seed}
 \usage{
 set_seed(seed = NULL)
 }
@@ -24,4 +24,3 @@ the obtained simulation results.
 
 This code is basically copied from \code{stats:::simulate.lm}.
 }
-\keyword{programming}
diff --git a/man/simulate-DesignGrouped-method.Rd b/man/simulate-DesignGrouped-method.Rd
new file mode 100644
index 000000000..54c2c4880
--- /dev/null
+++ b/man/simulate-DesignGrouped-method.Rd
@@ -0,0 +1,136 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/Design-methods.R
+\name{simulate,DesignGrouped-method}
+\alias{simulate,DesignGrouped-method}
+\alias{simulate-DesignGrouped}
+\title{Simulate Method for the \code{\link{DesignGrouped}} Class}
+\usage{
+\S4method{simulate}{DesignGrouped}(
+  object,
+  nsim = 1L,
+  seed = NULL,
+  truth,
+  combo_truth,
+  args = data.frame(),
+  firstSeparate = FALSE,
+  mcmcOptions = McmcOptions(),
+  parallel = FALSE,
+  nCores = min(parallelly::availableCores(), 5),
+  ...
+)
+}
+\arguments{
+\item{object}{(\code{DesignGrouped})\cr the design we want to simulate trials from.}
+
+\item{nsim}{(\code{number})\cr how many trials should be simulated.}
+
+\item{seed}{(\code{RNGstate})\cr generated with \code{\link[=set_seed]{set_seed()}}.}
+
+\item{truth}{(\code{function})\cr a function which takes as input a dose (vector) and
+returns the true probability (vector) for toxicity for the mono arm.
+Additional arguments can be supplied in \code{args}.}
+
+\item{combo_truth}{(\code{function})\cr same as \code{truth} but for the combo arm.}
+
+\item{args}{(\code{data.frame})\cr optional \code{data.frame} with arguments that work
+for both the \code{truth} and \code{combo_truth} functions. The column names correspond to
+the argument names, the rows to the values of the arguments. The rows are
+appropriately recycled in the \code{nsim} simulations.}
+
+\item{firstSeparate}{(\code{flag})\cr whether to enroll the first patient separately
+from the rest of the cohort and close the cohort in case a DLT occurs in this
+first patient.}
+
+\item{mcmcOptions}{(\code{McmcOptions})\cr MCMC options for each evaluation in the trial.}
+
+\item{parallel}{(\code{flag})\cr whether the simulation runs are parallelized across the
+cores of the computer.}
+
+\item{nCores}{(\code{number})\cr how many cores should be used for parallel computing.}
+
+\item{...}{not used.}
+}
+\value{
+A list of \code{mono} and \code{combo} simulation results as \code{\link{Simulations}} objects.
+}
+\description{
+\ifelse{html}{\href{https://lifecycle.r-lib.org/articles/stages.html#experimental}{\figure{lifecycle-experimental.svg}{options: alt='[Experimental]'}}}{\strong{[Experimental]}}
+
+A simulate method for \code{\link{DesignGrouped}} designs.
+}
+\examples{
+# Assemble ingredients for our group design.
+my_stopping <- StoppingTargetProb(target = c(0.2, 0.35), prob = 0.5) |
+  StoppingMinPatients(20) |
+  StoppingMissingDose()
+my_increments <- IncrementsDoseLevels(levels = 3L)
+my_next_best <- NextBestNCRM(
+  target = c(0.2, 0.3),
+  overdose = c(0.3, 1),
+  max_overdose_prob = 0.3
+)
+my_cohort_size <- CohortSizeConst(3)
+empty_data <- Data(doseGrid = c(0.1, 0.5, 1.5, 3, 6, seq(from = 10, to = 80, by = 2)))
+my_model <- LogisticLogNormalGrouped(
+  mean = c(-4, -4, -4, -4),
+  cov = diag(rep(6, 4)),
+  ref_dose = 0.1
+)
+
+# Put together the design. Note that if we only specify the mono arm,
+# then the combo arm is having the same settings.
+my_design <- DesignGrouped(
+  model = my_model,
+  mono = Design(
+    model = my_model,
+    stopping = my_stopping,
+    increments = my_increments,
+    nextBest = my_next_best,
+    cohort_size = my_cohort_size,
+    data = empty_data,
+    startingDose = 0.1
+  ),
+  first_cohort_mono_only = TRUE,
+  same_dose = TRUE
+)
+
+# Set up a realistic simulation scenario.
+my_truth <- function(x) plogis(-4 + 0.2 * log(x / 0.1))
+my_combo_truth <- function(x) plogis(-4 + 0.5 * log(x / 0.1))
+matplot(
+  x = empty_data@doseGrid,
+  y = cbind(
+    mono = my_truth(empty_data@doseGrid),
+    combo = my_combo_truth(empty_data@doseGrid)
+  ),
+  type = "l",
+  ylab = "true DLT prob",
+  xlab = "dose"
+)
+legend("topright", c("mono", "combo"), lty = c(1, 2), col = c(1, 2))
+
+# Start the simulations.
+set.seed(123)
+my_sims <- simulate(
+  my_design,
+  nsim = 4, # This should be at least 100 in actual applications.
+  seed = 123,
+  truth = my_truth,
+  combo_truth = my_combo_truth
+)
+
+# Looking at the summary of the simulations:
+mono_sims_sum <- summary(my_sims$mono, truth = my_truth)
+combo_sims_sum <- summary(my_sims$combo, truth = my_combo_truth)
+
+mono_sims_sum
+combo_sims_sum
+
+plot(mono_sims_sum)
+plot(combo_sims_sum)
+
+# Looking at specific simulated trials:
+trial_index <- 1
+plot(my_sims$mono@data[[trial_index]])
+plot(my_sims$combo@data[[trial_index]])
+}
diff --git a/man/v_design.Rd b/man/v_design.Rd
index c77e39a5c..8f00076d7 100644
--- a/man/v_design.Rd
+++ b/man/v_design.Rd
@@ -3,9 +3,12 @@
 \name{v_design}
 \alias{v_design}
 \alias{v_rule_design}
+\alias{v_design_grouped}
 \title{Internal Helper Functions for Validation of \code{\link{RuleDesign}} Objects}
 \usage{
 v_rule_design(object)
+
+v_design_grouped(object)
 }
 \arguments{
 \item{object}{(\code{RuleDesign})\cr object to validate.}
@@ -25,4 +28,7 @@ These functions are only used internally to validate the format of an input
 \item \code{v_rule_design()}: validates that the \code{\link{RuleDesign}} object
 contains valid \code{startingDose}.
 
+\item \code{v_design_grouped()}: validates that the \code{\link{DesignGrouped}} object
+contains valid flags.
+
 }}
diff --git a/tests/testthat/helper-model.R b/tests/testthat/helper-model.R
index e1725a900..b56e23665 100644
--- a/tests/testthat/helper-model.R
+++ b/tests/testthat/helper-model.R
@@ -386,22 +386,3 @@ h_get_fractional_crm <- function() {
     sigma2 = 2
   )
 }
-
-.NeedsExtraProbModel <- setClass("NeedsExtraProbModel", contains = "LogisticKadane")
-setMethod(
-  f = "prob",
-  signature = signature(
-    dose = "numeric",
-    model = "NeedsExtraProbModel",
-    samples = "Samples"
-  ),
-  definition = function(dose, model, samples, extra_argument) {
-    if (missing(extra_argument)) stop("we need extra_argument")
-    # We don't forward to LogisticKadane method here since that would
-    # not accept the extra argument.
-    rep(0.5, size(samples))
-  }
-)
-h_needs_extra_prob_model <- function() {
-  .NeedsExtraProbModel(h_get_logistic_kadane())
-}
diff --git a/tests/testthat/test-Design-class.R b/tests/testthat/test-Design-class.R
index adf469a4b..ea21ead52 100644
--- a/tests/testthat/test-Design-class.R
+++ b/tests/testthat/test-Design-class.R
@@ -424,3 +424,47 @@ test_that("DADesign user constructor arguments names are as expected", {
     ordered = TRUE
   )
 })
+
+# DesignGrouped ----
+
+test_that(".DesignGrouped works as expected", {
+  result <- .DesignGrouped()
+
+  expect_true(inherits(result, "CrmPackClass"))
+  expect_valid(result, "DesignGrouped")
+})
+
+test_that("DesignGrouped works as expected", {
+  empty_data <- Data(doseGrid = 2:50)
+  model <- .DefaultLogisticLogNormalGrouped()
+  stopping <- h_stopping_target_prob()
+  increments <- h_increments_relative()
+  placebo_cohort_size <- CohortSizeConst(0L)
+  next_best <- h_next_best_ncrm()
+  cohort_size <- CohortSizeRange(intervals = c(0, 30), cohort_size = c(1, 3))
+
+  result <- expect_silent(
+    DesignGrouped(
+      model = model,
+      mono = Design(
+        model,
+        stopping,
+        increments,
+        nextBest = next_best,
+        cohort_size = cohort_size,
+        data = empty_data,
+        startingDose = 3
+      )
+    )
+  )
+
+  expect_valid(result, "DesignGrouped")
+  expect_identical(result@mono, result@combo)
+  expect_true(result@first_cohort_mono_only)
+  expect_true(result@same_dose)
+})
+
+test_that(".DefaultDesignGrouped works as expected", {
+  result <- .DefaultDesignGrouped()
+  expect_valid(result, "DesignGrouped")
+})
diff --git a/tests/testthat/test-Design-methods.R b/tests/testthat/test-Design-methods.R
index 773875b4f..f5a7c01eb 100644
--- a/tests/testthat/test-Design-methods.R
+++ b/tests/testthat/test-Design-methods.R
@@ -1,51 +1,3 @@
-# helper functions ----
-
-## set_seed ----
-
-test_that("set_seed returns correct value if seed is a value", {
-  seed <- 1.909
-  seed_int <- 1
-
-  RNGkind("default")
-  rng_state <- set_seed(seed)
-  attr(seed_int, "kind") <- list("Mersenne-Twister", "Inversion", "Rejection")
-  expect_equal(rng_state, seed_int)
-
-  RNGkind("Super-Duper")
-  rng_state <- set_seed(seed)
-  attr(seed_int, "kind") <- list("Super-Duper", "Inversion", "Rejection")
-  expect_equal(rng_state, seed_int)
-
-  RNGkind("default")
-})
-
-test_that("set_seed returns correct value if seed is NULL", {
-  seed <- NULL
-
-  RNGkind("default")
-  rng_state <- set_seed(seed)
-  expect_equal(rng_state, .Random.seed)
-
-  RNGkind("Super-Duper")
-  rng_state <- set_seed(seed)
-  expect_equal(rng_state, .Random.seed)
-
-  RNGkind("default")
-})
-
-## get_result_list ----
-
-test_that("get_result_list returns correct value", {
-  res <- get_result_list(mean, 2, NULL, FALSE, 5)
-  expect_equal(res, list(1, 2))
-
-  res <- get_result_list(length, 2, NULL, FALSE, 5)
-  expect_equal(res, list(1, 1))
-
-  expect_error(get_result_list(length, 2, NULL, 5, 5))
-  expect_error(get_result_list(length, 2, NULL, FALSE, 0))
-})
-
 # simulate ----
 
 ## NextBestInfTheory ----
@@ -182,6 +134,151 @@ test_that("stop_reasons can be NA with certain stopping rule settings", {
   expect_identical(result, expected)
 })
 
+## DesignGrouped ----
+
+test_that("simulate for DesignGrouped works as expected", {
+  object <- DesignGrouped(
+    model = LogisticLogNormalGrouped(mean = rep(-1, 4), cov = diag(5, 4), ref_dose = 1),
+    mono = Design(
+      model = .LogisticNormal(),
+      nextBest = NextBestNCRM(target = c(0.3, 0.6), overdose = c(0.6, 1), max_overdose_prob = 0.7),
+      stopping = StoppingMinPatients(nPatients = 9),
+      increments = IncrementsDoseLevels(levels = 5),
+      cohort_size = CohortSizeConst(3),
+      data = Data(doseGrid = 1:100),
+      startingDose = 1
+    ),
+    same_dose = TRUE,
+    first_cohort_mono_only = TRUE
+  )
+  my_truth <- function(x) plogis(-4 + 0.2 * log(x / 0.1))
+  my_combo_truth <- function(x) plogis(-4 + 0.5 * log(x / 0.1))
+
+  result <- expect_silent(simulate(
+    object,
+    nsim = 2,
+    seed = 123,
+    truth = my_truth,
+    combo_truth = my_combo_truth,
+    mcmcOptions = h_get_mcmc_options()
+  ))
+
+  expect_list(result)
+  expect_names(names(result), identical.to = c("mono", "combo"))
+  expect_valid(result$mono, "Simulations")
+  expect_valid(result$combo, "Simulations")
+
+  mono_trial <- result$mono@data[[2L]]
+  combo_trial <- result$combo@data[[2L]]
+
+  # First cohort is only mono at the starting dose (lowest in dose grid).
+  expect_true(all(mono_trial@xLevel[1:3] == 1))
+
+  # We have the same dose for subsequent cohorts.
+  expect_true(all(mono_trial@xLevel[4:6] == combo_trial@xLevel[1:3]))
+  expect_true(all(mono_trial@xLevel[7:9] == combo_trial@xLevel[4:6]))
+})
+
+test_that("simulate for DesignGrouped works as expected with different doses, parallel first cohort", {
+  object <- DesignGrouped(
+    model = LogisticLogNormalGrouped(mean = rep(-1, 4), cov = diag(5, 4), ref_dose = 1),
+    mono = Design(
+      model = .LogisticNormal(),
+      nextBest = NextBestNCRM(target = c(0.3, 0.6), overdose = c(0.6, 1), max_overdose_prob = 0.7),
+      stopping = StoppingMinPatients(nPatients = 20),
+      increments = IncrementsDoseLevels(levels = 5),
+      cohort_size = CohortSizeConst(3),
+      data = Data(doseGrid = 1:100),
+      startingDose = 1
+    ),
+    same_dose = FALSE,
+    first_cohort_mono_only = FALSE
+  )
+  my_truth <- function(x) plogis(-4 + 0.2 * log(x / 0.1))
+  my_combo_truth <- function(x) plogis(-4 + 0.9 * log(x / 0.1))
+
+  result <- expect_silent(simulate(
+    object,
+    nsim = 1,
+    seed = 123,
+    truth = my_truth,
+    combo_truth = my_combo_truth,
+    mcmcOptions = h_get_mcmc_options()
+  ))
+
+  expect_list(result)
+  expect_names(names(result), identical.to = c("mono", "combo"))
+  expect_valid(result$mono, "Simulations")
+  expect_valid(result$combo, "Simulations")
+
+  mono_trial <- result$mono@data[[1L]]
+  combo_trial <- result$combo@data[[1L]]
+
+  # First cohort is joint at starting dose.
+  expect_true(all(mono_trial@xLevel[1:3] == combo_trial@xLevel[1:3]))
+
+  # We have different doses in subsequent cohorts.
+  expect_false(all(mono_trial@xLevel[4:20] == combo_trial@xLevel[4:20]))
+})
+
+test_that("simulate for DesignGrouped works when first patient is dosed separately, different combo design", {
+  object <- DesignGrouped(
+    model = LogisticLogNormalGrouped(mean = rep(-1, 4), cov = diag(5, 4), ref_dose = 1),
+    mono = Design(
+      model = .LogisticNormal(),
+      nextBest = NextBestNCRM(target = c(0.3, 0.6), overdose = c(0.6, 1), max_overdose_prob = 0.7),
+      stopping = StoppingMinPatients(nPatients = 10),
+      increments = IncrementsDoseLevels(levels = 3),
+      cohort_size = CohortSizeConst(2),
+      data = Data(doseGrid = 1:100),
+      startingDose = 10
+    ),
+    combo = Design(
+      model = .LogisticNormal(),
+      nextBest = NextBestNCRM(target = c(0.3, 0.6), overdose = c(0.6, 1), max_overdose_prob = 0.7),
+      stopping = StoppingMinPatients(nPatients = 20),
+      increments = IncrementsDoseLevels(levels = 5),
+      cohort_size = CohortSizeConst(3),
+      data = Data(doseGrid = 1:100),
+      startingDose = 1
+    ),
+    same_dose = FALSE,
+    first_cohort_mono_only = FALSE
+  )
+  my_truth <- function(x) plogis(-4 + 0.2 * log(x / 0.1))
+  my_combo_truth <- function(x) plogis(-2 + 0.9 * log(x / 0.1))
+
+  result <- expect_silent(simulate(
+    object,
+    nsim = 1,
+    seed = 123,
+    truth = my_truth,
+    firstSeparate = TRUE,
+    combo_truth = my_combo_truth,
+    mcmcOptions = h_get_mcmc_options()
+  ))
+
+  expect_list(result)
+  expect_names(names(result), identical.to = c("mono", "combo"))
+  expect_valid(result$mono, "Simulations")
+  expect_valid(result$combo, "Simulations")
+
+  mono_trial <- result$mono@data[[1L]]
+  combo_trial <- result$combo@data[[1L]]
+
+  # We expect at least one cohort with just one patient in the combo arm here
+  # because of the high toxicity.
+  expect_true(any(table(combo_trial@cohort) == 1))
+
+  # Check that we had the different cohort sizes between the two arms.
+  expect_true(max(table(mono_trial@cohort)) == 2)
+  expect_true(max(table(combo_trial@cohort)) == 3)
+
+  # Check that we had different starting doses in the two arms.
+  expect_true(mono_trial@x[1] == 10)
+  expect_true(combo_trial@x[1] == 1)
+})
+
 # examine ----
 
 ## DADesign ----
diff --git a/tests/testthat/test-Design-validity.R b/tests/testthat/test-Design-validity.R
index 5374ad254..651a1e3ba 100644
--- a/tests/testthat/test-Design-validity.R
+++ b/tests/testthat/test-Design-validity.R
@@ -61,3 +61,24 @@ test_that("v_rule_design returns message when startingDose is not on doseGrid",
   object@startingDose <- 6.5
   expect_equal(v_rule_design(object), err_msg)
 })
+
+## v_design_grouped ----
+
+test_that("v_design_grouped passes for valid object", {
+  object <- .DesignGrouped()
+  expect_true(v_design_grouped(object))
+})
+
+test_that("v_design_grouped messages wrong flag slots as expected", {
+  object <- .DesignGrouped()
+
+  object@same_dose <- c(NA, TRUE)
+  object@first_cohort_mono_only <- logical()
+
+  result <- v_design_grouped(object)
+  expected <- c(
+    "first_cohort_mono_only must be a flag",
+    "same_dose must be a flag"
+  )
+  expect_identical(result, expected)
+})
diff --git a/tests/testthat/test-Rules-methods.R b/tests/testthat/test-Rules-methods.R
index 890dd0f06..57d504c4c 100644
--- a/tests/testthat/test-Rules-methods.R
+++ b/tests/testthat/test-Rules-methods.R
@@ -81,13 +81,13 @@ test_that("nextBest-NextBestNCRM returns expected values of the objects (no dose
 })
 
 test_that("nextBest-NextBestNCRM can accept additional arguments and pass them to prob inside", {
-  my_data <- h_get_data()
-  my_model <- h_needs_extra_prob_model()
+  my_data <- h_get_data_grouped()
+  my_model <- .DefaultLogisticLogNormalGrouped()
   my_samples <- mcmc(my_data, my_model, h_get_mcmc_options(samples = 10, burnin = 10))
   nb_ncrm <- NextBestNCRM(
     target = c(0.2, 0.35), overdose = c(0.35, 1), max_overdose_prob = 0.25
   )
-  result <- nextBest(nb_ncrm, Inf, my_samples, my_model, my_data, extra_argument = "foo")
+  result <- nextBest(nb_ncrm, Inf, my_samples, my_model, my_data, group = "mono")
   expect_identical(result$value, NA_real_)
 })
 
@@ -1857,8 +1857,8 @@ test_that("StoppingTargetProb works correctly when above threshold", {
 })
 
 test_that("stopTrial-StoppingTargetProb can accept additional arguments and pass them to prob", {
-  my_data <- h_get_data()
-  my_model <- h_needs_extra_prob_model()
+  my_data <- h_get_data_grouped()
+  my_model <- .DefaultLogisticLogNormalGrouped()
   my_samples <- mcmc(my_data, my_model, h_get_mcmc_options(samples = 10, burnin = 10))
   stopping <- StoppingTargetProb(target = c(0.1, 0.4), prob = 0.3)
   result <- stopTrial(
@@ -1867,7 +1867,7 @@ test_that("stopTrial-StoppingTargetProb can accept additional arguments and pass
     samples = my_samples,
     model = my_model,
     data = my_data,
-    extra_argument = "bla"
+    group = "combo"
   )
   expect_false(result)
 })
diff --git a/tests/testthat/test-helpers.R b/tests/testthat/test-helpers.R
index ad687b3ba..90e0be660 100644
--- a/tests/testthat/test-helpers.R
+++ b/tests/testthat/test-helpers.R
@@ -473,7 +473,7 @@ test_that("h_find_interval works as expected for custom replacement", {
 test_that("default constructors exist for all subclasses of GeneralModel", {
   allModelSubclasses <- names(getClassDef("GeneralModel")@subclasses)
   # Exceptions.
-  classesNotToTest <- c("DualEndpoint", "NeedsExtraProbModel")
+  classesNotToTest <- "DualEndpoint"
   classesToTest <- setdiff(allModelSubclasses, classesNotToTest)
   lapply(
     classesToTest,
diff --git a/tests/testthat/test-helpers_design.R b/tests/testthat/test-helpers_design.R
new file mode 100644
index 000000000..2c8c18641
--- /dev/null
+++ b/tests/testthat/test-helpers_design.R
@@ -0,0 +1,100 @@
+# set_seed ----
+
+test_that("set_seed returns correct value if seed is a value", {
+  seed <- 1.909
+  seed_int <- 1
+
+  RNGkind("default")
+  rng_state <- set_seed(seed)
+  attr(seed_int, "kind") <- list("Mersenne-Twister", "Inversion", "Rejection")
+  expect_equal(rng_state, seed_int)
+
+  RNGkind("Super-Duper")
+  rng_state <- set_seed(seed)
+  attr(seed_int, "kind") <- list("Super-Duper", "Inversion", "Rejection")
+  expect_equal(rng_state, seed_int)
+
+  RNGkind("default")
+})
+
+test_that("set_seed returns correct value if seed is NULL", {
+  seed <- NULL
+
+  RNGkind("default")
+  rng_state <- set_seed(seed)
+  expect_equal(rng_state, .Random.seed)
+
+  RNGkind("Super-Duper")
+  rng_state <- set_seed(seed)
+  expect_equal(rng_state, .Random.seed)
+
+  RNGkind("default")
+})
+
+# get_result_list ----
+
+test_that("get_result_list returns correct value", {
+  res <- get_result_list(mean, 2, NULL, FALSE, 5)
+  expect_equal(res, list(1, 2))
+
+  res <- get_result_list(length, 2, NULL, FALSE, 5)
+  expect_equal(res, list(1, 1))
+
+  expect_error(get_result_list(length, 2, NULL, 5, 5))
+  expect_error(get_result_list(length, 2, NULL, FALSE, 0))
+})
+
+# h_add_dlts ----
+
+test_that("h_add_dlts works as expected", {
+  data <- h_get_data()
+  cohort_size <- CohortSizeConst(3)
+
+  set.seed(123)
+  result <- expect_silent(h_add_dlts(
+    data = data,
+    dose = data@doseGrid[3],
+    truth = plogis,
+    cohort_size = cohort_size,
+    first_separate = FALSE
+  ))
+  expect_valid(result, "Data")
+  expect_equal(tail(result@x, 3), rep(data@doseGrid[3], 3))
+  expect_true(data@nObs + 3 == result@nObs)
+})
+
+test_that("h_add_dlts works as expected when first separate patient has a DLT", {
+  data <- h_get_data()
+  cohort_size <- CohortSizeConst(3)
+
+  set.seed(123)
+  result <- expect_silent(h_add_dlts(
+    data = data,
+    dose = data@doseGrid[3],
+    truth = function(dose, ...) 1, # Make sure the first patient has a DLT.
+    cohort_size = cohort_size,
+    first_separate = TRUE
+  ))
+  expect_valid(result, "Data")
+  expect_true(tail(result@y, 1) == 1)
+  expect_equal(tail(result@x, 1), data@doseGrid[3])
+  expect_true(data@nObs + 1 == result@nObs)
+})
+
+test_that("h_add_dlts works as expected when first separate patient does not have a DLT", {
+  data <- h_get_data()
+  cohort_size <- CohortSizeConst(3)
+
+  set.seed(123)
+  result <- expect_silent(h_add_dlts(
+    data = data,
+    dose = data@doseGrid[3],
+    truth = function(dose, ...) 0, # Make sure the first patient does not have a DLT.
+    cohort_size = cohort_size,
+    first_separate = TRUE
+  ))
+  expect_valid(result, "Data")
+  expect_equal(tail(result@y, 3), rep(0, 3))
+  expect_equal(tail(result@x, 3), rep(data@doseGrid[3], 3))
+  expect_true(data@nObs + 3 == result@nObs)
+})