The code below generates 360 simulated participants. Randomization is
stratified by study site × baseline severity, with random block sizes of
4 or 6 within each stratum; in a real trial, block sizes should not be
disclosed to recruiters.
set.seed(20260811)
make_block_sequence <- function(n, block_sizes = c(4L, 6L)) {
allocation <- character(0)
while (length(allocation) < n) {
block_size <- sample(block_sizes, 1L)
block <- rep(c("Control", "Treatment"), each = block_size / 2L)
allocation <- c(allocation, sample(block))
}
allocation[seq_len(n)]
}
n <- 360L
trial <- data.frame(
id = sprintf("P%03d", seq_len(n)),
site = factor(sample(paste0("Site ", 1:4), n, replace = TRUE,
prob = c(0.28, 0.26, 0.24, 0.22))),
sex = factor(sample(c("Female", "Male"), n, replace = TRUE,
prob = c(0.58, 0.42))),
age = round(pmin(pmax(rnorm(n, 45, 13), 18), 75), 1)
)
site_baseline <- c("Site 1" = 0, "Site 2" = 0.8,
"Site 3" = -0.6, "Site 4" = 0.4)
trial$baseline <- round(
pmin(pmax(rnorm(n, 24, 5.5) + site_baseline[trial$site], 10), 40), 1
)
trial$severity <- factor(ifelse(trial$baseline >= 24, "High", "Low"),
levels = c("Low", "High"))
trial$stratum <- interaction(trial$site, trial$severity, drop = TRUE)
trial$arm <- NA_character_
for (stratum_i in levels(trial$stratum)) {
index <- which(trial$stratum == stratum_i)
trial$arm[index] <- make_block_sequence(length(index))
}
trial$arm <- factor(trial$arm, levels = c("Control", "Treatment"))
treatment <- as.integer(trial$arm == "Treatment")
# Week 12 continuous outcome: lower scores are better.
site_week12 <- c("Site 1" = 0, "Site 2" = 0.5,
"Site 3" = -0.5, "Site 4" = 0.2)
trial$week12_full <- 5 + 0.58 * trial$baseline - 2.5 * treatment +
site_week12[trial$site] + 0.25 * (trial$sex == "Male") + rnorm(n, 0, 5.2)
trial$week12_full <- round(pmin(pmax(trial$week12_full, 0), 45), 1)
# Poorer outcomes and the Treatment arm have higher missingness probabilities;
# complete values are used only as simulated truth for teaching purposes.
p_missing <- plogis(-2.25 + 0.35 * treatment +
0.085 * (trial$week12_full - 18) +
0.20 * (trial$severity == "High"))
trial$missing12 <- rbinom(n, 1, p_missing)
trial$week12 <- ifelse(trial$missing12 == 1, NA_real_, trial$week12_full)
trial$change <- trial$week12 - trial$baseline
# Response: at least a 30% reduction from baseline; treating missing observations
# as nonresponse is only an example of a composite strategy.
trial$responder_full <- as.integer(
(trial$baseline - trial$week12_full) / trial$baseline >= 0.30
)
trial$responder_nri <- ifelse(is.na(trial$week12), 0L, trial$responder_full)
# Safety outcome.
p_ae <- plogis(qlogis(0.13) + log(1.75) * treatment +
0.012 * (trial$age - 45))
trial$adverse_event <- rbinom(n, 1, p_ae)
# Recurrence/progression within one year, with independent censoring for loss to follow-up.
event_time <- rexp(n, rate = 0.0022 *
exp(log(0.68) * treatment + 0.020 * (trial$baseline - 24)))
dropout_time <- rexp(n, rate = 0.00075)
trial$time_days <- pmin(event_time, dropout_time, 365)
trial$status <- as.integer(event_time <= pmin(dropout_time, 365))
# Repeated scores and an overdispersed 12-week episode count; simulate these after
# the core outcomes to preserve the fixed primary results.
trial$week4 <- pmax(0, 3.5 + 0.78 * trial$baseline - 0.9 * treatment +
rnorm(n, 0, 4.7))
trial$week8 <- pmax(0, 4.0 + 0.66 * trial$baseline - 1.8 * treatment +
rnorm(n, 0, 5.0))
trial$week4[rbinom(n, 1, plogis(-3 + 0.2 * treatment)) == 1] <- NA
trial$week8[rbinom(n, 1, plogis(-2.7 + 0.25 * treatment)) == 1] <- NA
trial$exposure_weeks <- round(runif(n, 8, 12), 1)
individual_frailty <- rgamma(n, shape = 2, rate = 2)
episode_mean <- exp(log(0.16) + log(0.75) * treatment +
0.018 * (trial$baseline - 24)) *
trial$exposure_weeks * individual_frailty
trial$episode_count <- rpois(n, episode_mean)
stopifnot(nrow(trial) == n, !anyNA(trial$arm),
all(trial$status %in% 0:1), all(trial$time_days > 0))
allocation_table <- addmargins(table(trial$stratum, trial$arm))
allocation_table
