Clinical Trial Simulation: Operating Characteristics Under Dropout, Nonadherence, and Site Heterogeneity
Md Ashraful Islam Bhuiya
2026-07-07
1 Project purpose
This project demonstrates a clinical-trial simulation workflow for a time-to-event endpoint. The goal is to show how simulation can be used to evaluate trial operating characteristics before a study is conducted. If we designed a realistic clinical trial with randomization, staggered enrollment, censoring, dropout, site differences, and imperfect adherence, how often would the trial correctly detect a treatment effect?
The simulation focuses on questions that are common in clinical development and trial planning:
- What power is expected under a target treatment effect?
- How does dropout affect power?
- How does site-level heterogeneity affect performance?
- How does treatment nonadherence attenuate the observed treatment effect?
- What type I error is observed under the null scenario?
- How can simulation outputs be summarized for decision support?
This is an example using simulated data only. It is not based on any real patient data, or proprietary clinical-trial data.
2 Packages
required_packages <- c(
"survival", "dplyr", "purrr", "ggplot2", "broom", "knitr", "tidyr", "scales"
)
new_packages <- required_packages[!(required_packages %in% installed.packages()[, "Package"])]
if (length(new_packages) > 0) {
install.packages(new_packages, repos = "https://cloud.r-project.org")
}
invisible(lapply(required_packages, library, character.only = TRUE))3 Data-generating model
The simulated trial has:
- 1:1 randomization
- staggered enrollment over an accrual period
- administrative censoring at study end
- random dropout
- site-level heterogeneity using log-hazard site effects
- treatment nonadherence in the investigational arm
- a time-to-event endpoint analyzed by Cox proportional hazards regression
The primary analysis uses an intention-to-treat style Cox model based on randomized assignment. Nonadherence is introduced in the data-generating mechanism to show how operational factors can attenuate the observed randomized-arm effect.
simulate_one_trial <- function(
n_total = 400,
n_sites = 30,
accrual_months = 12,
study_duration_months = 24,
control_median_months = 11.5,
treatment_hr = 0.72,
annual_dropout_probability = 0.10,
nonadherence_probability = 0.00,
site_sd_log_hazard = 0.20,
biomarker_prevalence = 0.40,
biomarker_hr = 1.40,
age_hr_per_10yr = 1.12
) {
stopifnot(study_duration_months > accrual_months)
patient_id <- seq_len(n_total)
site_id <- sample(seq_len(n_sites), size = n_total, replace = TRUE)
site_effect <- rnorm(n_sites, mean = 0, sd = site_sd_log_hazard)
randomized_trt <- rbinom(n_total, size = 1, prob = 0.5)
adherent_to_active <- ifelse(
randomized_trt == 1,
rbinom(n_total, size = 1, prob = 1 - nonadherence_probability),
0
)
treatment_received <- ifelse(randomized_trt == 1 & adherent_to_active == 1, 1, 0)
age <- pmin(pmax(round(rnorm(n_total, mean = 62, sd = 9)), 35), 85)
biomarker_high <- rbinom(n_total, size = 1, prob = biomarker_prevalence)
enrollment_month <- runif(n_total, min = 0, max = accrual_months)
administrative_followup <- study_duration_months - enrollment_month
baseline_hazard <- log(2) / control_median_months
linear_predictor <-
log(treatment_hr) * treatment_received +
log(biomarker_hr) * biomarker_high +
log(age_hr_per_10yr) * ((age - 60) / 10) +
site_effect[site_id]
event_time <- rexp(n_total, rate = baseline_hazard * exp(linear_predictor))
monthly_dropout_rate <- -log(1 - annual_dropout_probability) / 12
dropout_time <- rexp(n_total, rate = monthly_dropout_rate)
observed_time <- pmin(event_time, dropout_time, administrative_followup)
event_status <- as.integer(event_time <= dropout_time & event_time <= administrative_followup)
data.frame(
patient_id = patient_id,
site_id = factor(site_id),
randomized_trt = randomized_trt,
treatment_received = treatment_received,
age = age,
biomarker_high = biomarker_high,
enrollment_month = enrollment_month,
time_months = pmax(observed_time, 0.001),
event = event_status
)
}
analyze_one_trial <- function(dat) {
fit <- tryCatch(
coxph(Surv(time_months, event) ~ randomized_trt, data = dat),
error = function(e) NULL
)
if (is.null(fit)) {
return(data.frame(
beta_hat = NA_real_,
se_hat = NA_real_,
hr_hat = NA_real_,
p_value = NA_real_,
ci_low = NA_real_,
ci_high = NA_real_,
n_events = sum(dat$event),
event_rate = mean(dat$event)
))
}
fit_summary <- broom::tidy(fit, conf.int = TRUE, exponentiate = FALSE)
row <- fit_summary[fit_summary$term == "randomized_trt", ]
data.frame(
beta_hat = row$estimate,
se_hat = row$std.error,
hr_hat = exp(row$estimate),
p_value = row$p.value,
ci_low = exp(row$conf.low),
ci_high = exp(row$conf.high),
n_events = sum(dat$event),
event_rate = mean(dat$event)
)
}
simulate_trial_scenario <- function(
scenario_name,
n_sim = 500,
n_total = 400,
treatment_hr = 0.72,
annual_dropout_probability = 0.10,
nonadherence_probability = 0.00,
site_sd_log_hazard = 0.20,
alpha = 0.05
) {
sim_results <- purrr::map_dfr(seq_len(n_sim), function(i) {
dat <- simulate_one_trial(
n_total = n_total,
treatment_hr = treatment_hr,
annual_dropout_probability = annual_dropout_probability,
nonadherence_probability = nonadherence_probability,
site_sd_log_hazard = site_sd_log_hazard
)
analyze_one_trial(dat) %>%
mutate(sim_id = i)
})
true_log_hr <- log(treatment_hr)
sim_results %>%
summarise(
scenario = scenario_name,
n_sim = n_sim,
n_total = n_total,
treatment_hr = treatment_hr,
annual_dropout_probability = annual_dropout_probability,
nonadherence_probability = nonadherence_probability,
site_sd_log_hazard = site_sd_log_hazard,
mean_events = mean(n_events, na.rm = TRUE),
mean_event_rate = mean(event_rate, na.rm = TRUE),
mean_beta_hat = mean(beta_hat, na.rm = TRUE),
mean_hr_hat = mean(hr_hat, na.rm = TRUE),
empirical_power_or_type1 =
mean(p_value < alpha & beta_hat < 0, na.rm = TRUE),
two_sided_rejection_rate =
mean(p_value < alpha, na.rm = TRUE),
monte_carlo_se =
sqrt(empirical_power_or_type1 * (1 - empirical_power_or_type1) / n_sim),
beta_error_vs_generating_log_hr =
mean(beta_hat - true_log_hr, na.rm = TRUE),
ci_coverage_generating_hr =
mean(ci_low <= treatment_hr & ci_high >= treatment_hr, na.rm = TRUE),
failed_fits = sum(is.na(beta_hat)),
.groups = "drop"
)
}4 Define simulation scenarios
n_sim <- 500
scenario_grid <- tibble::tribble(
~scenario_name, ~n_total, ~treatment_hr, ~annual_dropout_probability, ~nonadherence_probability, ~site_sd_log_hazard,
"Base design", 400, 0.72, 0.10, 0.00, 0.20,
"Smaller sample size", 250, 0.72, 0.10, 0.00, 0.20,
"Larger sample size", 650, 0.72, 0.10, 0.00, 0.20,
"High dropout", 400, 0.72, 0.30, 0.00, 0.20,
"High site heterogeneity", 400, 0.72, 0.10, 0.00, 0.50,
"20 percent nonadherence", 400, 0.72, 0.10, 0.20, 0.20,
"Null effect for type I error", 400, 1.00, 0.10, 0.00, 0.20
)
scenario_grid## # A tibble: 7 × 6
## scenario_name n_total treatment_hr annual_dropout_probability
## <chr> <dbl> <dbl> <dbl>
## 1 Base design 400 0.72 0.1
## 2 Smaller sample size 250 0.72 0.1
## 3 Larger sample size 650 0.72 0.1
## 4 High dropout 400 0.72 0.3
## 5 High site heterogeneity 400 0.72 0.1
## 6 20 percent nonadherence 400 0.72 0.1
## 7 Null effect for type I error 400 1 0.1
## # ℹ 2 more variables: nonadherence_probability <dbl>, site_sd_log_hazard <dbl>5 Run the clinical-trial simulation
simulation_summary <- pmap_dfr(
scenario_grid,
function(
scenario_name,
n_total,
treatment_hr,
annual_dropout_probability,
nonadherence_probability,
site_sd_log_hazard
) {
simulate_trial_scenario(
scenario_name = scenario_name,
n_sim = n_sim,
n_total = n_total,
treatment_hr = treatment_hr,
annual_dropout_probability = annual_dropout_probability,
nonadherence_probability = nonadherence_probability,
site_sd_log_hazard = site_sd_log_hazard
)
}
)
simulation_summary## scenario n_sim n_total treatment_hr
## 1 Base design 500 400 0.72
## 2 Smaller sample size 500 250 0.72
## 3 Larger sample size 500 650 0.72
## 4 High dropout 500 400 0.72
## 5 High site heterogeneity 500 400 0.72
## 6 20 percent nonadherence 500 400 0.72
## 7 Null effect for type I error 500 400 1.00
## annual_dropout_probability nonadherence_probability site_sd_log_hazard
## 1 0.1 0.0 0.2
## 2 0.1 0.0 0.2
## 3 0.1 0.0 0.2
## 4 0.3 0.0 0.2
## 5 0.1 0.0 0.5
## 6 0.1 0.2 0.2
## 7 0.1 0.0 0.2
## mean_events mean_event_rate mean_beta_hat mean_hr_hat
## 1 244.102 0.6102550 -0.3207943591 0.7312312
## 2 152.908 0.6116320 -0.3165870024 0.7373664
## 3 396.532 0.6100492 -0.3251641599 0.7263536
## 4 212.126 0.5303150 -0.3142571130 0.7378295
## 5 243.506 0.6087650 -0.2947592792 0.7509343
## 6 248.320 0.6208000 -0.2463417062 0.7879194
## 7 265.968 0.6649200 -0.0003237773 1.0073740
## empirical_power_or_type1 two_sided_rejection_rate monte_carlo_se
## 1 0.708 0.708 0.020334011
## 2 0.486 0.486 0.022351913
## 3 0.888 0.888 0.014103617
## 4 0.620 0.620 0.021707142
## 5 0.644 0.644 0.021413267
## 6 0.478 0.478 0.022339024
## 7 0.028 0.054 0.007377805
## beta_error_vs_generating_log_hr ci_coverage_generating_hr failed_fits
## 1 0.0077097079 0.954 0
## 2 0.0119170645 0.948 0
## 3 0.0033399070 0.940 0
## 4 0.0142469540 0.942 0
## 5 0.0337447877 0.936 0
## 6 0.0821623608 0.916 0
## 7 -0.0003237773 0.946 06 Summary table to support decision
display_table <- simulation_summary %>%
mutate(
empirical_power_or_type1 = round(empirical_power_or_type1, 3),
monte_carlo_se = round(monte_carlo_se, 3),
mean_hr_hat = round(mean_hr_hat, 3),
mean_events = round(mean_events, 1),
mean_event_rate = round(mean_event_rate, 3),
beta_error_vs_generating_log_hr = round(beta_error_vs_generating_log_hr, 3),
ci_coverage_generating_hr = round(ci_coverage_generating_hr, 3)
) %>%
select(
scenario,
n_total,
treatment_hr,
annual_dropout_probability,
nonadherence_probability,
site_sd_log_hazard,
mean_events,
mean_event_rate,
mean_hr_hat,
empirical_power_or_type1,
monte_carlo_se,
ci_coverage_generating_hr,
failed_fits
)
knitr::kable(
display_table,
caption = "Clinical-trial simulation operating characteristics across scenarios."
)| scenario | n_total | treatment_hr | annual_dropout_probability | nonadherence_probability | site_sd_log_hazard | mean_events | mean_event_rate | mean_hr_hat | empirical_power_or_type1 | monte_carlo_se | ci_coverage_generating_hr | failed_fits |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Base design | 400 | 0.72 | 0.1 | 0.0 | 0.2 | 244.1 | 0.610 | 0.731 | 0.708 | 0.020 | 0.954 | 0 |
| Smaller sample size | 250 | 0.72 | 0.1 | 0.0 | 0.2 | 152.9 | 0.612 | 0.737 | 0.486 | 0.022 | 0.948 | 0 |
| Larger sample size | 650 | 0.72 | 0.1 | 0.0 | 0.2 | 396.5 | 0.610 | 0.726 | 0.888 | 0.014 | 0.940 | 0 |
| High dropout | 400 | 0.72 | 0.3 | 0.0 | 0.2 | 212.1 | 0.530 | 0.738 | 0.620 | 0.022 | 0.942 | 0 |
| High site heterogeneity | 400 | 0.72 | 0.1 | 0.0 | 0.5 | 243.5 | 0.609 | 0.751 | 0.644 | 0.021 | 0.936 | 0 |
| 20 percent nonadherence | 400 | 0.72 | 0.1 | 0.2 | 0.2 | 248.3 | 0.621 | 0.788 | 0.478 | 0.022 | 0.916 | 0 |
| Null effect for type I error | 400 | 1.00 | 0.1 | 0.0 | 0.2 | 266.0 | 0.665 | 1.007 | 0.028 | 0.007 | 0.946 | 0 |
7 Power and type I error plot
In not null, empirical_power_or_type1 is empirical
power. For the null scenario, it is the one-sided false-positive rate in
the favorable direction. The two_sided_rejection_rate
column can be used to evaluate conventional two-sided type I error.
plot_df <- simulation_summary %>%
mutate(
scenario = factor(scenario, levels = scenario_grid$scenario_name),
metric_label = ifelse(treatment_hr == 1, "Type I error scenario", "Power scenario")
)
ggplot(plot_df, aes(x = scenario, y = empirical_power_or_type1)) +
geom_col() +
geom_errorbar(
aes(
ymin = pmax(0, empirical_power_or_type1 - 1.96 * monte_carlo_se),
ymax = pmin(1, empirical_power_or_type1 + 1.96 * monte_carlo_se)
),
width = 0.2
) +
coord_flip() +
scale_y_continuous(labels = scales::percent_format(accuracy = 1), limits = c(0, 1)) +
labs(
title = "Estimated Operating Characteristics by Trial Scenario",
subtitle = "Bars show estimated power or type I error with approximate Monte Carlo uncertainty",
x = "Scenario",
y = "Estimated probability"
) +
theme_minimal(base_size = 12)
8 Sensitivity analysis: sample size curve
This section estimates power across a simple range of sample sizes while holding other operational assumptions constant.
sample_size_grid <- tibble::tibble(
scenario_name = paste0("N = ", c(200, 300, 400, 500, 650, 800)),
n_total = c(200, 300, 400, 500, 650, 800),
treatment_hr = 0.72,
annual_dropout_probability = 0.10,
nonadherence_probability = 0.00,
site_sd_log_hazard = 0.20
)
sample_size_results <- pmap_dfr(
sample_size_grid,
function(
scenario_name,
n_total,
treatment_hr,
annual_dropout_probability,
nonadherence_probability,
site_sd_log_hazard
) {
simulate_trial_scenario(
scenario_name = scenario_name,
n_sim = n_sim,
n_total = n_total,
treatment_hr = treatment_hr,
annual_dropout_probability = annual_dropout_probability,
nonadherence_probability = nonadherence_probability,
site_sd_log_hazard = site_sd_log_hazard
)
}
)
sample_size_results %>%
select(scenario, n_total, mean_events, mean_hr_hat, empirical_power_or_type1, monte_carlo_se) %>%
knitr::kable(digits = 3, caption = "Power across sample-size scenarios.")| scenario | n_total | mean_events | mean_hr_hat | empirical_power_or_type1 | monte_carlo_se |
|---|---|---|---|---|---|
| N = 200 | 200 | 122.404 | 0.737 | 0.414 | 0.022 |
| N = 300 | 300 | 183.514 | 0.737 | 0.556 | 0.022 |
| N = 400 | 400 | 244.082 | 0.739 | 0.666 | 0.021 |
| N = 500 | 500 | 305.922 | 0.732 | 0.804 | 0.018 |
| N = 650 | 650 | 397.748 | 0.728 | 0.874 | 0.015 |
| N = 800 | 800 | 487.178 | 0.731 | 0.936 | 0.011 |
ggplot(sample_size_results, aes(x = n_total, y = empirical_power_or_type1)) +
geom_line() +
geom_point(size = 2) +
geom_errorbar(
aes(
ymin = pmax(0, empirical_power_or_type1 - 1.96 * monte_carlo_se),
ymax = pmin(1, empirical_power_or_type1 + 1.96 * monte_carlo_se)
),
width = 20
) +
scale_y_continuous(labels = scales::percent_format(accuracy = 1), limits = c(0, 1)) +
labs(
title = "Power Curve Across Sample Sizes",
subtitle = "Simulated time-to-event trial with Cox primary analysis",
x = "Total randomized sample size",
y = "Estimated power"
) +
theme_minimal(base_size = 12)
9 Interpretation template
This clinical-trial simulation demonstrates how operational assumptions can affect trial performance. In the base design, the simulation estimates the probability of detecting the target treatment effect under the assumed sample size, event rate, dropout, and site heterogeneity. The smaller and larger sample-size scenarios show how empirical power changes as the number of randomized participants increases. The high-dropout scenario demonstrates how loss to follow-up can reduce information by lowering the number of observed events. The high-site-heterogeneity scenario shows how unmodeled heterogeneity can increase variability. The nonadherence scenario shows attenuation of the estimated treatment effect when patients randomized to active treatment do not fully receive the active effect.
This type of simulation is useful because it connects design assumptions to operating characteristics such as power, type I error, bias, confidence-interval coverage, and expected event counts. In a real clinical-development setting, the data-generating assumptions would be calibrated using historical trial data, real-world evidence, operational benchmarks, disease knowledge, and stakeholder input.
10 References
Morris, T. P., White, I. R., & Crowther, M. J. (2019). Using simulation studies to evaluate statistical methods. Statistics in Medicine, 38(11), 2074-2102. https://doi.org/10.1002/sim.8086
Cox, D. R. (1972). Regression models and life-tables. Journal of the Royal Statistical Society: Series B, 34(2), 187-202. https://doi.org/10.1111/j.2517-6161.1972.tb00899.x
ICH. (2019). ICH E9(R1) addendum on estimands and sensitivity analysis in clinical trials to the guideline on statistical principles for clinical trials. https://database.ich.org/sites/default/files/E9-R1_Step4_Guideline_2019_1203.pdf
11 Reproducibility
sessionInfo()## R version 4.5.0 (2025-04-11 ucrt)
## Platform: x86_64-w64-mingw32/x64
## Running under: Windows 11 x64 (build 26200)
##
## Matrix products: default
## LAPACK version 3.12.1
##
## locale:
## [1] LC_COLLATE=English_United States.utf8
## [2] LC_CTYPE=English_United States.utf8
## [3] LC_MONETARY=English_United States.utf8
## [4] LC_NUMERIC=C
## [5] LC_TIME=English_United States.utf8
##
## time zone: America/Chicago
## tzcode source: internal
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] scales_1.4.0 tidyr_1.3.1 knitr_1.50 broom_1.0.8 ggplot2_3.5.2
## [6] purrr_1.2.2 dplyr_1.1.4 survival_3.8-3
##
## loaded via a namespace (and not attached):
## [1] Matrix_1.7-3 gtable_0.3.6 jsonlite_2.0.0 compiler_4.5.0
## [5] tidyselect_1.2.1 jquerylib_0.1.4 splines_4.5.0 yaml_2.3.10
## [9] fastmap_1.2.0 lattice_0.22-6 R6_2.6.1 labeling_0.4.3
## [13] generics_0.1.3 backports_1.5.0 tibble_3.2.1 bslib_0.9.0
## [17] pillar_1.10.2 RColorBrewer_1.1-3 rlang_1.2.0 utf8_1.2.4
## [21] cachem_1.1.0 xfun_0.52 sass_0.4.10 cli_3.6.4
## [25] withr_3.0.2 magrittr_2.0.3 digest_0.6.37 grid_4.5.0
## [29] rstudioapi_0.17.1 lifecycle_1.0.4 vctrs_0.7.3 evaluate_1.0.3
## [33] glue_1.8.0 farver_2.1.2 rmarkdown_2.29 tools_4.5.0
## [37] pkgconfig_2.0.3 htmltools_0.5.8.1