Bayesian Clinical Trial Design: From Prior to Approval

library(BOIN)

# Design parameters
target_dlt <- 0.30    # Target DLT rate (MTD definition)
ncohort <- 10         # Maximum cohorts
cohortsize <- 3       # Patients per cohort
n_doses <- 5          # Number of dose levels

# Step 1: Get decision boundaries
boundaries <- get.boundary(
  target = target_dlt,
  ncohort = ncohort,
  cohortsize = cohortsize
)

# Step 2: Simulate operating characteristics
true_dlt <- c(0.05, 0.10, 0.25, 0.35, 0.50)
oc <- get.oc(
  target = target_dlt,
  p.true = true_dlt,
  ncohort = ncohort,
  cohortsize = cohortsize,
  ntrial = 10000
)
# Dose 3 (25% DLT) selected ~55% of time
# <10% patients at doses above MTD

library(psborrow2)
library(cmdstanr)

analysis <- create_analysis_obj(
  data_matrix = combined_data,
  outcome = outcome_surv_exponential(
    time_var = "os_time",
    cens_var = "os_event",
    baseline_prior = prior_normal(0, 1000)
  ),
  borrowing = borrowing_hierarchical_commensurate(
    ext_flag_col = "is_external",
    tau_prior = prior_half_cauchy(0, 0.5)
  ),
  treatment = treatment_details(
    trt_flag_col = "treatment_arm",
    trt_prior = prior_normal(0, 2.5)
  )
)

result <- mcmc_sample(analysis,
  iter_warmup = 2000, iter_sampling = 4000, chains = 4
)
# Key outputs: HR posterior, P(HR < 1), ESS borrowed

# ============================================
# BAYESIAN CLINICAL TRIAL ANALYSIS TOOLKIT
# ============================================
library(BOIN)        # Dose-finding
library(psborrow2)   # Historical borrowing
library(gsDesign)    # Sample size
library(rstan)       # General Bayesian

# --- PRIOR ELICITATION ---
elicit_prior <- function(k, n, delta = 0.5) {
  alpha <- 1 + delta * k
  beta  <- 1 + delta * (n - k)
  list(
    alpha = alpha, beta = beta,
    ess = alpha + beta - 2,
    mean = alpha / (alpha + beta),
    ci = qbeta(c(0.025, 0.975), alpha, beta)
  )
}

# --- BAYESIAN SAMPLE SIZE ---
bayesian_sample_size <- function(
    prior_alpha, prior_beta,
    null_rate, alt_rate,
    threshold = 0.975,
    target_power = 0.80,
    n_sim = 10000
) {
  for (n in seq(20, 500, by = 5)) {
    null_successes <- rbinom(n_sim, n, null_rate)
    null_posts <- pbeta(null_rate,
      prior_alpha + null_successes,
      prior_beta + n - null_successes,
      lower.tail = FALSE)
    type1 <- mean(null_posts >= threshold)

    alt_successes <- rbinom(n_sim, n, alt_rate)
    alt_posts <- pbeta(null_rate,
      prior_alpha + alt_successes,
      prior_beta + n - alt_successes,
      lower.tail = FALSE)
    power <- mean(alt_posts >= threshold)

    if (type1 <= 0.05 && power >= target_power) {
      return(list(n = n, type1 = type1, power = power))
    }
  }
}

# --- PPoS ---
calculate_ppos <- function(
    current_successes, current_n, planned_n,
    prior_alpha, prior_beta,
    null_rate, threshold = 0.975,
    n_sim = 10000
) {
  post_alpha <- prior_alpha + current_successes
  post_beta  <- prior_beta + current_n - current_successes
  remaining  <- planned_n - current_n

  future_p <- rbeta(n_sim, post_alpha, post_beta)
  future_successes <- rbinom(n_sim, remaining, future_p)

  final_successes <- current_successes + future_successes
  final_alpha <- prior_alpha + final_successes
  final_beta  <- prior_beta + planned_n - final_successes

  final_posts <- pbeta(null_rate, final_alpha,
    final_beta, lower.tail = FALSE)
  mean(final_posts >= threshold)
}

# --- USAGE (REBYOTA) ---
prior <- elicit_prior(k = 25, n = 45, delta = 0.5)
cat("Prior: Beta(", prior$alpha, ",", prior$beta, ")\n")

ss <- bayesian_sample_size(
  prior$alpha, prior$beta,
  null_rate = 0.45, alt_rate = 0.65
)
cat("Required N:", ss$n, "\n")

ppos <- calculate_ppos(
  current_successes = 50, current_n = 90,
  planned_n = 180, prior$alpha, prior$beta,
  null_rate = 0.45
)
cat("PPoS:", round(ppos * 100, 1), "%\n")