Simulating Clinical Trial Data in R

Introduction

Clinical-trial statisticians often need data before a real trial dataset exists. You may need to evaluate a sample-size calculation, test an analysis model, validate a TLF program, investigate a difficult missing-data scenario, or demonstrate how an analysis behaves under a particular treatment effect.

Simulation solves this problem by allowing us to construct artificial trial datasets from an explicitly defined data-generating mechanism.

In R, simulation can range from a few lines of code that generate a normal endpoint to a complete virtual trial containing treatment assignment, baseline prognostic variables, repeated measurements, treatment discontinuation, rescue therapy, missing data, censoring, and time-to-event outcomes.

Why Simulate Clinical Trial Data?

There are several distinct reasons to simulate clinical-trial data. The simulation design should follow the question rather than the other way around.

The Core Simulation Workflow

Simulation Begins With a Data-Generating Mechanism

Suppose a two-arm trial has a continuous primary endpoint \(Y\). A simple model might be:

where \(TRT_i\) is treatment assignment, \(X_i\) is a baseline prognostic covariate, and \(\epsilon_i\) represents unexplained variability.

If the treatment effect is \(\beta_1=-5\), then the simulation has deliberately encoded an average treatment difference of 5 units in favor of the active arm. That known truth gives us something against which to evaluate the analysis.

Reproducibility and Random Seeds

Random-number generation means that a simulation can produce different data on different runs. That is desirable when exploring variability, but it can make debugging difficult unless the random-number state is controlled.

set.seed(20260910)

n <- 200

z <- rnorm(n, mean = 0, sd = 1)

head(z)

The same seed and same sequence of random-number calls reproduce the same results in a given R environment.

For production simulation work, record the seed, R version, package versions, simulation parameters, and analysis code. Reproducibility is part of the simulation specification rather than an afterthought.

Generating Subject IDs

n <- 200

dat <- data.frame(
  USUBJID = sprintf("SUBJ-%03d", seq_len(n))
)

head(dat)

A stable subject identifier becomes especially important once the simulation contains multiple records per subject.

Simulating Treatment Assignment

For a simple 1:1 randomized trial:

set.seed(20260910)

n <- 200

dat <- data.frame(
  USUBJID = sprintf("SUBJ-%03d", seq_len(n)),
  TRT01P = sample(
    c("Placebo", "Active"),
    size = n,
    replace = TRUE,
    prob = c(0.5, 0.5)
  )
)

table(dat$TRT01P)

For a 2:1 allocation:

dat$TRT01P <- sample(
  c("Placebo", "Active"),
  size = nrow(dat),
  replace = TRUE,
  prob = c(1/3, 2/3)
)

Simulating Stratified Randomization

Suppose randomization is stratified by region. One simple approach is to randomize within each stratum:

dat$REGION <- sample(
  c("North", "South", "West"),
  size = nrow(dat),
  replace = TRUE,
  prob = c(0.4, 0.35, 0.25)
)

dat$TRT01P <- ave(
  seq_len(nrow(dat)),
  dat$REGION,
  FUN = function(i) {
    sample(
      c("Placebo", "Active"),
      length(i),
      replace = TRUE
    )
  }
)

For a serious operating-characteristic study, the randomization algorithm should reproduce the planned trial procedure more faithfully than this simple illustration.

Simulating Demographics

Clinical-trial datasets contain a mixture of continuous and categorical baseline variables.

set.seed(20260910)

dat$AGE <- round(
  rnorm(nrow(dat), mean = 62, sd = 10)
)

dat$SEX <- sample(
  c("Female", "Male"),
  size = nrow(dat),
  replace = TRUE,
  prob = c(0.48, 0.52)
)

dat$RACE <- sample(
  c("White", "Black or African American", "Asian", "Other"),
  size = nrow(dat),
  replace = TRUE,
  prob = c(0.62, 0.12, 0.16, 0.10)
)

The values should not simply be generated because they "look realistic." Parameters should ideally be informed by historical trials, published literature, registry data, feasibility work, or explicit assumptions.

Simulating a Continuous Baseline Variable

Suppose baseline disease burden has an approximately normal distribution:

dat$BASEVAL <- rnorm(
  nrow(dat),
  mean = 100,
  sd = 20
)

If the variable cannot be negative, a normal distribution may be inappropriate. Possible alternatives include log-normal, gamma, truncated-normal, or empirical distributions.

Choosing a Distribution

Simulating a Continuous Primary Endpoint

Suppose the endpoint is change from baseline at Week 12. Assume:

• Placebo mean change = −5
• Active mean change = −10
• Residual standard deviation = 15

Then a simple simulation is:

dat$CHG12 <- rnorm(
  nrow(dat),
  mean = ifelse(
    dat$TRT01P == "Active",
    -10,
    -5
  ),
  sd = 15
)

The expected treatment difference is:

Adding Baseline Prognostic Effects

Clinical endpoints are rarely determined by treatment alone. Suppose patients with higher baseline disease burden tend to have larger changes:

dat$CHG12 <- with(
  dat,
  5 +
  ifelse(TRT01P == "Active", -5, 0) +
  0.20 * (BASEVAL - 100) +
  rnorm(nrow(dat), 0, 15)
)

Centering the baseline variable around 100 makes the intercept easier to interpret.

Simulating a Binary Response Endpoint

Suppose the primary endpoint is objective response. Let:

• Placebo response probability = 0.25
• Active response probability = 0.45

A Bernoulli simulation is:

dat$RESP <- rbinom(
  nrow(dat),
  size = 1,
  prob = ifelse(
    dat$TRT01P == "Active",
    0.45,
    0.25
  )
)

dat$RESPFL <- ifelse(
  dat$RESP == 1,
  "Y",
  "N"
)

The expected response-rate difference is:

or 20 percentage points.

Simulating Binary Outcomes With Covariates

A more realistic model can use a logistic probability:

eta <- with(
  dat,
  -1.10 +
  0.65 * (TRT01P == "Active") +
  0.30 * ((AGE - 60) / 10)
)

p <- plogis(eta)

dat$RESP <- rbinom(
  nrow(dat),
  size = 1,
  prob = p
)

This approach is preferable when response probability should depend on baseline prognostic variables.

Simulating Count Data

Suppose the endpoint is the number of exacerbations during follow-up. A Poisson model can be written:

lambda <- ifelse(
  dat$TRT01P == "Active",
  0.80,
  1.20
)

dat$EVENTS <- rpois(
  nrow(dat),
  lambda = lambda
)

The Poisson distribution assumes that the conditional mean and variance are equal. Clinical event counts frequently exhibit overdispersion, in which case a negative-binomial data-generating process may be more appropriate.

Simulating Negative-Binomial Counts

library(MASS)

dat$EVENTS <- rnbinom(
  nrow(dat),
  size = 2,
  mu = ifelse(
    dat$TRT01P == "Active",
    0.80,
    1.20
  )
)

The negative-binomial model provides an additional dispersion parameter and can better represent heterogeneous event rates.

Simulating Longitudinal Clinical Trial Data

Many clinical-trial endpoints are measured repeatedly. A common structure is:

Here \(b_i\) is a subject-specific random intercept. Because the same patient contributes multiple observations, the measurements are correlated.

Creating a Longitudinal Dataset

set.seed(20260910)

n <- 200
visits <- c(0, 4, 8, 12, 16)

subject <- data.frame(
  USUBJID = sprintf("SUBJ-%03d", 1:n),
  TRT01P = sample(
    c("Placebo", "Active"),
    n,
    replace = TRUE
  )
)

longdat <- merge(
  subject,
  data.frame(WEEK = visits)
)

longdat <- longdat[
  order(longdat$USUBJID, longdat$WEEK),
]

Each patient now has one record per nominal assessment.

Adding a Random Subject Effect

subject_effect <- rnorm(
  n,
  mean = 0,
  sd = 8
)

names(subject_effect) <- subject$USUBJID

longdat$RANDOM_INTERCEPT <- subject_effect[
  longdat$USUBJID
]

longdat$Y <- with(
  longdat,
  50 +
  ifelse(TRT01P == "Active", -5, 0) +
  -0.8 * WEEK +
  RANDOM_INTERCEPT +
  rnorm(nrow(longdat), 0, 6)
)

Patients with positive random intercepts tend to remain above patients with negative random intercepts across visits, inducing within-subject correlation.

Adding a Treatment-by-Time Interaction

Treatment effects may emerge gradually rather than being constant across all visits:

longdat$Y <- with(
  longdat,
  50 +
  -0.6 * WEEK +
  ifelse(
    TRT01P == "Active",
    -0.05 * WEEK^2,
    0
  ) +
  RANDOM_INTERCEPT +
  rnorm(nrow(longdat), 0, 6)
)

This is a deliberately simple nonlinear example. In a real simulation, the functional form should be selected because it reflects the intended clinical trajectory, not merely because it creates an attractive plot.

An Actual Simulated Longitudinal Trial

The interactive figure below illustrates one simulated trial in which the active treatment produces progressively greater mean improvement while individual patients retain substantial variability.

Simulating Correlation Directly

For multivariate repeated measurements, one approach is to generate correlated normal observations using a covariance matrix.

Suppose five visits have standard deviation 10 and a common correlation of 0.60. The covariance matrix is:

visits <- 5
sigma <- 10
rho <- 0.60

Sigma <- outer(
  1:visits,
  1:visits,
  FUN = function(j, k) {
    sigma^2 * rho^abs(j - k)
  }
)

A multivariate normal draw can then be used to generate a patient's repeated residuals. Packages such as MASS can be used for multivariate-normal generation.

Simulating Baseline and Post-Baseline Values Together

A common mistake is to simulate baseline independently of every post-baseline measurement. If the analysis uses change from baseline, the relationship between baseline and follow-up should be considered.

One simple model is:

This explicitly allows baseline disease burden to explain variation in the follow-up outcome.

Simulating Change From Baseline

Alternatively, the simulation can directly generate change from baseline:

dat$CHG12 <- with(
  dat,
  -5 +
  ifelse(TRT01P == "Active", -5, 0) +
  0.15 * (BASEVAL - 100) +
  rnorm(nrow(dat), 0, 12)
)

dat$WEEK12 <- dat$BASEVAL + dat$CHG12

Whether it is preferable to simulate the absolute post-baseline value or the change depends on the scientific question and the planned analysis.

Simulating Time-to-Event Data

Survival endpoints require both an event time and an observation process.

For an exponential survival distribution:

The event rate can depend on treatment:

set.seed(20260910)

dat$RATE <- ifelse(
  dat$TRT01P == "Active",
  0.07,
  0.10
)

dat$EVENT_TIME <- rexp(
  nrow(dat),
  rate = dat$RATE
)

Lower hazard produces longer event times on average.

Simulating a Weibull Survival Distribution

The Weibull distribution is often more flexible because its hazard can increase or decrease over time.

shape <- 1.4

scale_placebo <- 10
scale_active  <- 13

dat$SCALE <- ifelse(
  dat$TRT01P == "Active",
  scale_active,
  scale_placebo
)

dat$EVENT_TIME <- rweibull(
  nrow(dat),
  shape = shape,
  scale = dat$SCALE
)

Adding Administrative Censoring

study_end <- 18

dat$CENSOR_TIME <- study_end

dat$AVAL <- pmin(
  dat$EVENT_TIME,
  dat$CENSOR_TIME
)

dat$EVENT <- as.integer(
  dat$EVENT_TIME <= dat$CENSOR_TIME
)

Now:

• AVAL is the observed follow-up time.
• EVENT indicates whether the event was observed.
• Patients whose event occurs after study end are censored at study end.

Adding Random Dropout to a Survival Simulation

Clinical trials can have censoring caused by reasons other than administrative study closure. For example, a patient may withdraw consent.

dat$DROPOUT_TIME <- rexp(
  nrow(dat),
  rate = 0.03
)

dat$OBS_END <- pmin(
  dat$CENSOR_TIME,
  dat$DROPOUT_TIME
)

dat$AVAL <- pmin(
  dat$EVENT_TIME,
  dat$OBS_END
)

dat$EVENT <- as.integer(
  dat$EVENT_TIME <= dat$OBS_END
)

This is still a simplified mechanism. In a more sophisticated simulation, dropout may depend on treatment, prior outcome history, adverse events, or unobserved disease status.

Simulating Missing Data

Missing data are among the most important features to simulate when evaluating clinical-trial analyses.

MCAR Example

set.seed(20260910)

miss <- runif(nrow(dat)) < 0.10

dat$CHG12[miss] <- NA

Approximately 10% of observations are missing completely at random.

MAR Example

Suppose patients with worse observed baseline disease are more likely to miss the Week 12 visit:

logit_pmiss <-
  -2.2 +
  0.45 * ((dat$BASEVAL - 100) / 10) +
  0.25 * (dat$TRT01P == "Placebo")

p_miss <- plogis(logit_pmiss)

miss <- runif(nrow(dat)) < p_miss

dat$CHG12[miss] <- NA

This creates a missingness mechanism dependent on observed characteristics.

MNAR Sensitivity Simulation

MNAR scenarios are often implemented by specifying how unobserved outcomes would differ from outcomes predicted under an MAR model.

For example, a delta-adjustment scenario can be expressed as:

where \(\delta\) represents a clinically meaningful departure from the MAR assumption.

delta <- -5

imputed_mnar <- imputed_mar + delta

The exact implementation should match the estimand, imputation model, and sensitivity-analysis strategy.

Simulating Dropout in Longitudinal Data

Instead of randomly deleting individual observations, it is often more realistic to simulate a subject-level dropout time.

set.seed(20260910)

dropout_week <- rexp(
  n,
  rate = 0.04
)

subject$DROPOUT_WEEK <- pmin(
  dropout_week,
  max(visits)
)

longdat <- merge(
  longdat,
  subject[c("USUBJID", "DROPOUT_WEEK")],
  by = "USUBJID"
)

longdat$AVAILABLE <-
  longdat$WEEK <= longdat$DROPOUT_WEEK

longdat$Y[!longdat$AVAILABLE] <- NA

This creates monotone dropout after a simulated time point.

Simulating Noncompliance

Randomization and treatment actually received are not necessarily identical. A simulation can therefore distinguish assigned treatment from exposure.

set.seed(20260910)

dat$COMPLIANT <- rbinom(
  nrow(dat),
  1,
  ifelse(dat$TRT01P == "Active", 0.90, 0.94)
)

dat$TRT01A <- ifelse(
  dat$TRT01P == "Active" & dat$COMPLIANT == 1,
  "Active",
  ifelse(
    dat$TRT01P == "Placebo" & dat$COMPLIANT == 1,
    "Placebo",
    "Other"
  )
)

Simulating Rescue Treatment

Suppose some patients receive rescue treatment after worsening disease.

longdat$RESCUE <- FALSE

eligible <- longdat$WEEK >= 8 &
            longdat$Y > 55

longdat$RESCUE[eligible] <-
  runif(sum(eligible)) < 0.30

A more realistic simulation would make rescue treatment persist after its initiation and would define exactly how the primary endpoint handles rescue.

Simulating Adverse Events

Adverse events can be simulated as binary, count, or recurrent-event outcomes. For a simple treatment-emergent adverse-event indicator:

ae_prob <- ifelse(
  dat$TRT01P == "Active",
  0.35,
  0.25
)

dat$AE_ANY <- rbinom(
  nrow(dat),
  1,
  ae_prob
)

dat$AE_ANYFL <- ifelse(
  dat$AE_ANY == 1,
  "Y",
  "N"
)

You can extend this by simulating preferred terms, severity, seriousness, relatedness, onset time, resolution time, and discontinuation.

Simulating Laboratory Data

Laboratory values often have repeated measurements and reference ranges. A simple example is:

lab <- expand.grid(
  USUBJID = dat$USUBJID,
  WEEK = c(0, 4, 8, 12),
  PARAM = "ALT"
)

lab <- merge(
  lab,
  dat[c("USUBJID", "TRT01P")],
  by = "USUBJID"
)

lab$AVAL <- rlnorm(
  nrow(lab),
  meanlog = log(28),
  sdlog = 0.35
)

lab$TOXFL <- ifelse(
  lab$AVAL > 80,
  "Y",
  "N"
)

For regulatory-grade simulations, reference ranges, baseline flags, treatment emergence, units, and clinically relevant toxicity grades should be represented explicitly.

Simulating an Oncology Tumor-Response Endpoint

A simplified oncology simulation can generate baseline tumor burden followed by percentage changes over time.

set.seed(20260910)

n <- 120
weeks <- c(0, 4, 8, 12, 16, 20)

subject <- data.frame(
  USUBJID = sprintf("SUBJ-%03d", 1:n),
  TRT01P = sample(
    c("Control", "Experimental"),
    n,
    replace = TRUE
  ),
  BASE_SUM <- rlnorm(
    n,
    meanlog = log(100),
    sdlog = 0.35
  )
)

For each post-baseline visit, a treatment-specific mean trajectory can be defined and patient-level random variation added.

trajectory_mean <- function(trt, week) {

  ifelse(
    trt == "Experimental",
    -2.0 * week + 0.025 * week^2,
    -0.5 * week + 0.015 * week^2
  )

}

tumor <- merge(
  subject,
  data.frame(WEEK = weeks)
)

tumor$PCHG <- ifelse(
  tumor$WEEK == 0,
  0,
  trajectory_mean(
    tumor$TRT01P,
    tumor$WEEK
  ) +
  rnorm(nrow(tumor), 0, 8)
)

This produces an illustrative tumor-response dataset that can be used to test waterfall or spider-plot programs.

Simulating Subject Heterogeneity

A realistic trial should usually contain meaningful patient-to-patient variation. If every patient receives exactly the same treatment response plus small independent noise, the simulation may be too homogeneous.

One way to introduce heterogeneity is to give each patient a latent treatment response:

subject$RESPONSE_SENSITIVITY <-
  rnorm(
    nrow(subject),
    mean = 0,
    sd = 6
  )

tumor <- merge(
  tumor,
  subject[c(
    "USUBJID",
    "RESPONSE_SENSITIVITY"
  )],
  by = "USUBJID"
)

tumor$PCHG <- tumor$PCHG +
  tumor$RESPONSE_SENSITIVITY

The resulting trajectories vary around the average treatment effect.

Simulating Treatment Effect Heterogeneity

Treatment effects themselves can vary by baseline characteristics.

The interaction term \(\beta_3\) means that the treatment effect depends on \(X_i\).

dat$CHG12 <- with(
  dat,
  -5 +
  ifelse(TRT01P == "Active", -5, 0) +
  ifelse(
    TRT01P == "Active",
    -0.10 * (AGE - 60),
    0
  ) +
  rnorm(nrow(dat), 0, 12)
)

This can be useful when investigating subgroup effects, treatment-effect heterogeneity, or covariate-adjusted analyses.

Simulating Center Effects

Multicenter trials can contain site-level heterogeneity.

set.seed(20260910)

n_sites <- 20

site <- data.frame(
  SITEID = sprintf("SITE-%02d", 1:n_sites),
  SITE_EFFECT = rnorm(
    n_sites,
    mean = 0,
    sd = 3
  )
)

dat$SITEID <- sample(
  site$SITEID,
  nrow(dat),
  replace = TRUE
)

dat <- merge(
  dat,
  site,
  by = "SITEID"
)

dat$CHG12 <- dat$CHG12 +
  dat$SITE_EFFECT

This is a simple random-site-effect mechanism. Actual site behavior may include enrollment imbalance, site-by-treatment interactions, screening failure, investigator effects, or center-specific missingness.

Simulating Screening and Enrollment

A clinical trial does not begin with exactly \(N\) randomized subjects. Patients are screened, some fail eligibility criteria, and eligible patients are randomized.

set.seed(20260910)

n_screened <- 500

screen <- data.frame(
  SCRNID = sprintf(
    "SCRN-%04d",
    1:n_screened
  )
)

screen$ELIGIBLE <- rbinom(
  n_screened,
  1,
  0.62
)

randomized <- subset(
  screen,
  ELIGIBLE == 1
)

nrow(randomized)

This lets you simulate a screening-to-randomization flow rather than assuming that every screened patient enrolls.

Simulating Early Treatment Discontinuation

Treatment discontinuation can depend on adverse events, progression, or other factors.

p_disc <- plogis(
  -2.0 +
  0.60 * dat$AE_ANY +
  0.25 * (dat$BASEVAL > 120)
)

dat$DISC <- rbinom(
  nrow(dat),
  1,
  p_disc
)

dat$DISC_REASON <- ifelse(
  dat$DISC == 1 & dat$AE_ANY == 1,
  "Adverse Event",
  ifelse(
    dat$DISC == 1,
    "Other",
    NA
  )
)

The simulation now connects one clinical process to another: adverse events influence discontinuation.

Building a More Realistic Trial Simulation

A useful simulation can be thought of as a causal chain:

This structure is much more informative than independently generating a collection of unrelated variables.

Simulating ADaM-Like Analysis Data

Simulation can be especially useful for clinical programmers when the goal is to test derivations. You can create an SDTM-like source layer, then derive analysis variables rather than generating the final analysis dataset directly.

For example:

adsl <- dat[c(
  "USUBJID",
  "TRT01P",
  "AGE",
  "SEX",
  "BASEVAL"
)]

adsl$SAFFL <- "Y"

adsl$ITTFL <- "Y"

adsl$PPROTFL <- ifelse(
  adsl$COMPLIANT == 1,
  "Y",
  "N"
)

The exact flags and derivations in a real project must follow the protocol, SAP, and dataset specifications.

Testing an Analysis Program With Known Truth

Suppose the simulation truth is a treatment difference of −5. You can fit the planned analysis model:

fit <- lm(
  CHG12 ~ TRT01P + BASEVAL,
  data = dat
)

summary(fit)

coef(fit)

If treatment is coded so that the active-arm coefficient represents the active versus placebo difference, the estimate should fluctuate around −5 over repeated simulations.

One Simulation Is Not Enough

A single simulated trial demonstrates one possible realization of the data. It does not tell you the long-run performance of the analysis.

For operating characteristics, repeat the complete trial many times.

For each repetition, record the analysis result:

• Treatment-effect estimate
• Standard error
• Confidence interval
• p-value
• Convergence status
• Whether the confidence interval contains the true effect

A Basic Monte Carlo Power Simulation

simulate_trial <- function(
  n_per_arm = 100,
  effect = -5,
  sd = 15
) {

  trt <- rep(
    c("Placebo", "Active"),
    each = n_per_arm
  )

  y <- rnorm(
    2 * n_per_arm,
    mean = ifelse(
      trt == "Active",
      effect,
      0
    ),
    sd = sd
  )

  fit <- lm(y ~ trt)

  p_value <- summary(fit)$coefficients[
    "trtActive",
    "Pr(>|t|)"
  ]

  estimate <- coef(fit)[
    "trtActive"
  ]

  c(
    estimate = estimate,
    p_value = p_value
  )
}

Repeat it:

set.seed(20260910)

B <- 2000

results <- replicate(
  B,
  simulate_trial(),
  simplify = "matrix"
)

results <- as.data.frame(t(results))

power <- mean(
  results$p_value < 0.05
)

power

The estimated power is the proportion of simulated trials in which the prespecified hypothesis test rejects the null hypothesis.

Estimating Monte Carlo Error

Simulation results themselves are estimates. If the estimated power is \(\hat p\) from \(B\) independent simulations, its Monte Carlo standard error is approximately:

For example, if estimated power is 0.80 from 2,000 simulations:

Thus, simulation-based power estimates should not be interpreted as infinitely precise numbers.

Simulation of Type I Error

To evaluate type I error, set the true treatment effect to zero:

set.seed(20260910)

type1 <- replicate(
  5000,
  simulate_trial(
    n_per_arm = 100,
    effect = 0,
    sd = 15
  )[["p_value"]]
)

mean(type1 < 0.05)

If the test is correctly calibrated under the simulation assumptions, the result should be near the nominal significance level, such as 0.05, allowing for Monte Carlo variability.

Coverage Probability

Suppose each simulation produces a 95% confidence interval for a true effect \(\theta\). Coverage is estimated by:

In R:

coverage <-
  mean(
    results$lower <= true_effect &
    results$upper >= true_effect
  )

A well-calibrated 95% interval procedure should generally have coverage close to 0.95 under the simulated conditions.

Bias

If \(\hat\theta_b\) is the treatment-effect estimate from simulation \(b\):

Estimate it by:

bias <-
  mean(results$estimate) -
  true_effect

A small bias does not guarantee that the method is good; variance, coverage, power, convergence, and other properties also matter.

An Actual Operating-Characteristic Example

Vectorization Versus Loops

For a single trial, vectorized R code is often simple and fast. For thousands of complete trials, a well-structured function is usually more important than micro-optimizing individual lines.

A clean pattern is:

simulate_trial <- function(seed, n = 200) {

  set.seed(seed)

  # 1. Generate subjects
  # 2. Randomize
  # 3. Generate outcomes
  # 4. Add missingness
  # 5. Analyze

  list(
    data = trial_data,
    analysis = analysis_result
  )
}

Then:

seeds <- 1:1000

sim_results <- lapply(
  seeds,
  simulate_trial
)

For very large simulations, packages such as future and future.apply can support parallel execution, provided the random-number streams are handled reproducibly.

Parallel Simulation

library(future.apply)

plan(multisession)

results <- future_lapply(
  1:2000,
  function(i) {
    simulate_trial(seed = i)
  },
  future.seed = TRUE
)

Simulation Parameters Should Be Centralized

A common programming mistake is scattering assumptions throughout the code. Instead, put the major assumptions into one parameter object.

sim_parms <- list(
  n_per_arm = 150,
  allocation = c(0.5, 0.5),
  treatment_effect = -5,
  residual_sd = 15,
  dropout_rate = 0.15,
  alpha = 0.05,
  n_sim = 2000
)

The simulation function can then accept this object:

simulate_trial <- function(parms) {

  # use parms$n_per_arm
  # use parms$treatment_effect
  # use parms$residual_sd
  # ...

}

This makes scenario management and sensitivity analyses much easier.

Scenario-Based Simulation

Rather than evaluating one assumed world, define multiple scenarios.

Scenario-based simulation is usually more informative than reporting a single power estimate based on one uncertain assumption.

Validating a Simulation

Before using simulated data to make a methodological decision, verify that the generated data match the intended assumptions.

Distributional Validation

Suppose the intended baseline mean is 100 and SD is 20. A single simulated sample may produce a mean of 98.7 and SD of 20.9. That does not mean the simulation is wrong.

The question is whether the generated distribution behaves correctly over repeated samples.

B <- 1000

sample_means <- replicate(
  B,
  mean(rnorm(200, 100, 20))
)

mean(sample_means)

sd(sample_means)

The standard deviation of the sample means should be close to the theoretical standard error:

Empirical Versus Theoretical Treatment Effect

Suppose the simulation encodes a treatment difference of −5. Over a large number of trials, the average estimated effect should approach −5.

mean(results$estimate)

sd(results$estimate)

quantile(
  results$estimate,
  c(0.025, 0.50, 0.975)
)

This is one of the simplest and most powerful checks that the simulation and analysis are communicating correctly.

Extreme-Value and Boundary Checks

Clinical variables often have natural constraints:

• Age should not be negative.
• Probabilities must be between 0 and 1.
• Counts must be non-negative integers.
• Time-to-event values must be non-negative.
• Laboratory values may require positive support.
• Percentage changes may have scientifically meaningful bounds depending on the endpoint.

Use explicit checks:

stopifnot(
  all(dat$AGE > 0),
  all(dat$BASEVAL > 0),
  all(dat$RESP %in% c(0, 1))
)

Simulation Diagnostics

A useful simulation workflow should produce diagnostics automatically.

Simulating Edge Cases for Programming QC

Simulation is particularly useful for testing clinical-programming logic. Randomly generated "typical" patients are not enough. Deliberately create edge cases.

• A patient with no post-baseline assessment
• A patient with all assessments missing after Week 4
• A subject with an extreme baseline value
• A patient who discontinues on the first treatment day
• A patient with an event at exactly the analysis cutoff
• A patient with duplicate dates
• A subject with an unscheduled assessment
• A patient with a treatment switch
• A subject with a zero count
• A subgroup containing very few patients

Testing a TLF With Simulated Data

A useful simulated dataset can drive the same programs that will eventually process the real study.

# Example structure

ADSL  -> demographic table
ADLB  -> laboratory summaries
ADEFF -> efficacy summaries
ADTTE -> Kaplan-Meier analysis
ADAE  -> adverse-event tables

The advantage is that you can test the entire programming chain before the actual database is available.

Simulation for Sample-Size Planning

Closed-form power calculations are attractive when their assumptions are reasonable. Simulation becomes especially valuable when the design includes:

• Repeated measures
• Non-normal endpoints
• Unequal allocation
• Dropout
• Interim analyses
• Nonstandard estimands
• Complex covariance structures
• Time-to-event endpoints
• Response-adaptive or stratified procedures
• Multiple testing

The principle is straightforward: simulate the exact trial design, apply the planned analysis, and estimate the operating characteristic empirically.

Power Under Dropout

Suppose a trial requires 200 evaluable subjects but 15% of randomized subjects are expected to have missing primary endpoints. Simulation can reveal whether the resulting analysis retains adequate power.

Importantly, dropout can do more than reduce sample size. If dropout is informative, it can introduce bias or alter the effective estimand.

Power Curves

Rather than evaluating one treatment effect, evaluate a grid:

effects <- c(-2, -3, -4, -5, -6, -7)

power_results <- data.frame(
  effect = effects,
  power = NA_real_
)

for (i in seq_along(effects)) {

  pvals <- replicate(
    1000,
    simulate_trial(
      n_per_arm = 100,
      effect = effects[i]
    )[["p_value"]]
  )

  power_results$power[i] <-
    mean(pvals < 0.05)
}

This produces a more informative relationship between assumed effect size and power than a single point estimate.

Simulating Interim Analyses

Group-sequential designs require special care because repeated looks at the data can affect type I error.

A simulation should reproduce:

• Information fractions
• Interim timing
• Stopping boundaries
• Final analysis rules
• Multiplicity adjustments

The simulated trial should stop when the prespecified rule is triggered rather than continuing as if the interim analysis never happened.

Simulation and Estimands

Modern clinical-trial simulation should begin with the estimand. Treatment assignment, treatment adherence, intercurrent events, and outcomes should be generated in a way that supports the treatment effect being estimated.

The simulation mechanism should therefore be aligned with the estimand rather than treating the estimand as something chosen after the data are generated.

Simulation of Counterfactual Outcomes

For causal-methodology work, it may be necessary to simulate potential outcomes:

where \(Y_i(1)\) is the outcome that would be observed under active treatment and \(Y_i(0)\) is the outcome that would be observed under control.

Only one is observed for each randomized subject, but both can be generated inside the simulation because they define the data-generating truth.

y0 <- rnorm(
  n,
  mean = baseline_mean,
  sd = 15
)

y1 <- y0 + treatment_effect

TRT <- rbinom(
  n,
  1,
  0.5
)

Y <- ifelse(
  TRT == 1,
  y1,
  y0
)

Simulating Correlated Endpoints

Some trials have multiple related outcomes, such as several biomarkers or multiple domains of a composite endpoint.

Their joint distribution can be represented with a covariance matrix:

The covariance matrix must be positive semidefinite. An invalid covariance matrix will cause multivariate-normal generation to fail or produce mathematically inconsistent results.

Simulation of Recurrent Events

Recurrent-event endpoints require more than one event per subject. A simple Poisson-process simulation can generate event counts over exposure time, while more sophisticated models can introduce frailty, event dependence, or terminal events.

rate <- ifelse(
  dat$TRT01P == "Active",
  0.08,
  0.11
)

dat$FOLLOWUP <- rexp(
  nrow(dat),
  rate = 0.04
)

dat$N_EVENTS <- rpois(
  nrow(dat),
  lambda = rate * dat$FOLLOWUP
)

This is only a starting point. Recurrent events often require explicit consideration of event history and exposure time.

Simulating Clustered Data

If patients are nested within sites, physicians, regions, or other clusters, correlation can arise from shared cluster effects.

where \(u_j\) is a cluster-specific random effect.

This can be important when evaluating mixed models, generalized estimating equations, or robust standard errors.

Common Simulation Mistakes

1. Generating variables independently when they should be related. Baseline disease, treatment response, dropout, and adverse events may need explicit dependencies.
2. Using unrealistic distributions. A normal distribution is not appropriate for every clinical variable.
3. Ignoring within-subject correlation. Independent repeated measurements can badly misrepresent longitudinal data.
4. Simulating only the analysis endpoint. If the real analysis depends on treatment discontinuation, rescue, or intermediate outcomes, those processes may need to be simulated too.
5. Adding missing data as an afterthought. The missingness mechanism can materially change bias and power.
6. Running too few simulations. Monte Carlo error can make operating-characteristic estimates unstable.
7. Failing to validate the simulation. Code can execute perfectly while generating the wrong scientific scenario.
8. Changing the analysis between simulation and real study. The simulation should generally use the prespecified analysis it is intended to evaluate.
9. Hard-coding assumptions throughout the script. Centralized parameters make scenario management and auditing much easier.
10. Ignoring model convergence. Failed model fits should be counted and investigated, not silently discarded.

Handling Model Convergence Failures

Complex simulated datasets can produce model-fitting failures. For example, logistic regression can experience separation, mixed models can have singular fits, and survival models can encounter sparse-event problems.

fit <- tryCatch(
  glm(
    RESP ~ TRT01P + AGE,
    data = dat,
    family = binomial()
  ),
  error = function(e) NULL
)

if (is.null(fit)) {
  status <- "FAILED"
} else {
  status <- "SUCCESS"
}

Do not simply remove failed simulations from the denominator. If 4% of simulated trials fail to converge, that itself may be an important operating characteristic.

Simulation Results Should Be Tidy

A useful Monte Carlo result dataset might contain:

SIM     ESTIMATE    SE      LOWER     UPPER     PVALUE    CONVERGED
1       -4.82       2.10    -8.94     -0.70     0.022     Y
2       -5.31       2.25    -9.72     -0.90     0.019     Y
3       -3.88       2.01    -7.82      0.06     0.054     Y

This structure makes it straightforward to calculate bias, coverage, power, quantiles, failure rates, and scenario comparisons.

Simulation Output Should Be Auditable

For important work, save:

• Simulation scenario identifier
• Random-number seed or reproducible seed scheme
• Number of simulations
• Data-generating parameters
• Analysis parameters
• Software versions
• Convergence status
• Warnings and errors
• Summary operating characteristics

A simulation result without its assumptions is difficult to interpret or reproduce.

A Reusable Simulation Function

simulate_trial <- function(
  n_per_arm = 100,
  treatment_effect = -5,
  residual_sd = 15,
  dropout = 0.15,
  seed = NULL
) {

  if (!is.null(seed)) {
    set.seed(seed)
  }

  n <- 2 * n_per_arm

  trt <- rep(
    c("Placebo", "Active"),
    each = n_per_arm
  )

  baseline <- rnorm(
    n,
    mean = 100,
    sd = 20
  )

  y <- 5 +
    ifelse(
      trt == "Active",
      treatment_effect,
      0
    ) +
    0.20 * (baseline - 100) +
    rnorm(n, 0, residual_sd)

  miss <- runif(n) < dropout

  y[miss] <- NA_real_

  data.frame(
    USUBJID = sprintf(
      "SUBJ-%03d",
      seq_len(n)
    ),
    TRT01P = trt,
    BASE = baseline,
    CHG = y
  )
}

This function is intentionally compact. A production simulation should be expanded according to the trial's actual data-generating mechanism.

Calling the Function

trial1 <- simulate_trial(
  n_per_arm = 150,
  treatment_effect = -5,
  residual_sd = 15,
  dropout = 0.15,
  seed = 20260910
)

head(trial1)

Scenario Grid With expand.grid()

scenarios <- expand.grid(
  treatment_effect = c(-3, -5, -7),
  dropout = c(0.10, 0.20),
  residual_sd = c(12, 15)
)

scenarios$SCENARIO <-
  seq_len(nrow(scenarios))

scenarios

Each row represents a distinct simulated world.

Why Simulation Is So Important for Clinical Trial Statisticians

Simulation Versus Bootstrapping

Simulation and bootstrapping are related but answer different questions.

Simulation Versus Resampling

A simulation creates new observations from an assumed probability model. Resampling methods such as the bootstrap repeatedly reuse or perturb observed data. Both are valuable, but they rely on different assumptions.

Simulation Should Be Calibrated to Real Data When Possible

If historical clinical-trial data are available, they can inform:

• Baseline distributions
• Correlation structures
• Dropout rates
• Visit completion
• Event rates
• Treatment discontinuation
• Adverse-event frequencies
• Subgroup prevalence

Historical data should inform assumptions without blindly reproducing one specific study's quirks.

Posterior Predictive and Bayesian Simulation

In Bayesian work, a fitted model can generate posterior predictive datasets. The same general principle applies: generate data from a probability model, apply the analysis, and examine whether the resulting behavior is plausible.

Simulation therefore plays a central role in Bayesian design evaluation, assurance, adaptive designs, and decision-theoretic analyses.

Simulation for Adaptive Designs

Adaptive trials can require simulation because the operating characteristics depend on decisions made during the trial.

For example:

A Simulation Specification

Before writing code, create a simulation specification.

Quality-Control Checklist

What a Good Simulation Report Should Contain

A simulation report should make it possible for another statistician to understand exactly what was simulated and why.

• Scientific objective
• Trial design
• Estimand
• Data-generating distributions
• Parameter values
• Dependencies and correlation structures
• Missing-data mechanisms
• Intercurrent-event assumptions
• Analysis method
• Number of simulation repetitions
• Random-number strategy
• Performance metrics
• Monte Carlo error or precision
• Software and package versions

Common Reporting Language

Key Differences Between a Simple and a Clinical-Trial Simulation

The Most Important Concept

The most important skill in clinical-trial simulation is not knowing every R random-number function. It is learning to translate a clinical-trial design into an explicit data-generating process.

If treatment is randomized, encode randomization.

If patients are repeatedly measured, encode correlation.

If treatment effects change over time, encode the time interaction.

If patients discontinue treatment, encode discontinuation.

If missingness depends on observed outcomes, encode that dependence.

If the analysis has a known target parameter, make that parameter explicit in the simulation.

Then repeatedly generate, analyze, and compare the results with the known truth.

References

Morris, T.P., White, I.R., and Crowther, M.J. (2019). Using simulation studies to evaluate statistical methods. Statistics in Medicine, 38, 2074–2102.

Burton, A., Altman, D.G., Royston, P., and Holder, R.L. (2006). The design of simulation studies in medical statistics. Statistics in Medicine, 25, 4279–4292.

Kahan, B.C., and Morris, T.P. (2018). Analysis of multicentre trials with continuous outcomes: when and how should we account for centre effects? Statistical Methods in Medical Research.

Little, R.J.A., D'Agostino, R., Cohen, M.L., et al. (2012). The prevention and treatment of missing data in clinical trials. New England Journal of Medicine, 367, 1355–1360.

ICH E9. Statistical Principles for Clinical Trials. International Council for Harmonisation.

ICH E9(R1). Addendum on Estimands and Sensitivity Analysis in Clinical Trials. International Council for Harmonisation.