Skip to contents
# Core simulation framework
library(rxsim)

# Bayesian utilities
library(RBesT)
#> This is RBesT version 1.8.2 (released 2025-04-25, git-sha b9dab00)
library(BayesianMCPMod)
library(DoseFinding)

# Helpers
library(dplyr)
#> 
#> Attaching package: 'dplyr'
#> The following objects are masked from 'package:stats':
#> 
#>     filter, lag
#> The following objects are masked from 'package:base':
#> 
#>     intersect, setdiff, setequal, union

set.seed(8808)

Seamless Phase IIa/IIb designs allow a single trial to serve both exploratory and confirmatory objectives, starting with a limited arm set and expanding enrollment based on interim evidence. This example simulates a seamless design where the trial begins as a two-arm study (placebo + one active dose), conducts a Bayesian Go/No-Go interim at time 4, then expands to a full multi-arm dose-finding phase in Phase IIb, concluding with a Bayesian MCP-Mod analysis across all dose groups.

Scenario

A seamless Ph2a/2b trial with one continuous endpoint:

  • Ph2a (early period): 2 arms (placebo + one active dose).
  • Ph2b (later period): expansion to multi-arm dose finding.
  • Interim: Bayesian Go/No-Go based on treatment effect for the Ph2a dose.
  • Final: Bayesian MCP-Mod across all doses.

dose_map defines five arms: d0 (placebo, 0 mg), d5 (5 mg), d10 (10 mg), d20 (20 mg), and d40 (40 mg). The Emax truth has emax = 0.9 and ed50 = 15, meaning the dose-response curve saturates around 30–40 mg. The interim Go/No-Go at time 4 uses delta_interim = 0.1 and gamma_interim = 0.8 (same convention as Example 7), while alpha_final = 0.05 controls the Type I error at the final Bayesian MCP-Mod step.

# Dose levels and arm labels
dose_map <- c(d0 = 0, d5 = 5, d10 = 10, d20 = 20, d40 = 40)
arm_names <- names(dose_map)
dose_levels <- as.numeric(dose_map)

# Truth for data generation (Emax shape)
e0    <- 0.0
emax  <- 0.9
ed50  <- 15
sigma <- 1.0
mean_fun <- function(d) e0 + emax * d / (ed50 + d)

# Interim Go/No-Go thresholds
delta_interim <- 0.1
gamma_interim <- 0.8

# Multiplicity level for final Bayesian MCP step
alpha_final <- 0.05

# Planned analysis times
interim_time <- 4
final_time   <- 12

Seamless enrollment schedule (2-arm -> multi-arm)

In this implementation, all arms are pre-defined, but enrollment switches from two-arm to multi-arm by time.

The enroll_plan data frame specifies integer subject counts (enroller) per arm per time unit. During time 1–4 (Ph2a), only d0 and d20 enroll (8 per time each). From time 5 onward (Ph2b), three new doses (d5, d10, d40) begin enrolling while the d20 arm tapers off; d0 continues at 8 per time throughout. dropper is set to 0 (no dropout). tr_timer$get_end_timepoint() confirms the last time point at which any enrollment occurs.

time <- 1:12

# Ph2a: only placebo + d20 active
# Ph2b: expand to d5, d10, d40 while placebo/d20 continue per plan
enroll_plan <- data.frame(
  time = rep(time, times = length(arm_names)),
  arm  = rep(arm_names, each = length(time)),
  enroller = c(
    # d0 placebo (total 80)
    8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 0, 0,
    # d5 (total 20)
    0, 0, 0, 0, 3, 3, 4, 3, 3, 4, 0, 0,
    # d10 (total 20)
    0, 0, 0, 0, 3, 3, 4, 3, 3, 4, 0, 0,
    # d20 (total 40)
    8, 8, 8, 8, 2, 2, 2, 2, 0, 0, 0, 0,
    # d40 (total 20)
    0, 0, 0, 0, 3, 3, 4, 3, 3, 4, 0, 0
  ),
  dropper = 0L
)

enroll_plan$enroller <- as.integer(enroll_plan$enroller)
enroll_plan$dropper  <- as.integer(enroll_plan$dropper)

tr_timer <- Timer$new(name = "timer_example_8")
add_timepoints(tr_timer, enroll_plan)

tr_timer$get_end_timepoint()
#> [1] 12

Populations

mk_pop() creates Population objects directly from pre-generated data (not from generator functions), computing each arm’s n from sum(enroller) in the enrollment plan. This approach pre-generates the full cohort up-front; Trial$new() then draws from each population’s pool according to the time-specific enrollment counts.

# Planned max N by arm from enrollment schedule
n_by_arm <- enroll_plan |>
  group_by(arm) |>
  summarise(n = sum(enroller), .groups = "drop")

mk_pop <- function(arm, n) {
  d <- dose_map[[arm]]
  Population$new(
    name = arm,
    data = data.frame(
      id = seq_len(n),
      dose = d,
      y = rnorm(n, mean = mean_fun(d), sd = sigma),
      readout_time = 1
    )
  )
}

pops <- lapply(seq_len(nrow(n_by_arm)), function(i) {
  mk_pop(n_by_arm$arm[i], n_by_arm$n[i])
})

Priors

For interim Go/No-Go, we use placebo borrowing from small historical placebo data and a non-informative prior for active.

For final Bayesian MCP-Mod, we pass arm-wise priors to BayesianMCPMod.

The same historical borrowing approach as Example 7 is used for the interim and final placebo priors. prior_list_final has five entries — Ctrl (the robustified historical placebo prior) and DG_1 through DG_4 (non-informative priors for d5, d10, d20, d40) — matching the five dose groups expected by BayesianMCPMod::assessDesign().

# Historical placebo summary data
hist_placebo <- data.frame(
  study = c("H1", "H2", "H3"),
  n     = c(20, 24, 18),
  mean  = c(-0.06, 0.00, 0.03)
)

n_hist <- sum(hist_placebo$n)
m_hist <- weighted.mean(hist_placebo$mean, hist_placebo$n)

# Base non-informative prior
prior_noninf <- RBesT::mixnorm(c(1, 0, 1e-6), sigma = sigma, param = "mn")

# Placebo borrowing prior for interim and final control arm
prior_placebo_inf <- RBesT::postmix(prior_noninf, n = n_hist, m = m_hist)
#> Using default prior reference scale 1
prior_placebo <- RBesT::robustify(
  prior_placebo_inf,
  weight = 0.2,
  mean = 0,
  n = 1,
  sigma = sigma
)

# Prior list for Bayesian MCP-Mod (Ctrl + 4 active groups)
prior_list_final <- list(
  Ctrl = prior_placebo,
  DG_1 = prior_noninf,
  DG_2 = prior_noninf,
  DG_3 = prior_noninf,
  DG_4 = prior_noninf
)

Analyses

We add two trigger conditions:

  1. time == interim_time: Bayesian Go/No-Go on d20 vs placebo.
  2. time == final_time: Bayesian MCP-Mod over all dose groups.

trigger_by_calendar(interim_time, tr_timer, ...) registers the interim analysis to fire when calendar time reaches interim_time = 4. The <<- operator in go_interim <<- as.integer(...) writes the interim Go/No-Go decision to the enclosing environment, making it available to downstream analyses or conditional logic in the final step. The final analysis calls BayesianMCPMod::assessDesign() with arm-level sample sizes (n_patients), the candidate dose-response models (mods), the five-arm prior list (prior_list), and arm-level posterior mean estimates (estimates_sim); alpha_crit_val = alpha_final sets the significance threshold for the Bayesian MCP contrast step.

# shared state between interim and final trigger
go_interim <- NA_integer_

# Candidate dose-response models for Bayesian MCP-Mod
mods <- DoseFinding::Mods(
  linear      = NULL,
  emax        = 20,
  exponential = 20,
  sigEmax     = c(20, 3),
  doses       = dose_levels
)

# Interim: Bayesian Go/No-Go (d20 vs placebo)
trigger_by_calendar(interim_time, tr_timer, analysis = function(df, current_time) {
  dat <- df |>
    dplyr::filter(!is.na(enroll_time), arm %in% c("d0", "d20"))

  y_p <- dat$y[dat$arm == "d0"]
  y_t <- dat$y[dat$arm == "d20"]

  n_p <- length(y_p)
  n_t <- length(y_t)
  m_p <- mean(y_p)
  m_t <- mean(y_t)

  post_p <- RBesT::postmix(prior_placebo, n = n_p, m = m_p)
  post_t <- RBesT::postmix(prior_noninf,  n = n_t, m = m_t)

  prob_delta <- RBesT::pmixdiff(post_t, post_p, delta_interim, lower.tail = FALSE)
  delta_q <- RBesT::qmixdiff(post_t, post_p, c(0.025, 0.5, 0.975))

  go_interim <<- as.integer(prob_delta >= gamma_interim)

  list(
    n_placebo = n_p,
    n_d20 = n_t,
    post_prob_delta = prob_delta,
    delta_q025 = delta_q[1],
    delta_q500 = delta_q[2],
    delta_q975 = delta_q[3],
    decision_go = go_interim
  )
})

# Final: Bayesian MCP-Mod across all dose groups
tr_timer$add_condition(
  time == final_time,
  analysis = function(df, time) {
    dat <- df |>
      dplyr::filter(!is.na(enroll_time))

    # Arm-wise summary at final
    summ <- dat |>
      dplyr::group_by(arm) |>
      dplyr::summarise(
        dose = unique(dose)[1],
        n = dplyr::n(),
        mu_hat = mean(y),
        .groups = "drop"
      ) |>
      dplyr::arrange(dose)

    # Order expected by BayesianMCPMod: Ctrl, DG_1 ...
    # with dose order 0, 5, 10, 20, 40.
    n_vec <- summ$n |> as.double()
    mu_vec <- summ$mu_hat

    # covariance of mean estimates (independent groups here)
    S_hat <- diag((sigma / sqrt(n_vec))^2)

    # Prepare estimates from observed data for assessDesign
    estimates_sim <- list(
      mu_hats = list(mu_vec),
      S_hats = list(S_hat)
    )

    # Run Bayesian MCP-Mod assessment with observed data
    bmcp <- BayesianMCPMod::assessDesign(
      n_patients = n_vec,
      mods = mods,
      prior_list = prior_list_final,
      estimates_sim = estimates_sim,
      alpha_crit_val = alpha_final,
      n_sim = 1
    )

    # Extract key results into row-bindable format
    data.frame(
      n_total_final = nrow(dat),
      max_eff = attr(bmcp, "maxEff")
    )
  },
  name = "final_bayesian_mcpmod"
)

Trial

trial <- Trial$new(
  name = "example_8_seamless_ph2a2b",
  seed = 8808,
  timer = tr_timer,
  population = pops
)

Simulate

trial$run()
#> Using default prior reference scale 1
#> Using default prior reference scale 1
#> Consider to provide 'contr' for your custom simulated data or analysis results.

Results

Interim and final analyses are stored under the corresponding time keys.

trial$results is a named list indexed by time-point keys ("time_4", "time_12"); each entry holds the named results from analyses registered for that time. The interim results (under "time_4") contain the Bayesian Go/No-Go metrics for d20 vs placebo. The final results (under "time_12") include a final_bayesian_mcpmod entry whose sub-list components can be inspected with names(trial$results[[final_key]]$final_bayesian_mcpmod). Note that this single-run demonstration shows the structural output; in practice you would replicate using clone_trial() + run_trials() and then collect_results() to aggregate operating characteristics across many simulated trials.

trial$results
#> $time_4
#> $time_4$cal_time_4
#> $time_4$cal_time_4$n_placebo
#> [1] 32
#> 
#> $time_4$cal_time_4$n_d20
#> [1] 32
#> 
#> $time_4$cal_time_4$post_prob_delta
#> [1] 0.9448993
#> 
#> $time_4$cal_time_4$delta_q025
#> [1] 0.01981831
#> 
#> $time_4$cal_time_4$delta_q500
#> [1] 0.4387862
#> 
#> $time_4$cal_time_4$delta_q975
#> [1] 0.8449165
#> 
#> $time_4$cal_time_4$decision_go
#> [1] 1
#> 
#> 
#> 
#> $time_12
#> $time_12$final_bayesian_mcpmod
#>   n_total_final max_eff
#> 1           180       1
final_key <- paste0("time_", final_time)

# top-level components returned by final analysis
names(trial$results[[final_key]]$final_bayesian_mcpmod)
#> [1] "n_total_final" "max_eff"

# Bayesian MCP-Mod object class
class(trial$results[[final_key]]$final_bayesian_mcpmod$bayesian_mcpmod)
#> [1] "NULL"

This example demonstrates the seamless setup and analysis flow in one simulation run: 2-arm Ph2a enrollment, interim Bayesian Go/No-Go check, then final Bayesian MCP-Mod after multi-arm expansion.