# Specify the package name
packages <- c("gsDesign", "stats", "survival")
# Function to check, install and call the packages
install_if_needed <- function(packages) {
for (package in packages) {
if (!require(package, character.only = TRUE)) {
install.packages(package, dependencies = TRUE)
library(package, character.only = TRUE)
}
}
}
# Execute the function
install_if_needed(packages)
## Loading required package: gsDesign
## Warning: package 'gsDesign' was built under R version 4.2.3
## Loading required package: survival
## Warning: package 'survival' was built under R version 4.2.3
# 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 <- 2 # Two-sided symmetric
rate1 <- 0.02 # Hazard rate for group 1
rate2 <- 0.01 # Hazard rate for group 2
t1 <- 0.5 # Information fraction
z1 <- seq(0.5, 3.00, by = 0.01)
# 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]
# Get b=2.006085
# 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*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 d=88, d1=d2=44, dmax=176
# Use literature d assumptions to match literature's results
d <- 90
d1 <- 45
d2 <- 45
dmax <- 180
# Step (i) and (ii): 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 new d2* needed and round to next integer
calc_d2_star <- calc_d2_star(beta, d1, z1, b, d, dmax)
d2_star <- ceiling(calc_d2_star)
# Step (iii): 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 critical value b_star
b_star <- calc_b_star(d1, d2_star, d, b, z1)
# Step (iv) and (v): Calculate the new adjusted CP based on d2*, and original CP based on d2
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)*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)*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))
}
# Calculate the adjusted CP and original CP
calc_cp <- calc_cp(d1, d2_star, d, b_star, b, z1)
cp_star <- calc_cp$cp_star
cp_org <- calc_cp$cp_org
# Step (vi): Iterate steps (i)–(v) for z1 ∈ [0.5, 3.00] by increment of 0.01 (already done by parameter settings)
# Combine all results from Part 1 CP and critical value calculation
results <- data.frame(z1 = z1, d2_star = d2_star, b_star = b_star, cp_star = cp_star, cp_org = cp_org)
# Part 2: Identify the promising zone
# Step (vii) and (viii): Plotting b*(z1, d2*) versus z1, and b as reference line
# Adjust margins to make space for the secondary y-axis label
par(mar = c(4, 5, 5, 4) + 0.1)
plot(results$z1, results$b_star, type = 'l', lty = 'longdash', col = 4, yaxt="n", ylim = c(1, 3.5), ylab = "Critical Value", xlab = "z1 (Log-rank test statistic at interim)")
abline(h = b, col = 3, lty = "dashed")
# Add the first y-axis on the left
axis(side = 2, at=seq(1, 3.5, by=0.5))
# Step (x) and (xi): Define Promising zone
pz <- results[results$b_star <= b, ]
cp_star_min <- min(pz$cp_org)
cp_star_max <- 1 - beta
# Step (xii): Define three zones
unfavorable_zone <- results[results$cp_org < cp_star_min, ]
promising_zone <- results[results$cp_org >= cp_star_min & results$cp_org <= cp_star_max, ]
favorable_zone <- results[results$cp_org > cp_star_max, ]
# Step (xiii): Set d2_star = d2 and cp_star = cp_org when z1 is located in either unfavorable or favorable zones
unfavorable_zone$d2_star <- d2
unfavorable_zone$cp_star <- unfavorable_zone$cp_org
favorable_zone$d2_star <- d2
favorable_zone$cp_star <- favorable_zone$cp_org
results <- rbind(unfavorable_zone, promising_zone, favorable_zone)
# Step (ix): Plotting the curve of original CP against z1, using original d2
# Allows the addition of a new plot on top of the existing one.
par(new = TRUE)
plot(results$z1, results$cp_org, type = 'l', lty = "dotted", col = 6, yaxt="n", ylim = c(0, 1), ylab = " ", xlab = " ")
# Add the secondary y-axis on the right
axis(side = 4, at=seq(0, 1.0, by=0.1))
mtext("Conditional Power", side = 4, line=2.5)
# Step (xiv): Plotting the curve of new CP against z1, using adjusted d2*
# Allows the addition of a new plot on top of the existing one.
par(new = TRUE)
plot(results$z1, results$cp_star, type = 'l', lty = "dotdash", col = 8, yaxt="n", ylim = c(0, 1), ylab = " ", xlab = " ")
axis(side = 4)
legend("topleft", legend = c("b*", paste("b =", round(b, 6)), "CP", "CP*"), col = c(4, 3, 6, 8), lty = c("longdash", "dashed", "dotted", "dotdash"), bty = "n", text.font = 1, cex=0.8)
# Print out some results
example <- subset(results, z1 >= 1.23 & z1 <= 2.08)
example
## z1 d2_star b_star cp_star cp_org
## 74 1.23 45 2.006731 0.3530745 0.3530745
## 75 1.24 135 2.003071 0.7090837 0.3605335
## 76 1.25 135 1.999411 0.7183782 0.3680460
## 77 1.26 135 1.995751 0.7275271 0.3756093
## 78 1.27 135 1.992090 0.7365259 0.3832206
## 79 1.28 135 1.988430 0.7453705 0.3908774
## 80 1.29 135 1.984770 0.7540571 0.3985766
## 81 1.30 135 1.981110 0.7625821 0.4063156
## 82 1.31 135 1.977449 0.7709422 0.4140913
## 83 1.32 135 1.973789 0.7791347 0.4219008
## 84 1.33 135 1.970129 0.7871567 0.4297412
## 85 1.34 135 1.966469 0.7950061 0.4376094
## 86 1.35 135 1.962808 0.8026807 0.4455023
## 87 1.36 135 1.959148 0.8101790 0.4534169
## 88 1.37 135 1.955488 0.8174994 0.4613501
## 89 1.38 135 1.951828 0.8246409 0.4692986
## 90 1.39 135 1.948167 0.8316026 0.4772594
## 91 1.40 135 1.944507 0.8383840 0.4852293
## 92 1.41 135 1.940847 0.8449849 0.4932051
## 93 1.42 135 1.937187 0.8514052 0.5011836
## 94 1.43 135 1.933526 0.8576452 0.5091617
## 95 1.44 135 1.929866 0.8637056 0.5171361
## 96 1.45 135 1.926206 0.8695869 0.5251036
## 97 1.46 135 1.922546 0.8752904 0.5330611
## 98 1.47 135 1.918885 0.8808172 0.5410054
## 99 1.48 135 1.915225 0.8861688 0.5489334
## 100 1.49 135 1.911565 0.8913469 0.5568419
## 101 1.50 135 1.907905 0.8963534 0.5647278
## 102 1.51 135 1.904244 0.9011903 0.5725880
## 103 1.52 132 1.903771 0.9008251 0.5804196
## 104 1.53 129 1.903478 0.9003320 0.5882194
## 105 1.54 127 1.902226 0.9014895 0.5959845
## 106 1.55 124 1.902270 0.9007691 0.6037120
## 107 1.56 122 1.901291 0.9017482 0.6113989
## 108 1.57 119 1.901695 0.9007892 0.6190425
## 109 1.58 117 1.901006 0.9015858 0.6266399
## 110 1.59 114 1.901794 0.9003734 0.6341883
## 111 1.60 112 1.901414 0.9009809 0.6416852
## 112 1.61 110 1.901175 0.9015019 0.6491279
## 113 1.62 108 1.901081 0.9019357 0.6565137
## 114 1.63 105 1.902664 0.9002247 0.6638403
## 115 1.64 103 1.902919 0.9004527 0.6711052
## 116 1.65 101 1.903334 0.9005880 0.6783060
## 117 1.66 99 1.903912 0.9006292 0.6854404
## 118 1.67 97 1.904660 0.9005740 0.6925062
## 119 1.68 95 1.905582 0.9004203 0.6995014
## 120 1.69 93 1.906683 0.9001655 0.7064237
## 121 1.70 92 1.906071 0.9021148 0.7132713
## 122 1.71 90 1.907488 0.9016953 0.7200422
## 123 1.72 88 1.909100 0.9011680 0.7267347
## 124 1.73 86 1.910913 0.9005287 0.7333469
## 125 1.74 85 1.910795 0.9022309 0.7398772
## 126 1.75 83 1.912962 0.9014135 0.7463241
## 127 1.76 81 1.915349 0.9004736 0.7526861
## 128 1.77 80 1.915626 0.9019888 0.7589618
## 129 1.78 78 1.918398 0.9008562 0.7651499
## 130 1.79 77 1.918938 0.9022460 0.7712492
## 131 1.80 75 1.922117 0.9009105 0.7772584
## 132 1.81 74 1.922933 0.9021736 0.7831767
## 133 1.82 72 1.926542 0.9006222 0.7890031
## 134 1.83 71 1.927650 0.9017561 0.7947366
## 135 1.84 70 1.928874 0.9028306 0.8003765
## 136 1.85 68 1.933131 0.9009740 0.8059222
## 137 1.86 67 1.934670 0.9019148 0.8113729
## 138 1.87 66 1.936335 0.9027961 0.8167283
## 139 1.88 64 1.941292 0.9006035 0.8219878
## 140 1.89 63 1.943298 0.9013436 0.8271511
## 141 1.90 62 1.945441 0.9020228 0.8322179
## 142 1.91 61 1.947722 0.9026410 0.8371881
## 143 1.92 60 1.950147 0.9031981 0.8420616
## 144 1.93 58 1.956331 0.9004432 0.8468383
## 145 1.94 57 1.959142 0.9008421 0.8515182
## 146 1.95 56 1.962108 0.9011763 0.8561015
## 147 1.96 55 1.965233 0.9014449 0.8605885
## 148 1.97 54 1.968520 0.9016471 0.8649792
## 149 1.98 53 1.971973 0.9017818 0.8692742
## 150 1.99 52 1.975596 0.9018479 0.8734738
## 151 2.00 51 1.979394 0.9018439 0.8775785
## 152 2.01 50 1.983370 0.9017683 0.8815889
## 153 2.02 49 1.987528 0.9016193 0.8855054
## 154 2.03 48 1.991875 0.9013951 0.8893289
## 155 2.04 47 1.996413 0.9010935 0.8930600
## 156 2.05 46 2.001148 0.9007121 0.8966995
## 157 2.06 45 2.006085 0.9002482 0.9002482
## 158 2.07 45 2.006085 0.9037070 0.9037070
## 159 2.08 45 2.011309 0.9070768 0.9070768