I did a CP/PZ simulation based on an unfavorable case, where the initial assumptions of hazard rates were 0.9 and 0.3 for two groups respectively, while the actual rates were 0.3 for both groups. In that case, I would expect almost all the simulations to turn out to be in the unfavorable zone. However, from Liu and Lim (2017) method, all 100 simulations turn out to be in either promising zone or favorable zone.
I think the major reason is because this paper used z1 (log-rank test statistic at interim) or all the calculations. When the interim is unblinded, it’s fine, z1 is based on observed real data. However, when the method is applied to blinded interim, with only pooled hazard rate observed, and applying the initial assumptions of either difference of hazard rate or hazard ratio, we assume there is still difference between groups. The logrank test is just testing the difference of the two survival curves. Therefore, the method of using log-rank test statistics at interim is not good for blinded interim.
# Load required libraries
library(survival)## Warning: package 'survival' was built under R version 4.2.3library(gsDesign)## Warning: package 'gsDesign' was built under R version 4.2.3set.seed(12345)
# Set parameters
k <- 2 # Number of analyses (including final)
alpha <- 0.025 # One-sided Type I error rate
beta <- 0.1 # Type II error rate (1 - power)
shape <- 0.15 # Parameter for Wang & Tsiatis boundaries (rho = 0.5 is equivalent to Pocock's boundary)
test.type <- 1 # Two-sided symmetric
rate1 <- 0.9 # Assumed Hazard rate for group 1
rate2 <- 0.3 # Assumed Hazard rate for group 2
t1 <- 0.5 # Information fraction at interim
real_rate1 <- 0.3 # Actual Hazard rate for group 1
real_rate2 <- 0.3 # Actual Hazard rate for group 2
num_simulations <- 100
# Calculate Wang & Tsiatis upper boundaries
calc_WT_boundary <- function(k, alpha, beta, shape, test.type) {
gs_result <- gsDesign(k = k, alpha = alpha, beta = beta, sfupar = shape, test.type = test.type, sfu = "WT")
return(gs_result$upper$bound)
}
# Calculate the WT upper boundaries
upper_boundaries <- calc_WT_boundary(k, alpha, beta, shape, test.type)
# Extract the specific upper boundary value (for the final analysis)
b <- upper_boundaries[2]
# Calculate initial d (number of events needed)
calc_d <- function(alpha, beta, rate1, rate2) {
z_alpha <- qnorm(1 - alpha)
z_beta <- qnorm(1 - beta)
log_hazard_ratio <- log(rate1 / rate2)
d <- 2 * ((z_alpha + z_beta)^2) / (log_hazard_ratio^2)
d1 <- t1 * d
d2 <- d - d1
return(list(d = d, d1 = d1, d2 = d2))
}
# Calculate initial d and d1
initial_d <- calc_d(alpha, beta, rate1, rate2)
# Extract initial d and d1, and dmax for increase
d <- ceiling(initial_d$d)
d1 <- ceiling(initial_d$d1)
d2 <- ceiling(initial_d$d2)
dmax <- 2 * ceiling(initial_d$d)
# Get number of subjects needed at initial planning (up round to next even integer)
n <- ceiling(d / ((rate1 + rate2) / 2) / 2) * 2
# Simulate z1 for num_simulations times at blinded interim
z1_values <- numeric(num_simulations)
for (i in 1:num_simulations) {
# Generate survival times
event_time1 <- rexp(n / 2, rate = real_rate1)
event_time2 <- rexp(n / 2, 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)
# Use full survival data to calculate interim time
cum_events <- cumsum(fit_full$n.event)
# Find the interim time at which half of the events is reached
interim_time <- fit_full$time[cum_events >= (d / 2)][1]
# 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
est.p <- 1 - min(a.pooled$surv)
# Calculate event rate in each treatment at interim, using assumed hazard ratio
p1 <- est.p * rate1 / (rate1 + rate2)
p2 <- est.p * rate2 / (rate1 + rate2)
# Perform simulation with calculated p1 and p2
simu_event_time1 <- rexp(n / 2, rate = p1)
simu_event_time2 <- rexp(n / 2, rate = p2)
simu_event_time <- c(simu_event_time1, simu_event_time2)
simu_time <- ifelse(simu_event_time <= interim_time, simu_event_time, interim_time)
simu_cens <- ifelse(simu_event_time <= interim_time, 1, 0)
simu_group <- c(rep(1, n / 2), rep(2, n / 2))
simu_fit <- survfit(Surv(simu_time, simu_cens) ~ simu_group)
simu_logrank <- survdiff(Surv(simu_time, simu_cens) ~ simu_group)
z1_values[i] <- sqrt(simu_logrank$chisq)
}## Warning in pchisq(chi, df, lower.tail = FALSE): NaNs produced## Warning in pchisq(chi, df, lower.tail = FALSE): NaNs produced# Calculate new number of events needed in stage two d2*
calc_d2_star <- function(beta, d1, z1, b, d, dmax) {
z_beta <- qnorm(1 - beta)
d2_pound <- d1 / (z1^2) * ((b * sqrt(d) - z1 * sqrt(d1)) / sqrt(d - d1) + z_beta)^2
d2_star <- ifelse(d2_pound <= dmax - d1, d2_pound, dmax - d1)
return(d2_star)
}
# Calculate the adjusted critical value b_star
calc_b_star <- function(d1, d2_star, d, b, z1) {
b_star <- 1 / sqrt(d1 + d2_star) * ((sqrt(d2_star / (d - d1)) * (b * sqrt(d) - z1 * sqrt(d1)) + z1 * sqrt(d1)))
return(b_star)
}
# Calculate the adjusted CP and original CP
calc_cp <- function(d1, d2_star, d, b_star, b, z1) {
val_new <- (b_star * sqrt((d1 + d2_star) / 4) - z1 * sqrt(d1 / 4) - sqrt(d2_star / 4) * z1 / sqrt(d1 / 4)) / sqrt(d2_star / 4)
cp_star <- 1 - pnorm(val_new)
val_org <- (b * sqrt((d1 + d2) / 4) - z1 * sqrt(d1 / 4) - sqrt(d2 / 4) * z1 / sqrt(d1 / 4)) / sqrt(d2 / 4)
cp_org <- 1 - pnorm(val_org)
return(list(cp_star = cp_star, cp_org = cp_org))
}
# Store results
results <- data.frame()
for (i in 1:num_simulations) {
z1 <- z1_values[i]
d2_star_value <- ceiling(calc_d2_star(beta, d1, z1, b, d, dmax))
d_star_value <- d2_star_value + d1
b_star_value <- calc_b_star(d1, d2_star_value, d, b, z1)
cp_values <- calc_cp(d1, d2_star_value, d, b_star_value, b, z1)
results <- rbind(results, data.frame(z1 = z1, d = d, d_star = d_star_value, b_star = b_star_value, cp_star = cp_values$cp_star, cp_org = cp_values$cp_org))
}
# Define Promising zone
pz <- results[results$b_star <= b, ]
cp_star_min_i <- ifelse(nrow(pz) > 0, min(pz$cp_org, na.rm = TRUE), 0)
cp_star_min <- min(cp_star_min_i, 0)
cp_star_max <- 1 - beta
# Define three zones
unfavorable_zone <- results[results$cp_org < cp_star_min, ]
if(nrow(unfavorable_zone) > 0) {
unfavorable_zone$zone <- 1
unfavorable_zone$d_star <- d
unfavorable_zone$cp_star <- unfavorable_zone$cp_org
}
promising_zone <- results[results$cp_org >= cp_star_min & results$cp_org <= cp_star_max, ]
if(nrow(promising_zone) > 0) {
promising_zone$zone <- 2
}
favorable_zone <- results[results$cp_org > cp_star_max, ]
if(nrow(favorable_zone) > 0) {
favorable_zone$zone <- 3
favorable_zone$d_star <- d
favorable_zone$cp_star <- favorable_zone$cp_org
}
results_final <- rbind(unfavorable_zone, promising_zone, favorable_zone)
results_final$sample_ratio <- results_final$d_star / results_final$d
results_final## z1 d d_star b_star cp_star cp_org zone sample_ratio
## 1 0.83207229 18 36 2.152383 0.073494980 0.073494980 2 2.0000000
## 2 0.55755214 18 36 2.252864 0.028209862 0.028209862 2 2.0000000
## 3 0.95748366 18 36 2.106479 0.107260960 0.107260960 2 2.0000000
## 4 0.80390747 18 36 2.162692 0.067172268 0.067172268 2 2.0000000
## 5 2.18584005 18 17 2.037147 0.789888290 0.789888290 2 0.9444444
## 6 1.47672763 18 36 1.916423 0.353525235 0.353525235 2 2.0000000
## 7 1.84211492 18 23 1.928547 0.592180826 0.592180826 2 1.2777778
## 8 1.50567948 18 36 1.905826 0.371619446 0.371619446 2 2.0000000
## 9 1.79508190 18 25 1.910610 0.561498727 0.561498727 2 1.3888889
## 10 0.74784183 18 36 2.183214 0.055846434 0.055846434 2 2.0000000
## 11 1.07250198 18 36 2.064380 0.146967184 0.146967184 2 2.0000000
## 12 1.03609160 18 36 2.077707 0.133453473 0.133453473 2 2.0000000
## 13 1.43902439 18 36 1.930223 0.330454296 0.330454296 2 2.0000000
## 14 1.00000000 18 36 2.090917 0.120926779 0.120926779 2 2.0000000
## 15 1.29467446 18 36 1.983059 0.248492179 0.248492179 2 2.0000000
## 16 1.07250198 18 36 2.064380 0.146967184 0.146967184 2 2.0000000
## 17 0.58009615 18 36 2.244613 0.030727599 0.030727599 2 2.0000000
## 18 0.58009615 18 36 2.244613 0.030727599 0.030727599 2 2.0000000
## 19 1.43902439 18 36 1.930223 0.330454296 0.330454296 2 2.0000000
## 20 0.89417934 18 36 2.129650 0.089032934 0.089032934 2 2.0000000
## 21 0.63916874 18 36 2.222991 0.038217881 0.038217881 2 2.0000000
## 22 0.83207229 18 36 2.152383 0.073494980 0.073494980 2 2.0000000
## 23 1.20849998 18 36 2.014601 0.205291981 0.205291981 2 2.0000000
## 24 1.77589072 18 25 1.914448 0.548861252 0.548861252 2 1.3888889
## 25 0.42342758 18 36 2.301957 0.016531377 0.016531377 2 2.0000000
## 26 1.00000000 18 36 2.090917 0.120926779 0.120926779 2 2.0000000
## 27 0.55755214 18 36 2.252864 0.028209862 0.028209862 2 2.0000000
## 28 1.00159780 18 36 2.090332 0.121463217 0.121463217 2 2.0000000
## 29 2.11344043 18 18 2.006085 0.753444031 0.753444031 2 1.0000000
## 30 0.58009615 18 36 2.244613 0.030727599 0.030727599 2 2.0000000
## 31 0.86232340 18 36 2.141310 0.080782507 0.080782507 2 2.0000000
## 32 1.43902439 18 36 1.930223 0.330454296 0.330454296 2 2.0000000
## 33 1.79508190 18 25 1.910610 0.561498727 0.561498727 2 1.3888889
## 34 0.86232340 18 36 2.141310 0.080782507 0.080782507 2 2.0000000
## 35 1.20370916 18 36 2.016355 0.203028869 0.203028869 2 2.0000000
## 36 0.00000000 18 36 2.456943 0.002276745 0.002276745 2 2.0000000
## 37 0.82356203 18 36 2.155498 0.071538294 0.071538294 2 2.0000000
## 38 0.48316675 18 36 2.280091 0.021089226 0.021089226 2 2.0000000
## 39 1.79508190 18 25 1.910610 0.561498727 0.561498727 2 1.3888889
## 40 1.24343737 18 36 2.001813 0.222243738 0.222243738 2 2.0000000
## 41 2.11344043 18 18 2.006085 0.753444031 0.753444031 2 1.0000000
## 43 1.79508190 18 25 1.910610 0.561498727 0.561498727 2 1.3888889
## 44 0.48316675 18 36 2.280091 0.021089226 0.021089226 2 2.0000000
## 45 0.02439024 18 36 2.448015 0.002583909 0.002583909 2 2.0000000
## 46 1.37835320 18 36 1.952430 0.294675167 0.294675167 2 2.0000000
## 47 0.04881744 18 36 2.439074 0.002928597 0.002928597 2 2.0000000
## 48 0.03577135 18 36 2.443849 0.002739666 0.002739666 2 2.0000000
## 49 1.03609160 18 36 2.077707 0.133453473 0.133453473 2 2.0000000
## 50 1.03609160 18 36 2.077707 0.133453473 0.133453473 2 2.0000000
## 51 2.41215184 18 14 2.187939 0.881639020 0.881639020 2 0.7777778
## 52 0.00000000 18 36 2.456943 0.002276745 0.002276745 2 2.0000000
## 53 1.43902439 18 36 1.930223 0.330454296 0.330454296 2 2.0000000
## 54 1.43902439 18 36 1.930223 0.330454296 0.330454296 2 2.0000000
## 55 0.04881744 18 36 2.439074 0.002928597 0.002928597 2 2.0000000
## 57 0.02439024 18 36 2.448015 0.002583909 0.002583909 2 2.0000000
## 58 1.43902439 18 36 1.930223 0.330454296 0.330454296 2 2.0000000
## 59 2.41215184 18 14 2.187939 0.881639020 0.881639020 2 0.7777778
## 60 1.07250198 18 36 2.064380 0.146967184 0.146967184 2 2.0000000
## 61 0.93066773 18 36 2.116295 0.099235462 0.099235462 2 2.0000000
## 62 1.26435129 18 36 1.994158 0.232762062 0.232762062 2 2.0000000
## 63 1.40922170 18 36 1.941132 0.312657979 0.312657979 2 2.0000000
## 64 0.55755214 18 36 2.252864 0.028209862 0.028209862 2 2.0000000
## 65 1.79508190 18 25 1.910610 0.561498727 0.561498727 2 1.3888889
## 66 0.02439024 18 36 2.448015 0.002583909 0.002583909 2 2.0000000
## 67 1.39083689 18 36 1.947861 0.301890285 0.301890285 2 2.0000000
## 68 1.44192976 18 36 1.929160 0.332210745 0.332210745 2 2.0000000
## 69 0.58009615 18 36 2.244613 0.030727599 0.030727599 2 2.0000000
## 70 1.67178985 18 29 1.899013 0.479775659 0.479775659 2 1.6111111
## 71 1.29708956 18 36 1.982175 0.249768916 0.249768916 2 2.0000000
## 72 1.00159780 18 36 2.090332 0.121463217 0.121463217 2 2.0000000
## 73 1.11110266 18 36 2.050251 0.162264080 0.162264080 2 2.0000000
## 74 1.51391656 18 36 1.902811 0.376821741 0.376821741 2 2.0000000
## 75 2.15136921 18 17 2.035713 0.772945499 0.772945499 2 0.9444444
## 76 1.43902439 18 36 1.930223 0.330454296 0.330454296 2 2.0000000
## 77 0.48316675 18 36 2.280091 0.021089226 0.021089226 2 2.0000000
## 78 2.11344043 18 18 2.006085 0.753444031 0.753444031 2 1.0000000
## 79 1.87898132 18 22 1.938261 0.615851218 0.615851218 2 1.2222222
## 80 0.59737730 18 36 2.238287 0.032781172 0.032781172 2 2.0000000
## 81 2.18584005 18 17 2.037147 0.789888290 0.789888290 2 0.9444444
## 82 2.41215184 18 14 2.187939 0.881639020 0.881639020 2 0.7777778
## 83 1.02173245 18 36 2.082963 0.128366734 0.128366734 2 2.0000000
## 84 1.07250198 18 36 2.064380 0.146967184 0.146967184 2 2.0000000
## 85 0.04206097 18 36 2.441547 0.002829333 0.002829333 2 2.0000000
## 86 1.47672763 18 36 1.916423 0.353525235 0.353525235 2 2.0000000
## 87 1.03609160 18 36 2.077707 0.133453473 0.133453473 2 2.0000000
## 88 1.37835320 18 36 1.952430 0.294675167 0.294675167 2 2.0000000
## 89 1.40922170 18 36 1.941132 0.312657979 0.312657979 2 2.0000000
## 90 1.55973808 18 34 1.897750 0.406146253 0.406146253 2 1.8888889
## 92 0.03577135 18 36 2.443849 0.002739666 0.002739666 2 2.0000000
## 93 2.11344043 18 18 2.006085 0.753444031 0.753444031 2 1.0000000
## 94 0.63916874 18 36 2.222991 0.038217881 0.038217881 2 2.0000000
## 95 2.11344043 18 18 2.006085 0.753444031 0.753444031 2 1.0000000
## 96 1.91873006 18 21 1.950275 0.640894350 0.640894350 2 1.1666667
## 97 0.63916874 18 36 2.222991 0.038217881 0.038217881 2 2.0000000
## 98 1.44192976 18 36 1.929160 0.332210745 0.332210745 2 2.0000000
## 99 1.43902439 18 36 1.930223 0.330454296 0.330454296 2 2.0000000
## 100 1.07250198 18 36 2.064380 0.146967184 0.146967184 2 2.0000000
## 42 3.27522012 18 18 2.968580 0.995624958 0.995624958 3 1.0000000
## 56 2.84910453 18 18 2.571965 0.972028185 0.972028185 3 1.0000000
## 91 3.14451287 18 18 2.885913 0.991887658 0.991887658 3 1.0000000