I am looking to simulate a time to event study in R from the exponential distribution. I want to explore different scenarios to test the power of the log-rank test when cure exists. I want to look at different cut-offs for event driven analysis and assess the power also be able to look at different calendar time analysis
I want the study to have N=400, 200 per arm, enrolment rate of 20 patients per month from uniform distribution, random censoring. Also look at different cure proportions between control and treatment to see how this effects the power of log-rank
I am unsure how to set the initial code up to run the different scenarios
I have tried to simulate a standard study in R first but I am not sure I have even done this correctly.
# Load required libraries
library(survival)
library(survminer)
# Set seed for reproducibility
set.seed(42)
# Parameters
n_months <- 60 # Total study duration in months
n_patients_per_month <- 20 # Number of patients enrolled per month
total_patients <- n_months * n_patients_per_month # Total number of patients
median_pfs_control <- 15 # Median PFS in control arm in months
hazard_ratio <- 0.8 # Hazard ratio
# Generate time to event for control and treatment arms
# Control arm
pfs_control <- rexp(total_patients, rate = 1/median_pfs_control)
# Treatment arm
pfs_treatment <- rexp(total_patients, rate = 1/(median_pfs_control / hazard_ratio))
# Random censoring variable
censoring_control <- runif(total_patients, min = 0, max = median_pfs_control * 2)
censoring_treatment <- runif(total_patients, min = 0, max = median_pfs_control * 2 / hazard_ratio)
# Apply censoring
observed_control <- pmin(pfs_control, censoring_control)
observed_treatment <- pmin(pfs_treatment, censoring_treatment)
# Create survival objects
surv_control <- Surv(observed_control)
surv_treatment <- Surv(observed_treatment)
# Fit Kaplan-Meier models
fit_control <- survfit(surv_control ~ 1)
fit_treatment <- survfit(surv_treatment ~ 1)
# Plot Kaplan-Meier curves
plot(fit_control, col = "blue", lwd = 2, main = "Kaplan-Meier Curves", xlab = "Time (months)", ylab = "Survival Probability")
lines(fit_treatment, col = "red", lwd = 2)
legend("topright", legend = c("Control Arm", "Treatment Arm"), col = c("blue", "red"), lwd = 2)
# Median survival times
median_pfs_control_estimate <- median(fit_control$time)
median_pfs_treatment_estimate <- median(fit_treatment$time)
cat("Median PFS in Control Arm:", median_pfs_control_estimate, "months\n")
cat("Median PFS in Treatment Arm:", median_pfs_treatment_estimate, "months\n")
I want to ensure this code is correct and add the other part into this