R simulation to calculate CP

library(survival)

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

library(stats)
set.seed(12345)

# Set parameters
alpha <- 0.05
beta <- 0.2
z_alpha <- qnorm(1 - alpha/2)
z_beta <- qnorm(1 - beta)
# Assumed Hazard rate for controlled group 1 which is always fixed
rate1 <- 0.7
# Assumed Hazard rate for investigational group 2
rate2 <- 0.3
# Information fraction at interim
t1 <- 0.5
# Actual Hazard rate for group 1
real_rate1 <- rate1
# Actual Hazard rate for group 2
real_rate2 <- 0.5
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)

# Assume interim observation is at half of the planned events are observed
interim_time <- log(1 / t1)/((rate1+ rate2) / 2)
# Compare with interim time
time <- ifelse(event_time <= interim_time, event_time, interim_time)
# Censoring 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
p1 <- real_rate1
p2 <- p * 2 - p1
# Calculate the difference in proportions
delta <- p1 - p2

# Calculate observed test statistic at interim
SE <- sqrt(p * (1 - p) * (1/n1 + 1/n2))
Z_t <-(p1-p2)/SE
# Type I error rate
alpha_D <- 0.05
Z_alpha_D <- qnorm(1 -alpha_D/2)
# Calculate B-value
B_t <- Z_t * sqrt(t1)
# Estimate drift parameter
delta_hat_t <- B_t / t1
# Calculate conditional power
Z_CP_t <- (delta_hat_t - Z_alpha_D) / sqrt(1 - t1)
cp <- pnorm(Z_CP_t)

# 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): 110.614

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

## Proportion of Zones: Unfavorable, Promissing, Favorable: 0.661 0.247 0.092