1 Project purpose
2 Simulate a randomized time-to-event trial
3 Summary
4 Kaplan-Meier estimation
5 Median survival by arm
6 Log-rank test
7 Cox proportional hazards model
8 Proportional-hazards diagnostics
9 Adjusted survival curves
10 Interpretation
11 References
12 Reproducibility

1 Project purpose

This project demonstrates a complete survival-analysis workflow using a simulated randomized clinical-trial dataset. The goal is to show working knowledge of time-to-event endpoints, censoring, Kaplan-Meier estimation, log-rank testing, Cox proportional hazards modeling, proportional-hazards diagnostics, adjusted effect estimation, and clear reporting.

The dataset is simulated for portfolio purpose. It contains no patient-level data and should not be interpreted as evidence about any real drug, disease, company, or clinical trial.

2 Simulate a randomized time-to-event trial

The example mimics a two-arm randomized trial with an overall-survival-like endpoint. The simulation includes:

1:1 randomization to control vs investigational treatment
staggered enrollment
administrative censoring
random dropout
baseline prognostic factors
site-level heterogeneity

required_packages <- c(
  "survival", "dplyr", "ggplot2", "broom", "knitr", "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))

#####################################################
#####################################################

simulate_survival_trial <- function(
  n = 420,
  n_sites = 24,
  accrual_months = 12,
  study_duration_months = 24,
  baseline_median_months = 11.5,
  treatment_hr = 0.70,
  biomarker_hr = 1.45,
  age_hr_per_10yr = 1.15,
  site_sd_log_hazard = 0.25,
  annual_dropout_probability = 0.10
) {
  stopifnot(study_duration_months > accrual_months)

  site_id <- sample(seq_len(n_sites), size = n, replace = TRUE)
  site_effect <- rnorm(n_sites, mean = 0, sd = site_sd_log_hazard)

  arm <- rbinom(n, size = 1, prob = 0.5)
  arm_label <- factor(ifelse(arm == 1, "Investigational", "Control"),
                      levels = c("Control", "Investigational"))

  age <- round(rnorm(n, mean = 62, sd = 9))
  age <- pmin(pmax(age, 35), 85)

  biomarker_high <- rbinom(n, size = 1, prob = 0.40)
  biomarker_label <- factor(ifelse(biomarker_high == 1, "High", "Low"),
                            levels = c("Low", "High"))

  region <- factor(sample(c("North", "South", "East", "West"), size = n, replace = TRUE))

  enrollment_month <- runif(n, min = 0, max = accrual_months)
  administrative_followup <- study_duration_months - enrollment_month

  baseline_hazard <- log(2) / baseline_median_months

  linear_predictor <-
    log(treatment_hr) * arm +
    log(biomarker_hr) * biomarker_high +
    log(age_hr_per_10yr) * ((age - 60) / 10) +
    site_effect[site_id]

  event_time <- rexp(n, rate = baseline_hazard * exp(linear_predictor))

  monthly_dropout_rate <- -log(1 - annual_dropout_probability) / 12
  dropout_time <- rexp(n, 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 = seq_len(n),
    site_id = factor(site_id),
    region = region,
    arm = arm_label,
    arm_numeric = arm,
    age = age,
    biomarker = biomarker_label,
    biomarker_high = biomarker_high,
    enrollment_month = enrollment_month,
    time_months = observed_time,
    event = event_status
  )
}

trial_dat <- simulate_survival_trial()

dplyr::glimpse(trial_dat)

## Rows: 420
## Columns: 11
## $ patient_id       <int> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16…
## $ site_id          <fct> 1, 6, 13, 15, 12, 4, 16, 5, 12, 2, 24, 19, 15, 14, 17…
## $ region           <fct> West, South, North, East, East, West, West, West, Nor…
## $ arm              <fct> Investigational, Control, Investigational, Control, I…
## $ arm_numeric      <int> 1, 0, 1, 0, 1, 1, 1, 1, 1, 1, 0, 0, 1, 0, 1, 1, 0, 1,…
## $ age              <dbl> 64, 63, 56, 69, 60, 62, 54, 55, 65, 58, 69, 70, 49, 5…
## $ biomarker        <fct> Low, Low, High, Low, High, Low, High, Low, Low, High,…
## $ biomarker_high   <int> 0, 0, 1, 0, 1, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0,…
## $ enrollment_month <dbl> 2.8782932, 0.8360315, 9.6757741, 4.2033894, 3.1663988…
## $ time_months      <dbl> 7.0039140, 0.5214482, 14.3242259, 6.0279233, 12.37876…
## $ event            <int> 1, 1, 0, 1, 1, 1, 0, 1, 0, 1, 0, 1, 0, 1, 1, 0, 1, 1,…

3 Summary

endpoint_summary <- trial_dat %>%
  group_by(arm) %>%
  summarise(
    n = n(),
    events = sum(event),
    censored = n - events,
    event_rate = events / n,
    median_followup_months = median(time_months),
    .groups = "drop"
  )

knitr::kable(
  endpoint_summary,
  digits = 3,
  caption = "Endpoint summary by randomized treatment arm."
)

Endpoint summary by randomized treatment arm.
arm	n	events	censored	event_rate	median_followup_months
Control	193	119	74	0.617	10.236
Investigational	227	112	115	0.493	12.244

4 Kaplan-Meier estimation

km_fit <- survfit(Surv(time_months, event) ~ arm, data = trial_dat)

km_fit

## Call: survfit(formula = Surv(time_months, event) ~ arm, data = trial_dat)
## 
##                       n events median 0.95LCL 0.95UCL
## arm=Control         193    119   10.9    9.24    13.7
## arm=Investigational 227    112   16.1   12.38    20.3

km_df <- broom::tidy(km_fit) %>%
  mutate(
    arm = sub("arm=", "", strata),
    lower = conf.low,
    upper = conf.high
  )

ggplot(km_df, aes(x = time, y = estimate, group = arm, linetype = arm)) +
  geom_step(linewidth = 0.9) +
  geom_step(aes(y = lower), linewidth = 0.3, alpha = 0.5) +
  geom_step(aes(y = upper), linewidth = 0.3, alpha = 0.5) +
  scale_y_continuous(labels = scales::percent_format(accuracy = 1), limits = c(0, 1)) +
  labs(
    title = "Kaplan-Meier Survival Curves",
    subtitle = "Simulated randomized clinical trial with administrative censoring and dropout",
    x = "Months since randomization",
    y = "Estimated survival probability",
    linetype = "Arm"
  ) +
  theme_minimal(base_size = 12)

5 Median survival by arm

median_table <- broom::tidy(km_fit) %>%
  group_by(strata) %>%
  summarise(
    median_survival_months = min(time[estimate <= 0.5], na.rm = TRUE),
    .groups = "drop"
  ) %>%
  mutate(
    arm = sub("arm=", "", strata),
    median_survival_months = ifelse(is.infinite(median_survival_months), NA, median_survival_months)
  ) %>%
  select(arm, median_survival_months)

knitr::kable(
  median_table,
  digits = 2,
  caption = "Approximate median survival by treatment arm."
)

Approximate median survival by treatment arm.
arm	median_survival_months
Control	10.94
Investigational	16.08

6 Log-rank test

The log-rank test compares survival curves between treatment arms without adjusting for additional covariates.

logrank_fit <- survdiff(Surv(time_months, event) ~ arm, data = trial_dat)

logrank_fit

## Call:
## survdiff(formula = Surv(time_months, event) ~ arm, data = trial_dat)
## 
##                       N Observed Expected (O-E)^2/E (O-E)^2/V
## arm=Control         193      119      104      2.23      4.06
## arm=Investigational 227      112      127      1.82      4.06
## 
##  Chisq= 4.1  on 1 degrees of freedom, p= 0.04

logrank_p <- 1 - pchisq(logrank_fit$chisq, df = length(logrank_fit$n) - 1)

logrank_p

## [1] 0.04395138

7 Cox proportional hazards model

The Cox model estimates the hazard ratio while accounting for censoring. The first model is unadjusted and includes only randomized treatment arm.

cox_unadjusted <- coxph(Surv(time_months, event) ~ arm, data = trial_dat)

summary(cox_unadjusted)

## Call:
## coxph(formula = Surv(time_months, event) ~ arm, data = trial_dat)
## 
##   n= 420, number of events= 231 
## 
##                       coef exp(coef) se(coef)      z Pr(>|z|)  
## armInvestigational -0.2646    0.7675   0.1317 -2.009   0.0446 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##                    exp(coef) exp(-coef) lower .95 upper .95
## armInvestigational    0.7675      1.303    0.5929    0.9936
## 
## Concordance= 0.531  (se = 0.017 )
## Likelihood ratio test= 4.03  on 1 df,   p=0.04
## Wald test            = 4.03  on 1 df,   p=0.04
## Score (logrank) test = 4.06  on 1 df,   p=0.04

broom::tidy(cox_unadjusted, exponentiate = TRUE, conf.int = TRUE) %>%
  transmute(
    term,
    hazard_ratio = estimate,
    conf_low = conf.low,
    conf_high = conf.high,
    p_value = p.value
  ) %>%
  knitr::kable(digits = 3, caption = "Unadjusted Cox proportional hazards model.")

Unadjusted Cox proportional hazards model.
term	hazard_ratio	conf_low	conf_high	p_value
armInvestigational	0.768	0.593	0.994	0.045

The adjusted model includes treatment arm, age, biomarker status, and region. In a real clinical trial, the adjustment set would be prespecified in the statistical analysis plan.

cox_adjusted <- coxph(
  Surv(time_months, event) ~ arm + age + biomarker + region,
  data = trial_dat
)

summary(cox_adjusted)

## Call:
## coxph(formula = Surv(time_months, event) ~ arm + age + biomarker + 
##     region, data = trial_dat)
## 
##   n= 420, number of events= 231 
## 
##                         coef exp(coef)  se(coef)      z Pr(>|z|)   
## armInvestigational -0.289884  0.748351  0.133087 -2.178  0.02939 * 
## age                 0.013792  1.013887  0.007328  1.882  0.05984 . 
## biomarkerHigh       0.371886  1.450468  0.135681  2.741  0.00613 **
## regionNorth         0.360998  1.434761  0.189294  1.907  0.05651 . 
## regionSouth         0.191874  1.211518  0.199138  0.964  0.33528   
## regionWest          0.048105  1.049280  0.193661  0.248  0.80383   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##                    exp(coef) exp(-coef) lower .95 upper .95
## armInvestigational    0.7484     1.3363    0.5765    0.9714
## age                   1.0139     0.9863    0.9994    1.0286
## biomarkerHigh         1.4505     0.6894    1.1118    1.8923
## regionNorth           1.4348     0.6970    0.9900    2.0793
## regionSouth           1.2115     0.8254    0.8200    1.7899
## regionWest            1.0493     0.9530    0.7179    1.5337
## 
## Concordance= 0.587  (se = 0.02 )
## Likelihood ratio test= 18.62  on 6 df,   p=0.005
## Wald test            = 18.9  on 6 df,   p=0.004
## Score (logrank) test = 18.99  on 6 df,   p=0.004

adjusted_hr_table <- broom::tidy(cox_adjusted, exponentiate = TRUE, conf.int = TRUE) %>%
  transmute(
    term,
    hazard_ratio = estimate,
    conf_low = conf.low,
    conf_high = conf.high,
    p_value = p.value
  )

knitr::kable(
  adjusted_hr_table,
  digits = 3,
  caption = "Adjusted Cox proportional hazards model."
)

Adjusted Cox proportional hazards model.
term	hazard_ratio	conf_low	conf_high	p_value
armInvestigational	0.748	0.577	0.971	0.029
age	1.014	0.999	1.029	0.060
biomarkerHigh	1.450	1.112	1.892	0.006
regionNorth	1.435	0.990	2.079	0.057
regionSouth	1.212	0.820	1.790	0.335
regionWest	1.049	0.718	1.534	0.804

8 Proportional-hazards diagnostics

The proportional-hazards assumption is evaluated using scaled Schoenfeld residuals.

ph_test <- cox.zph(cox_adjusted)

ph_test

##            chisq df    p
## arm       0.3421  1 0.56
## age       0.0002  1 0.99
## biomarker 0.0129  1 0.91
## region    1.7567  3 0.62
## GLOBAL    2.3163  6 0.89

plot(ph_test)

A non-significant test does not prove that the assumption is true, but it provides no strong evidence against proportional hazards in this simulated example. In real clinical-trial analysis, diagnostic plots and clinical plausibility should be evaluated together.

9 Adjusted survival curves

library(dplyr)
library(tidyr)
library(ggplot2)

new_patients <- data.frame(
  arm = factor(
    c("Control", "Investigational"),
    levels = c("Control", "Investigational")
  ),
  age = median(trial_dat$age, na.rm = TRUE),
  biomarker = factor("Low", levels = c("Low", "High")),
  region = factor("North", levels = levels(trial_dat$region))
)

adj_surv <- survfit(cox_adjusted, newdata = new_patients)

surv_mat <- adj_surv$surv

# If there are multiple newdata rows, survfit usually stores survival as a matrix
if (is.null(dim(surv_mat))) {
  surv_mat <- matrix(
    surv_mat,
    ncol = nrow(new_patients),
    byrow = FALSE
  )
}

colnames(surv_mat) <- as.character(new_patients$arm)

adj_df <- as.data.frame(surv_mat) %>%
  mutate(time = adj_surv$time) %>%
  pivot_longer(
    cols = -time,
    names_to = "arm",
    values_to = "estimate"
  ) %>%
  mutate(
    arm = factor(arm, levels = levels(new_patients$arm))
  )

ggplot(adj_df, aes(x = time, y = estimate, group = arm, linetype = arm)) +
  geom_step(linewidth = 0.9) +
  scale_y_continuous(
    labels = scales::percent_format(accuracy = 1),
    limits = c(0, 1)
  ) +
  labs(
    title = "Adjusted Survival Curves from Cox Model",
    subtitle = "Curves standardized to median age, biomarker-low status, and North region",
    x = "Months since randomization",
    y = "Adjusted survival probability",
    linetype = "Arm"
  ) +
  theme_minimal(base_size = 12)

10 Interpretation

In this simulated trial, the Kaplan-Meier curves visually compare time-to-event experience between randomized arms. The log-rank test provides an unadjusted comparison of the survival distributions. The Cox proportional hazards model estimates the treatment hazard ratio while accounting for censoring, and the adjusted model additionally accounts for age, biomarker status, and region.

The key points are:

Survival endpoints require special methods because some participants are censored.
Kaplan-Meier curves describe survival probability over time.
The log-rank test compares survival curves between groups.
The Cox model estimates hazard ratios and allows covariate adjustment.
The proportional-hazards assumption should be checked before relying on a standard Cox-model interpretation.

11 References

Kaplan, E. L., & Meier, P. (1958). Nonparametric estimation from incomplete observations. Journal of the American Statistical Association, 53(282), 457-481. https://doi.org/10.1080/01621459.1958.10501452

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

Therneau, T. M. (2024). A Package for Survival Analysis in R. R package version documentation. https://CRAN.R-project.org/package=survival

12 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] tidyr_1.3.1    scales_1.4.0   knitr_1.50     broom_1.0.8    ggplot2_3.5.2 
## [6] 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        cachem_1.1.0      
## [21] xfun_0.52          sass_0.4.10        cli_3.6.4          withr_3.0.2       
## [25] magrittr_2.0.3     digest_0.6.37      grid_4.5.0         rstudioapi_0.17.1 
## [29] lifecycle_1.0.4    vctrs_0.7.3        evaluate_1.0.3     glue_1.8.0        
## [33] farver_2.1.2       rmarkdown_2.29     purrr_1.2.2        tools_4.5.0       
## [37] pkgconfig_2.0.3    htmltools_0.5.8.1

Survival Analysis in a Simulated Randomized Clinical Trial

Md Ashraful Islam Bhuiya

2026-07-16