library(survival)
## Warning: package 'survival' was built under R version 4.2.3
library(stats)
library(dplyr)
## Warning: package 'dplyr' was built under R version 4.2.3
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
set.seed(12345)
simu <- function(alpha, beta, rate1, rate2, t1, real_rate1, real_rate2, rec, CL, CU, num_simulations){
z_alpha <- qnorm(1 - alpha/2)
z_beta <- qnorm(1 - beta)
rate <- (rate1+rate2)/2
diff <- rate1-rate2
# n <- ceiling(((z_alpha + z_beta)^2 * 2 * rate * (1 - rate)) / diff^2) * 2
n <- ceiling (((z_alpha + z_beta)^2 * (rate1 * (1 - rate1) + rate2 * (1 - rate2))) / diff^2) * 2
n1 <- n/2
n2 <- n/2
# Calculate initial events needed
d <- ceiling((n1 * rate1) + (n2 * rate2))
# Actual Hazard rate for group 1 is fixed
real_rate1 <- rate1
final_results <- data.frame()
scen_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)
group1 <- rep(1, n1)
group2 <- rep(2, n2)
group <- c(group1, group2)
subj <-1:n
# Total follow-up time with d events occurred
follow_up <- log(1-d/n)/(- (rate1+ rate2) / 2)
# Assume recruitment time is rec rate of total follow up time (final time)
rec_length <- follow_up*rec
# Assume uniform distribution for recruitment
rec_time <- runif(n, 0, rec_length)
# Study time for each subject at event occurence is recruitment time plus event time
study_time <- rec_time + event_time
# Time to interim analysis (half of final events observed)
sorted_sequence <- sort(study_time)
interim_time <- sorted_sequence[d*t1]
all_data <- data.frame(study_time = study_time, event_time = event_time, rec_time = rec_time, group = group, interim_time=interim_time, subj=subj)
surv_data <- all_data %>% filter(rec_time<interim_time)
left_data <- all_data %>% filter(rec_time>=interim_time)
# Compare with interim time
surv_data$time <- ifelse(surv_data$study_time <= surv_data$interim_time, surv_data$event_time, surv_data$interim_time-surv_data$rec_time)
# Censoring indicator: 1 = event, 0 = censored
surv_data$cens <- ifelse(surv_data$study_time <= surv_data$interim_time, 1, 0)
# Pooled arms for blinded interim
a.pooled <- survfit(Surv(surv_data$time, surv_data$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)
# Decision rules
if (cp > CL & cp < CU) {
cat="Promising"
# Calculate the required sample size
n_star <- ceiling (((z_alpha + z_beta)^2 * (p1 * (1 - p1) + p2 * (1 - p2))) / delta^2) * 2
# Apply max N
n_new <- ifelse(n_star<=2*n, n_star, 2*n)
# add n
n_add <- n_new-n
add_event_time1 <- rexp(n_add/2, rate = real_rate1)
add_event_time2 <- rexp(n_add/2, rate = real_rate2)
add_event_time <- c(add_event_time1, add_event_time2)
add_group1 <- rep(1, n_add/2)
add_group2 <- rep(2, n_add/2)
add_group <- c(add_group1, add_group2)
add_rec_length <- rec_length*n_add/n
add_rec_time <- runif(n_add, 0, add_rec_length) + interim_time
add_study_time <- add_rec_time + add_event_time
add_subj <- (n+1):n_new
add_data <- data.frame(study_time = add_study_time, event_time = add_event_time, rec_time = add_rec_time, group = add_group, interim_time=interim_time, subj=add_subj)
comb_data <- rbind(all_data, add_data)
# Time to final analysis (d events observed)
final_sorted_sequence <- sort(comb_data$study_time)
final_time <- final_sorted_sequence[d]
final_data <- comb_data %>% subset(select=c(study_time, event_time, rec_time, group, interim_time, subj))
final_data$final_time <- final_time
# event time compare with final time
final_data$time <- ifelse(final_data$study_time <= final_data$final_time, final_data$event_time, final_data$final_time-final_data$rec_time)
# Censoring indicator: 1 = event, 0 = censored
final_data$cens <- ifelse(final_data$study_time <= final_data$final_time, 1, 0)
}
else {
n_new=n
final_sorted_sequence <- sort(all_data$study_time)
final_time <- final_sorted_sequence[d]
final_data <- all_data %>% subset(select=c(study_time, event_time, rec_time, group, interim_time, subj))
final_data$final_time <- final_time
# event time compare with final time
final_data$time <- ifelse(final_data$study_time <= final_data$final_time, final_data$event_time, final_data$final_time-final_data$rec_time)
# Censoring indicator: 1 = event, 0 = censored
final_data$cens <- ifelse(final_data$study_time <= final_data$final_time, 1, 0)
if (cp <= CL) {cat="Unfavorable"}
else if (cp >= CU) {cat="Favorable"}
}
# final and interim analysis results
final_table <- table(final_data$group, final_data$cens)
count_group1 <- final_table["1", "1"]
count_group2 <- final_table["2", "1"]
final_p1 <- count_group1/(n_new/2)
final_p2 <- count_group2/(n_new/2)
final_p <- (count_group1 + count_group2) / n_new
final_se <- sqrt(final_p * (1 - final_p) * (1 / (n_new/2) + 1 / (n_new/2)))
final_z <- (final_p1 - final_p2) / final_se
p_value <- 2 * (1 - pnorm(abs(final_z)))
hypo <- ifelse(p_value < alpha, 1, 0)
# Evaluate how much percent of subjects' time is beyond follow-up time at interim
eval_data <- merge(x=surv_data, y=final_data, by = "subj", all.x = TRUE)
eval <- sum(eval_data$time.x == eval_data$time.y, na.rm = TRUE)/n
final_results <- rbind(final_results, data.frame(cp = cp, cat=cat, eval=eval, n_new = n_new, p_value = p_value, hypo = hypo))
}
# interim analysis results printouts
est_n <- mean(final_results$n_new)
cp_mean <- mean(final_results$cp)
eval_mean <- mean(final_results$eval)
# Calculate the frequency and proportion of each category
cat_count <- table(final_results$cat)
prop_cat <- prop.table(cat_count)
re_est <- prop_cat["Promising"]
cat("Subjects Reached Final Endpoint at Interim", eval_mean, "\n")
cat("Estimation of Sample Size (Mean):", est_n, "\n")
cat("New Sample Size / Original Sample Size:", est_n/n, "\n")
cat("Estimation of Conditional Power:", cp_mean, "\n")
cat("Proportion of Zones: Favorable, Promissing, Unfavorable:", prop_cat, "\n")
# final analysis results printouts
sig <- sum(final_results$hypo)/num_simulations
nosig <- 1 - sig
cat("Probability of claim success:", sig, "\n")
cat("Proportion of not claim success:", nosig, "\n")
# scen_results <- rbind(scen_results, data.frame(real_rate1<-real_rate1, real_rate2<-real_rate2, t1<-t1, recru_vs_study_length<-rec, CL<-CL, CU<-CU, N_new<-est_n, N_new_vs_old <- est_n/n, CP<-cp_mean, Increase<-re_est, Final_at_Interim<-eval_mean, H1<-sig, H0<-nosig))
}
### under alternative hypothesis and type II error
simu(alpha=0.05, beta=0.2, rate1=0.7, rate2=0.3, t1=0.5, real_rate1=0.7, real_rate2=0.2, rec=1, CL=0.2, CU=0.8, num_simulations=1000)
## Subjects Reached Final Endpoint at Interim 0.2380952
## Estimation of Sample Size (Mean): 44.31
## New Sample Size / Original Sample Size: 1.055
## Estimation of Conditional Power: 0.8462952
## Proportion of Zones: Favorable, Promissing, Unfavorable: 0.823 0.055 0.122
## Probability of claim success: 0.804
## Proportion of not claim success: 0.196
### when real rate is same as assumed rate
simu(alpha=0.05, beta=0.2, rate1=0.7, rate2=0.3, t1=0.5, real_rate1=0.7, real_rate2=0.3, rec=1, CL=0.2, CU=0.8, num_simulations=1000)
## Subjects Reached Final Endpoint at Interim 0.2380952
## Estimation of Sample Size (Mean): 45.486
## New Sample Size / Original Sample Size: 1.083
## Estimation of Conditional Power: 0.813716
## Proportion of Zones: Favorable, Promissing, Unfavorable: 0.774 0.083 0.143
## Probability of claim success: 0.471
## Proportion of not claim success: 0.529
### under null hypothesis and type I error (lower boundary at 20%)
simu(alpha=0.05, beta=0.2, rate1=0.7, rate2=0.3, t1=0.5, real_rate1=0.7, real_rate2=0.7, rec=3, CL=0.2, CU=0.8, num_simulations=1000)
## Subjects Reached Final Endpoint at Interim 0.2380952
## Estimation of Sample Size (Mean): 47.67
## New Sample Size / Original Sample Size: 1.135
## Estimation of Conditional Power: 0.1245066
## Proportion of Zones: Favorable, Promissing, Unfavorable: 0.054 0.135 0.811
## Probability of claim success: 0.06
## Proportion of not claim success: 0.94
### under null hypothesis and type I error (lower boundary at 40%)
simu(alpha=0.05, beta=0.2, rate1=0.7, rate2=0.3, t1=0.5, real_rate1=0.7, real_rate2=0.7, rec=3, CL=0.4, CU=0.8, num_simulations=1000)
## Subjects Reached Final Endpoint at Interim 0.2380952
## Estimation of Sample Size (Mean): 45.066
## New Sample Size / Original Sample Size: 1.073
## Estimation of Conditional Power: 0.1322873
## Proportion of Zones: Favorable, Promissing, Unfavorable: 0.073 0.073 0.854
## Probability of claim success: 0.056
## Proportion of not claim success: 0.944