Skip to contents
# Core simulation framework
library(rxsim)

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

# 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(777)

Early development trials often face small sample sizes where borrowing information from historical placebo data can increase efficiency and reduce required sample size. This example demonstrates a Bayesian Go/No-Go decision rule for a Phase IIa trial, where the placebo posterior is informed by historical data via a robust mixture prior, and the treatment posterior uses a non-informative prior. A “Go” decision is issued when the posterior probability of a meaningful treatment effect (Δ > δ) exceeds a predefined threshold γ.

Scenario

A two-arm, fixed design trial (placebo vs treatment) with a single continuous endpoint and a single final analysis.

  • Placebo arm uses historical borrowing through an informative prior derived from historical placebo data and robustified with a non-informative component.
  • Treatment arm uses a non-informative prior.
  • Go/No-Go rule at final analysis:

Go if P(Δ>δdata)γ, \text{Go if } P(\Delta > \delta \mid \text{data}) \ge \gamma,

where Δ=μTμP\Delta = \mu_T - \mu_P, with δ=0.1\delta = 0.1 and γ=0.8\gamma = 0.8.

delta = 0.1 defines the minimum treatment-placebo difference of clinical relevance — the effect must exceed this threshold for the drug to be worth advancing. gamma = 0.8 requires 80% posterior probability of exceeding that threshold before a Go decision is issued, balancing false-positive risk against the cost of a missed opportunity. The true simulated effect mu_treat = 0.30 lies above delta, so a well-powered trial should frequently result in a Go.

# Trial design
sample_size <- 80
allocation  <- c(1, 1)
arms        <- c("placebo", "treatment")

# Data-generating truth for this simulation example
mu_placebo_true <- 0.00
delta_true      <- 0.30  # true treatment - placebo mean difference
sigma_known     <- 1.0

# Enrollment/dropout profile (fixed design timeline)
enrollment <- list(
  end_time = c(4, 8),
  rate     = c(6, 10)
)
dropout <- list(
  end_time = c(4, 8),
  rate     = c(0, 0)
)

# Bayesian decision thresholds
delta <- 0.1
gamma <- 0.8

scenario <- tidyr::expand_grid(
  sample_size = sample_size,
  allocation  = list(allocation),
  delta_true  = delta_true,
  delta       = delta,
  gamma       = gamma
)

Historical placebo dataset and priors

We include a small historical placebo dataset as study-level means and sample sizes.

Studies H1, H2, and H3 contribute a combined 62 historical placebo observations and are summarised into an informative prior via postmix() — a Bayesian update of the non-informative starting prior with the pooled historical data. robustify() blends this informative component with a vague component (20% weight), producing a mixture prior that protects against prior-data conflict: if the current trial’s placebo behaves unexpectedly, the vague component limits how strongly the historical data pulls the posterior. The treatment arm uses prior_noninf (effectively a flat prior), so the treatment posterior is driven entirely by the current data.

hist_placebo <- data.frame(
  study = c("H1", "H2", "H3"),
  n     = c(24, 18, 20),
  mean  = c(-0.05, 0.02, 0.00)
)

hist_placebo
#>   study  n  mean
#> 1    H1 24 -0.05
#> 2    H2 18  0.02
#> 3    H3 20  0.00

# Pooled historical summary used to derive an informative placebo prior
n_hist_total <- sum(hist_placebo$n)
m_hist_pooled <- weighted.mean(hist_placebo$mean, hist_placebo$n)

# Non-informative prior (very small pseudo-sample size in mn parametrization)
# Used as treatment prior and as base for historical update.
prior_noninf <- mixnorm(
  c(1, 0, 1e-6),
  sigma = sigma_known,
  param = "mn"
)

# Informative placebo prior from historical placebo summary
prior_placebo_inf <- postmix(
  prior_noninf,
  n = n_hist_total,
  m = m_hist_pooled
)
#> Using default prior reference scale 1

# Robustified placebo prior to protect against prior-data conflict
prior_placebo <- robustify(
  prior_placebo_inf,
  weight = 0.2,
  mean   = 0,
  n      = 1,
  sigma  = sigma_known
)

# Treatment prior remains non-informative
prior_treatment <- prior_noninf

summary(prior_placebo)
#>        mean          sd        2.5%       50.0%       97.5% 
#> -0.01083871  0.46144620 -1.15034971 -0.01312827  1.15034975

Populations

population_generators <- list(
  placebo = function(n) {
    data.frame(
      id = 1:n,
      y = rnorm(n, mean = mu_placebo_true, sd = sigma_known),
      readout_time = 1
    )
  },
  treatment = function(n) {
    data.frame(
      id = 1:n,
      y = rnorm(n, mean = mu_placebo_true + delta_true, sd = sigma_known),
      readout_time = 1
    )
  }
)

Trial Parameters

enrollment_fn and dropout_fn supply the inter-arrival and time-to-dropout distributions for individual subjects. These are passed to replicate_trial(), which uses them when building each replicate’s enrollment schedule.

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

Final analysis trigger (Bayesian Go/No-Go)

Only one trigger is added at final analysis.

RBesT::pmixdiff(post_treat, post_placebo, delta, lower.tail = FALSE) evaluates P(μ_T − μ_P > δ | data) by numerical integration over the mixture posterior distributions for treatment and placebo. qmixdiff() computes quantiles of the same posterior difference distribution, yielding a 95% credible interval for Δ = μ_T − μ_P.

analysis_generators <- list(
  final = list(
    trigger = rlang::exprs(
      sum(!is.na(enroll_time)) >= !!sample_size
    ),
    analysis = function(df, timer) {
      dat <- df |>
        dplyr::filter(!is.na(enroll_time)) |>
        dplyr::mutate(arm = factor(arm, levels = c("placebo", "treatment")))

      y_p <- dat$y[dat$arm == "placebo"]
      y_t <- dat$y[dat$arm == "treatment"]

      n_p <- length(y_p)
      n_t <- length(y_t)

      m_p <- mean(y_p)
      m_t <- mean(y_t)

      # Posterior for each arm mean
      post_placebo <- RBesT::postmix(prior_placebo, n = n_p, m = m_p)
      post_treat   <- RBesT::postmix(prior_treatment, n = n_t, m = m_t)

      # Posterior probability for treatment effect exceeding delta
      prob_delta <- RBesT::pmixdiff(post_treat, post_placebo, delta, lower.tail = FALSE)

      # Posterior summaries for Delta = mu_T - mu_P
      delta_ci <- RBesT::qmixdiff(post_treat, post_placebo, c(0.025, 0.5, 0.975))

      data.frame(
        scenario,
        n_placebo = n_p,
        n_treatment = n_t,
        mean_placebo = m_p,
        mean_treatment = m_t,
        post_prob_delta = prob_delta,
        delta_q025 = delta_ci[1],
        delta_q500 = delta_ci[2],
        delta_q975 = delta_ci[3],
        decision_go = as.integer(prob_delta >= gamma),
        stringsAsFactors = FALSE
      )
    }
  )
)

Trial

trials <- replicate_trial(
  trial_name = "example_7_bayes_two_arm_fixed",
  sample_size = sample_size,
  arms = arms,
  allocation = allocation,
  enrollment = enrollment_fn,
  dropout = dropout_fn,
  analysis_generators = analysis_generators,
  population_generators = population_generators,
  n = 3
)

Simulate

run_trials(trials)
#> Using default prior reference scale 1
#> Using default prior reference scale 1
#> Using default prior reference scale 1
#> Using default prior reference scale 1
#> Using default prior reference scale 1
#> Using default prior reference scale 1
#> [[1]]
#> <Trial>
#>   Public:
#>     clone: function (deep = FALSE) 
#>     conditions: list
#>     initialize: function (name, seed = NULL, timer = NULL, population = list(), 
#>     locked_data: list
#>     name: example_7_bayes_two_arm_fixed_1
#>     population: list
#>     results: list
#>     run: function () 
#>     seed: NULL
#>     timer: Timer, R6
#> 
#> [[2]]
#> <Trial>
#>   Public:
#>     clone: function (deep = FALSE) 
#>     conditions: list
#>     initialize: function (name, seed = NULL, timer = NULL, population = list(), 
#>     locked_data: list
#>     name: example_7_bayes_two_arm_fixed_2
#>     population: list
#>     results: list
#>     run: function () 
#>     seed: NULL
#>     timer: Timer, R6
#> 
#> [[3]]
#> <Trial>
#>   Public:
#>     clone: function (deep = FALSE) 
#>     conditions: list
#>     initialize: function (name, seed = NULL, timer = NULL, population = list(), 
#>     locked_data: list
#>     name: example_7_bayes_two_arm_fixed_3
#>     population: list
#>     results: list
#>     run: function () 
#>     seed: NULL
#>     timer: Timer, R6

Results

collect_results() row-binds analysis outputs across all replicates and prepends replicate (integer index), timepoint (calendar time at which the analysis fired), and analysis (the analysis name) to each row. post_prob_delta is the key decision metric: values at or above gamma = 0.8 yield decision_go = 1 (Go) and values below yield decision_go = 0 (No-Go). The credible interval columns delta_q025, delta_q500, and delta_q975 characterise the posterior uncertainty about the treatment effect Δ; with n=80 and σ=1, the posterior is still relatively wide. In practice you would run thousands of replicates and report the Go probability under both the null (Δ = 0) and the alternative (e.g., Δ = 0.30) to assess the design’s operating characteristics. See Example 8 for a seamless design that builds on this Bayesian decision rule.

replicate_results <- collect_results(trials)
replicate_results
#>   replicate timepoint analysis sample_size allocation delta_true delta gamma
#> 1         1  78.01039    final          80       1, 1        0.3   0.1   0.8
#> 2         2  72.27821    final          80       1, 1        0.3   0.1   0.8
#> 3         3  89.63501    final          80       1, 1        0.3   0.1   0.8
#>   n_placebo n_treatment mean_placebo mean_treatment post_prob_delta  delta_q025
#> 1        40          40    0.1581546      0.3864720       0.8826182 -0.05183437
#> 2        40          40    0.1669418      0.1654918       0.5048670 -0.27784403
#> 3        40          40   -0.4034483      0.3106430       0.9825414  0.12871903
#>   delta_q500 delta_q975 decision_go
#> 1  0.3272105  0.6968002           1
#> 2  0.1023216  0.4720979           0
#> 3  0.5184149  0.9892527           1

decision_go = 1 indicates Go, and decision_go = 0 indicates No-Go under the criterion P(Δ>0.1data)0.8P(\Delta > 0.1 \mid \text{data}) \ge 0.8.