R simulation to Calculate Binomial Rate at Blinded Interim

library(survival)

## Warning: package 'survival' was built under R version 4.2.3

set.seed(12345)

# Set parameters
alpha <- 0.05
beta <- 0.2
# Assumed Hazard rate for group 1
rate1 <- 0.7
# Assumed Hazard rate for group 2 which is always fixed  
rate2 <- 0.3
# Information fraction at interim
t1 <- 0.5
# Actual Hazard rate for group 1
real_rate1 <- 0.5
# Actual Hazard rate for group 2
real_rate2 <-  rate2
n <- 60
n1 <- n/2
n2 <- n/2
num_simulations <- 1000
results <- data.frame()

for (i in 1:num_simulations) {
# Generate survival times
event_time1 <- rexp(n1, rate = real_rate1)
event_time2 <- rexp(n2, rate = real_rate2)
event_time <- c(event_time1, event_time2)

# Generate full survival data for calculation of interim time
cens_full <- rep(1, n)
fit_full <- survfit(Surv(event_time, cens_full) ~ 1)

# Assume follow-up length is half of the events are observed, and interim observation is at half of the planned events observed
interim_time <-  mean(event_time) / 2 * t1

# Compare with interim time
time <- ifelse(event_time <= interim_time, event_time, interim_time)
# Event indicator: 1 = event, 0 = censored
cens <- ifelse(event_time <= interim_time, 1, 0)

# Pooled arms for blinded interim
a.pooled <- survfit(Surv(time, cens) ~ 1)

# Get the event rate at the last event time point at interim
p <- 1 - min(a.pooled$surv)

# Calculate event rate in each treatment at interim, with control group fixed
p2 <- real_rate2
p1 <- p * 2 - p2

# Calculate the standard error of the difference in proportions
sigma <- sqrt(p * (1 - p) * (1 / n1 + 1 / n2))

# Calculate the difference in proportions
delta <- p1 - p2

# Determine the critical value from the normal distribution
z_alpha <- qnorm(1 - alpha / 2)
z_beta <- qnorm(1 - beta)

# Calculate the test statistic
z <- delta / sigma

# Calculate the power
cp <- 1 - pnorm(z_alpha - z) + pnorm(-z_alpha - z)

# Define bounds for stopping for futility and continuing without resizing
CL <- 0.2
CU <- 0.8

# Decision rules
if (cp <= CL) {
   cat="Unfavorable"
   n_new=n
} else if (cp >= CU) {
   cat="Favorable"
   n_new=n
} else {
   cat="Promising"
   # Calculate the required sample size  
   n_new <- ceiling(((z_alpha + z_beta)^2 * 2 * p * (1 - p)) / delta^2) * 2
}

results <- rbind(results, data.frame(cp = cp, cat=cat, n_new=n_new))

}

est_n <- mean(results$n_new)

# Calculate the frequency of each category 
cat_count <- table(results$cat) 
# Calculate the proportion of each category 
prop_cat <- prop.table(cat_count)

cat("Estimation of Sample Size (Mean):", est_n, "\n")

## Estimation of Sample Size (Mean): 126.95

cat("Proportion of Zones: Unfavorable, Promissing, Favorable:", prop_cat, "\n")

## Proportion of Zones: Unfavorable, Promissing, Favorable: 0.07 0.55 0.38