Skip to contents
# Core simulation framework (Timer, Population, Trial, deterministic_schedule, add_timepoints, ...)
library(rxsim)

# Analyses
library(multcomp)     # Dunnett
#> Loading required package: mvtnorm
#> Loading required package: survival
#> Loading required package: TH.data
#> Loading required package: MASS
#> 
#> Attaching package: 'TH.data'
#> The following object is masked from 'package:MASS':
#> 
#>     geyser
library(DoseFinding)  # MCP-Mod

# Helpers
library(dplyr)
#> 
#> Attaching package: 'dplyr'
#> The following object is masked from 'package:MASS':
#> 
#>     select
#> The following objects are masked from 'package:stats':
#> 
#>     filter, lag
#> The following objects are masked from 'package:base':
#> 
#>     intersect, setdiff, setequal, union
set.seed(4566)

Dose-finding trials aim to identify the dose-response relationship and establish an effective dose level. This example simulates a multi-arm fixed design across a placebo and four active doses, with the endpoint generated from an Emax dose-response model. At the final analysis, two methods are applied: a Dunnett test (comparing each active dose against placebo while controlling family-wise error rate) and MCP-Mod (a model-based contrast approach that selects the best-fitting dose-response shape from a candidate set).

Unique focus: multi-arm design, Emax dose-response data generation, Dunnett test vs MCP-Mod - comparing detection rates across effect sizes.

The simulation skeleton is the same as Example 1: define the scenario, create population generators, register conditions, and run replicate trials. This vignette focuses on what changes in a dose-finding setting.

Scenario

A multi-arm, fixed design (placebo + 4 active doses) with a single continuous endpoint. At the final analysis, we perform:

  • Dunnett test: all active doses vs placebo using a normal-theory linear model.
  • MCP-Mod: model-based multiple contrast test + model fitting on a candidate set of dose–response shapes.

The Emax model generates mean responses via E(d) = e0 + emax × d / (ed50 + d), where e0 = 0 (placebo baseline), emax = 1 (maximum achievable effect), and ed50 = 20 (dose at half-maximum effect). arm_names = paste0("d", doses) creates human-readable labels (d0, d5, d10, d20, d50) that propagate through trial data and analysis results. delta = 0.1 is the minimum effect size of clinical relevance used by MCP-Mod, and alpha = 0.05 controls the family-wise Type I error rate for both testing procedures.

# Dose levels (placebo + 4 actives)
doses <- c(0, 5, 10, 20, 50)

# Total N and allocation (balanced across arms)
sample_size <- 150
allocation  <- rep(1, length(doses))
arm_names   <- paste0("d", doses)  # e.g., d0, d5, ... used as arm labels

# Data-generating model: Emax with homoscedastic noise
e0    <- 0.0
emax  <- 1.0
ed50  <- 20.0
sigma <- 1.0

mean_fun <- function(d) e0 + emax * d / (ed50 + d)

# Operating characteristics
alpha <- 0.05
delta <- 0.1

scenario <- tidyr::expand_grid(
  sample_size = sample_size,
  allocation  = list(allocation),
  n_arms      = length(doses),
  alpha       = alpha,
  delta       = delta
)

enrollment_fn <- function(n) rexp(n, rate = 1)
dropout_fn    <- function(n) rexp(n, rate = 0.01)

Populations

mk_pop_gen is a closure factory: it captures the dose value d and returns a generator function that, when called with n, draws n responses from N(E(d), sigma^2). Each arm’s generator stores dose as a numeric column because MCP-Mod requires actual dose values - not just arm labels - to fit dose-response models and evaluate contrasts at analysis time.

mk_pop_gen <- function(d) {
  function(n) {
    mu <- mean_fun(d)
    y  <- rnorm(n, mean = mu, sd = sigma)
    data.frame(
      id = seq_len(n),
      dose = d,
      y = y,
      readout_time = 1
    )
  }
}

population_generators <- lapply(seq_along(doses), function(i) {
  mk_pop_gen(doses[i])
})
names(population_generators) <- arm_names

Conditions

At the final time, run:

  1. Dunnett test (active vs placebo) using multcomp::glht on lm(y ~ arm).

  2. MCP-Mod using DoseFinding::MCPMod with a candidate model set (Mods).

The Dunnett test (multcomp::glht with mcp(arm = "Dunnett")) simultaneously compares each active dose against placebo while controlling the family-wise error rate at alpha. MCP-Mod is run with five candidate dose-response shapes; selModel = "aveAIC" selects the best model by AIC-weighted averaging, and Delta = delta sets the minimum effect size of clinical relevance for the contrast step. mct_min_p extracts the minimum p-value across all candidate model contrast tests - a small value indicates that at least one candidate model detects a dose-response signal.

# Candidate model set for MCP-Mod
mods <- Mods(
  linear     = NULL,
  emax       = 20,
  exponential= 50,
  sigEmax    = c(20, 3),
  quadratic  = -0.2,
  doses      = doses
)

analysis_generators <- list(
  final = list(
    trigger = enroll_trigger(1.0, sample_size),
    analysis = function(df, current_time){
      df_e <- df |>
        dplyr::filter(!is.na(enroll_time)) |>
        dplyr::mutate(
          arm  = factor(arm, levels = paste0("d", doses)),
          dose = as.numeric(dose)
        )

      # 1) Dunnett (active vs placebo)
      fit <- lm(y ~ arm, data = df_e)
      dun  <- multcomp::glht(fit, linfct = multcomp::mcp(arm = "Dunnett"))
      summ <- summary(dun)

      # 2) MCP-Mod (one-step)
      mm <- DoseFinding::MCPMod(
          dose   = df_e$dose,
          resp   = df_e$y,
          models = mods,
          type   = "normal",
          Delta  = delta,
          alpha  = alpha,
          selModel = "aveAIC"
      )

      mct_min_p <- min(attr(mm$MCTtest$tStat, "pVal"), na.rm = TRUE)

      data.frame(
        scenario,
        n_total      = nrow(df_e),
        dunn_min_p   = summ$test$pvalues |> as.numeric() |> min(),
        mcpmod_min_p = mct_min_p,
        stringsAsFactors = FALSE
      )
    }
  )
)

Simulate

set.seed(5)
trials <- replicate_trial(
  trial_name            = "multiarm_dunnett_mcpmod",
  sample_size           = sample_size,
  arms                  = arm_names,
  allocation            = allocation,
  enrollment            = enrollment_fn,
  dropout               = dropout_fn,
  analysis_generators   = analysis_generators,
  population_generators = population_generators,
  n                     = 3
)
run_trials(trials)
#> Warning in MCTpval(contMat, corMat, df, tStat, alternative, mvtcontrol):
#> Warning from mvtnorm::pmvt: Completion with error > abseps.

Results

collect_results() row-binds analysis outputs across all replicates and prepends replicate, timepoint, and analysis columns. dunn_min_p is the smallest Dunnett-adjusted p-value across the four active-dose comparisons - a value below alpha indicates that at least one dose significantly differs from placebo after multiplicity correction. mcpmod_min_p is the minimum MCP-Mod contrast test p-value across candidate models; MCP-Mod is generally more powerful when the true dose-response shape is well-captured by one of the candidate models.

replicate_results <- collect_results(trials)
replicate_results
#>   replicate timepoint analysis sample_size    allocation n_arms alpha delta
#> 1         1  163.6457    final         150 1, 1, 1, 1, 1      5  0.05   0.1
#> 2         2  154.0418    final         150 1, 1, 1, 1, 1      5  0.05   0.1
#> 3         3  148.2880    final         150 1, 1, 1, 1, 1      5  0.05   0.1
#>   n_total   dunn_min_p mcpmod_min_p
#> 1     150 1.255875e-01 3.571107e-02
#> 2     150 2.087216e-06 3.870043e-08
#> 3     150 1.129478e-03 2.088616e-05

Power curve

How does detection rate scale with the true maximum effect (emax)? We sweep emax over four values, keeping ed50 = 20 and n = 150 fixed, and compare the Dunnett and MCP-Mod rejection rates.

set.seed(55)
n_reps_pw <- 100
emax_vals <- c(0.5, 1.0, 1.5, 2.0)

pw_df5 <- do.call(rbind, lapply(emax_vals, function(em) {
  mf <- function(d) 0 + em * d / (20 + d)
  pop_pw <- lapply(seq_along(doses), function(i) {
    d <- doses[i]
    local({
      mu_d <- mf(d)
      function(n) {
        data.frame(
          id = 1:n, dose = d,
          y  = rnorm(n, mu_d, sigma),
          readout_time = 1
        )
      }
    })
  })
  names(pop_pw) <- arm_names

  an_pw <- list(final = list(
    trigger  = enroll_trigger(1.0, sample_size),
    analysis = function(df, ct) {
      df_e <- dplyr::filter(df, !is.na(enroll_time)) |>
        dplyr::mutate(
          arm  = factor(arm, levels = arm_names),
          dose = as.numeric(dose)
        )
      fit  <- lm(y ~ arm, data = df_e)
      dun  <- multcomp::glht(fit, linfct = multcomp::mcp(arm = "Dunnett"))
      mm   <- DoseFinding::MCPMod(
        dose = df_e$dose, resp = df_e$y,
        models = mods, type = "normal",
        Delta = delta, alpha = alpha, selModel = "aveAIC"
      )
      data.frame(
        dunn_sig   = as.integer(min(summary(dun)$test$pvalues) < alpha),
        mcpmod_sig = as.integer(
          min(attr(mm$MCTtest$tStat, "pVal"), na.rm = TRUE) < alpha
        )
      )
    }
  ))
  tr <- replicate_trial(
    "pw5", sample_size, arm_names, allocation,
    enrollment_fn, dropout_fn, an_pw, pop_pw, n_reps_pw
  )
  invisible(run_trials(tr))
  res <- collect_results(tr)
  data.frame(
    emax          = em,
    power_dunnett = mean(res$dunn_sig),
    power_mcpmod  = mean(res$mcpmod_sig)
  )
}))
#> Warning in MCTpval(contMat, corMat, df, tStat, alternative, mvtcontrol):
#> Warning from mvtnorm::pmvt: Completion with error > abseps.
#> Warning in MCTpval(contMat, corMat, df, tStat, alternative, mvtcontrol):
#> Warning from mvtnorm::pmvt: Completion with error > abseps.
#> Warning in MCTpval(contMat, corMat, df, tStat, alternative, mvtcontrol):
#> Warning from mvtnorm::pmvt: Completion with error > abseps.
#> Warning in MCTpval(contMat, corMat, df, tStat, alternative, mvtcontrol):
#> Warning from mvtnorm::pmvt: Completion with error > abseps.
#> Warning in MCTpval(contMat, corMat, df, tStat, alternative, mvtcontrol):
#> Warning from mvtnorm::pmvt: Completion with error > abseps.
#> Warning in MCTpval(contMat, corMat, df, tStat, alternative, mvtcontrol):
#> Warning from mvtnorm::pmvt: Completion with error > abseps.
#> Warning in MCTpval(contMat, corMat, df, tStat, alternative, mvtcontrol):
#> Warning from mvtnorm::pmvt: Completion with error > abseps.
#> Warning in MCTpval(contMat, corMat, df, tStat, alternative, mvtcontrol):
#> Warning from mvtnorm::pmvt: Completion with error > abseps.
#> Warning in MCTpval(contMat, corMat, df, tStat, alternative, mvtcontrol):
#> Warning from mvtnorm::pmvt: Completion with error > abseps.
#> Warning in MCTpval(contMat, corMat, df, tStat, alternative, mvtcontrol):
#> Warning from mvtnorm::pmvt: Completion with error > abseps.
#> Warning in MCTpval(contMat, corMat, df, tStat, alternative, mvtcontrol):
#> Warning from mvtnorm::pmvt: Completion with error > abseps.
#> Warning in MCTpval(contMat, corMat, df, tStat, alternative, mvtcontrol):
#> Warning from mvtnorm::pmvt: Completion with error > abseps.
#> Warning in MCTpval(contMat, corMat, df, tStat, alternative, mvtcontrol):
#> Warning from mvtnorm::pmvt: Completion with error > abseps.
#> Warning in MCTpval(contMat, corMat, df, tStat, alternative, mvtcontrol):
#> Warning from mvtnorm::pmvt: Completion with error > abseps.
#> Warning in MCTpval(contMat, corMat, df, tStat, alternative, mvtcontrol):
#> Warning from mvtnorm::pmvt: Completion with error > abseps.
#> Warning in MCTpval(contMat, corMat, df, tStat, alternative, mvtcontrol):
#> Warning from mvtnorm::pmvt: Completion with error > abseps.
#> Warning in MCTpval(contMat, corMat, df, tStat, alternative, mvtcontrol):
#> Warning from mvtnorm::pmvt: Completion with error > abseps.
#> Warning in MCTpval(contMat, corMat, df, tStat, alternative, mvtcontrol):
#> Warning from mvtnorm::pmvt: Completion with error > abseps.
#> Warning in MCTpval(contMat, corMat, df, tStat, alternative, mvtcontrol):
#> Warning from mvtnorm::pmvt: Completion with error > abseps.
#> Warning in MCTpval(contMat, corMat, df, tStat, alternative, mvtcontrol):
#> Warning from mvtnorm::pmvt: Completion with error > abseps.
#> Warning in MCTpval(contMat, corMat, df, tStat, alternative, mvtcontrol):
#> Warning from mvtnorm::pmvt: Completion with error > abseps.
#> Warning in MCTpval(contMat, corMat, df, tStat, alternative, mvtcontrol):
#> Warning from mvtnorm::pmvt: Completion with error > abseps.
#> Warning in MCTpval(contMat, corMat, df, tStat, alternative, mvtcontrol):
#> Warning from mvtnorm::pmvt: Completion with error > abseps.
#> Warning in MCTpval(contMat, corMat, df, tStat, alternative, mvtcontrol):
#> Warning from mvtnorm::pmvt: Completion with error > abseps.
#> Warning in MCTpval(contMat, corMat, df, tStat, alternative, mvtcontrol):
#> Warning from mvtnorm::pmvt: Completion with error > abseps.
#> Warning in MCTpval(contMat, corMat, df, tStat, alternative, mvtcontrol):
#> Warning from mvtnorm::pmvt: Completion with error > abseps.
#> Warning in MCTpval(contMat, corMat, df, tStat, alternative, mvtcontrol):
#> Warning from mvtnorm::pmvt: Completion with error > abseps.
#> Warning in MCTpval(contMat, corMat, df, tStat, alternative, mvtcontrol):
#> Warning from mvtnorm::pmvt: Completion with error > abseps.
#> Warning in MCTpval(contMat, corMat, df, tStat, alternative, mvtcontrol):
#> Warning from mvtnorm::pmvt: Completion with error > abseps.
#> Warning in MCTpval(contMat, corMat, df, tStat, alternative, mvtcontrol):
#> Warning from mvtnorm::pmvt: Completion with error > abseps.
#> Warning in MCTpval(contMat, corMat, df, tStat, alternative, mvtcontrol):
#> Warning from mvtnorm::pmvt: Completion with error > abseps.
#> Warning in MCTpval(contMat, corMat, df, tStat, alternative, mvtcontrol):
#> Warning from mvtnorm::pmvt: Completion with error > abseps.
#> Warning in MCTpval(contMat, corMat, df, tStat, alternative, mvtcontrol):
#> Warning from mvtnorm::pmvt: Completion with error > abseps.
#> Warning in MCTpval(contMat, corMat, df, tStat, alternative, mvtcontrol):
#> Warning from mvtnorm::pmvt: Completion with error > abseps.
#> Warning in MCTpval(contMat, corMat, df, tStat, alternative, mvtcontrol):
#> Warning from mvtnorm::pmvt: Completion with error > abseps.
#> Warning in MCTpval(contMat, corMat, df, tStat, alternative, mvtcontrol):
#> Warning from mvtnorm::pmvt: Completion with error > abseps.
#> Warning in MCTpval(contMat, corMat, df, tStat, alternative, mvtcontrol):
#> Warning from mvtnorm::pmvt: Completion with error > abseps.
#> Warning in MCTpval(contMat, corMat, df, tStat, alternative, mvtcontrol):
#> Warning from mvtnorm::pmvt: Completion with error > abseps.
#> Warning in MCTpval(contMat, corMat, df, tStat, alternative, mvtcontrol):
#> Warning from mvtnorm::pmvt: Completion with error > abseps.
#> Warning in MCTpval(contMat, corMat, df, tStat, alternative, mvtcontrol):
#> Warning from mvtnorm::pmvt: Completion with error > abseps.
#> Warning in MCTpval(contMat, corMat, df, tStat, alternative, mvtcontrol):
#> Warning from mvtnorm::pmvt: Completion with error > abseps.
#> Warning in MCTpval(contMat, corMat, df, tStat, alternative, mvtcontrol):
#> Warning from mvtnorm::pmvt: Completion with error > abseps.
#> Warning in MCTpval(contMat, corMat, df, tStat, alternative, mvtcontrol):
#> Warning from mvtnorm::pmvt: Completion with error > abseps.
#> Warning in MCTpval(contMat, corMat, df, tStat, alternative, mvtcontrol):
#> Warning from mvtnorm::pmvt: Completion with error > abseps.
#> Warning in MCTpval(contMat, corMat, df, tStat, alternative, mvtcontrol):
#> Warning from mvtnorm::pmvt: Completion with error > abseps.
#> Warning in MCTpval(contMat, corMat, df, tStat, alternative, mvtcontrol):
#> Warning from mvtnorm::pmvt: Completion with error > abseps.
#> Warning in MCTpval(contMat, corMat, df, tStat, alternative, mvtcontrol):
#> Warning from mvtnorm::pmvt: Completion with error > abseps.
#> Warning in MCTpval(contMat, corMat, df, tStat, alternative, mvtcontrol):
#> Warning from mvtnorm::pmvt: Completion with error > abseps.
#> Warning in MCTpval(contMat, corMat, df, tStat, alternative, mvtcontrol):
#> Warning from mvtnorm::pmvt: Completion with error > abseps.
#> Warning in MCTpval(contMat, corMat, df, tStat, alternative, mvtcontrol):
#> Warning from mvtnorm::pmvt: Completion with error > abseps.
#> Warning in MCTpval(contMat, corMat, df, tStat, alternative, mvtcontrol):
#> Warning from mvtnorm::pmvt: Completion with error > abseps.
#> Warning in MCTpval(contMat, corMat, df, tStat, alternative, mvtcontrol):
#> Warning from mvtnorm::pmvt: Completion with error > abseps.
#> Warning in MCTpval(contMat, corMat, df, tStat, alternative, mvtcontrol):
#> Warning from mvtnorm::pmvt: Completion with error > abseps.
#> Warning in MCTpval(contMat, corMat, df, tStat, alternative, mvtcontrol):
#> Warning from mvtnorm::pmvt: Completion with error > abseps.
#> Warning in MCTpval(contMat, corMat, df, tStat, alternative, mvtcontrol):
#> Warning from mvtnorm::pmvt: Completion with error > abseps.
#> Warning in MCTpval(contMat, corMat, df, tStat, alternative, mvtcontrol):
#> Warning from mvtnorm::pmvt: Completion with error > abseps.
#> Warning in MCTpval(contMat, corMat, df, tStat, alternative, mvtcontrol):
#> Warning from mvtnorm::pmvt: Completion with error > abseps.
#> Warning in MCTpval(contMat, corMat, df, tStat, alternative, mvtcontrol):
#> Warning from mvtnorm::pmvt: Completion with error > abseps.
#> Warning in MCTpval(contMat, corMat, df, tStat, alternative, mvtcontrol):
#> Warning from mvtnorm::pmvt: Completion with error > abseps.
#> Warning in MCTpval(contMat, corMat, df, tStat, alternative, mvtcontrol):
#> Warning from mvtnorm::pmvt: Completion with error > abseps.
#> Warning in MCTpval(contMat, corMat, df, tStat, alternative, mvtcontrol):
#> Warning from mvtnorm::pmvt: Completion with error > abseps.
#> Warning in MCTpval(contMat, corMat, df, tStat, alternative, mvtcontrol):
#> Warning from mvtnorm::pmvt: Completion with error > abseps.
#> Warning in MCTpval(contMat, corMat, df, tStat, alternative, mvtcontrol):
#> Warning from mvtnorm::pmvt: Completion with error > abseps.
#> Warning in MCTpval(contMat, corMat, df, tStat, alternative, mvtcontrol):
#> Warning from mvtnorm::pmvt: Completion with error > abseps.
#> Warning in MCTpval(contMat, corMat, df, tStat, alternative, mvtcontrol):
#> Warning from mvtnorm::pmvt: Completion with error > abseps.
#> Warning in MCTpval(contMat, corMat, df, tStat, alternative, mvtcontrol):
#> Warning from mvtnorm::pmvt: Completion with error > abseps.
#> Warning in MCTpval(contMat, corMat, df, tStat, alternative, mvtcontrol):
#> Warning from mvtnorm::pmvt: Completion with error > abseps.
#> Warning in MCTpval(contMat, corMat, df, tStat, alternative, mvtcontrol):
#> Warning from mvtnorm::pmvt: Completion with error > abseps.
#> Warning in MCTpval(contMat, corMat, df, tStat, alternative, mvtcontrol):
#> Warning from mvtnorm::pmvt: Completion with error > abseps.
#> Warning in MCTpval(contMat, corMat, df, tStat, alternative, mvtcontrol):
#> Warning from mvtnorm::pmvt: Completion with error > abseps.
#> Warning in MCTpval(contMat, corMat, df, tStat, alternative, mvtcontrol):
#> Warning from mvtnorm::pmvt: Completion with error > abseps.
#> Warning in MCTpval(contMat, corMat, df, tStat, alternative, mvtcontrol):
#> Warning from mvtnorm::pmvt: Completion with error > abseps.
#> Warning in MCTpval(contMat, corMat, df, tStat, alternative, mvtcontrol):
#> Warning from mvtnorm::pmvt: Completion with error > abseps.
#> Warning in MCTpval(contMat, corMat, df, tStat, alternative, mvtcontrol):
#> Warning from mvtnorm::pmvt: Completion with error > abseps.
#> Warning in MCTpval(contMat, corMat, df, tStat, alternative, mvtcontrol):
#> Warning from mvtnorm::pmvt: Completion with error > abseps.
#> Warning in MCTpval(contMat, corMat, df, tStat, alternative, mvtcontrol):
#> Warning from mvtnorm::pmvt: Completion with error > abseps.
#> Warning in MCTpval(contMat, corMat, df, tStat, alternative, mvtcontrol):
#> Warning from mvtnorm::pmvt: Completion with error > abseps.
#> Warning in MCTpval(contMat, corMat, df, tStat, alternative, mvtcontrol):
#> Warning from mvtnorm::pmvt: Completion with error > abseps.
#> Warning in MCTpval(contMat, corMat, df, tStat, alternative, mvtcontrol):
#> Warning from mvtnorm::pmvt: Completion with error > abseps.
#> Warning in MCTpval(contMat, corMat, df, tStat, alternative, mvtcontrol):
#> Warning from mvtnorm::pmvt: Completion with error > abseps.
#> Warning in MCTpval(contMat, corMat, df, tStat, alternative, mvtcontrol):
#> Warning from mvtnorm::pmvt: Completion with error > abseps.
#> Warning in MCTpval(contMat, corMat, df, tStat, alternative, mvtcontrol):
#> Warning from mvtnorm::pmvt: Completion with error > abseps.
#> Warning in MCTpval(contMat, corMat, df, tStat, alternative, mvtcontrol):
#> Warning from mvtnorm::pmvt: Completion with error > abseps.
matplot(
  pw_df5$emax,
  pw_df5[, c("power_dunnett", "power_mcpmod")],
  type = "b", pch = 19, lty = 1,
  col  = c("steelblue", "tomato"),
  xlab = "True emax (maximum effect)",
  ylab = "Empirical detection rate (alpha = 0.05)",
  main = "Dunnett vs MCP-Mod: detection rate across emax",
  ylim = c(0, 1)
)
legend("bottomright",
       legend = c("Dunnett", "MCP-Mod"),
       col = c("steelblue", "tomato"), lty = 1, pch = 19)
abline(h = 0.80, lty = 2, col = "grey50")

Next steps

  • Example 4 - Bayesian Go/No-Go with historical placebo borrowing
  • Population - endpoint data setup for continuous, binary, and time-to-event outcomes