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_rate2, rec, CL, CU, num_simulations){
z_alpha <- qnorm(1 - alpha/2)
z_beta <- qnorm(1 - beta)
# Rate is the probability to develop disease event in 76 weeks (for example)
rate <- (rate1+rate2)/2
diff <- rate1-rate2
# Calculate total Z-score required
z_total <- z_alpha + z_beta
# Calculate required sample size
n <- ceiling((z_total^2 * (rate * (1 - rate) * 2)) / diff^2)*2
n1 <- n/2
n2 <- n/2
# Actual event rate for group 1 is fixed
real_rate1 <- rate1
final_results <- data.frame()
scen_results <- data.frame()
count_simu <- 0
for (i in 1:num_simulations) {
# Calculate hazard rate from event rate
hr1 <- -log(1-real_rate1) / 76
hr2 <- -log(1-real_rate2) / 76
# Generate survival times
event_time1 <- rexp(n1, rate = hr1)
event_time2 <- rexp(n2, rate = hr2)
event_time <- c(event_time1, event_time2)
group1 <- rep(1, n1)
group2 <- rep(2, n2)
group <- c(group1, group2)
subj <-1:n
# Follow-up time is 76 weeks
follow_up <- 76
# 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 occurrence is recruitment time plus event time
study_time <- rec_time + event_time
# Final time for each subject (if no event happened) is recruitment time plus follow-up time
final_time <- rec_time + follow_up
# Time to interim analysis (n*t1 subjects reach 76 weeks or have events)
sorted_interim <- sort(final_time)
interim_time <- sorted_interim[n*t1]
all_data <- data.frame(study_time = study_time, event_time = event_time, rec_time = rec_time, group = group, interim_time = interim_time, final_time = final_time, subj=subj)
surv_data <- all_data %>% filter(final_time <= interim_time)
# Compare with interim time for event indicator: 1 = event, 0 = censored
surv_data$cens <- ifelse(surv_data$study_time <= surv_data$interim_time, 1, 0)
n_int <- nrow(surv_data)
p <- sum(surv_data$cens)/n_int
# Calculate event rate in each treatment at interim, with control group fixed
if (p==1 | p==0 | p*2<=real_rate1 | p*2-real_rate1>=1) {
next
}
else{
p1 <- real_rate1
p2 <- p * 2 - p1
delta <- p1 - p2
# Calculate interim Z test statistic and CP using Lachin's formulas
# If use n_int/2 in the SE formula, the CP would be lower than use n1 n2
# SE <- sqrt(p * (1 - p) * (1/(n_int/2) + 1/(n_int/2)))
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)
# If use formulas below to calculate se and cp with n1 and n2, the results are very similar to Lachin's formulas above with n_int/2
#se <- sqrt(p * (1 - p) * (1 / n1 + 1 / n2))
#z_observed <- delta / se
#cp <- pnorm(z_observed - z_alpha) + pnorm(-z_observed - z_alpha)
# Decision rules
if (is.na(cp)) {next}
else {
if (cp > CL & cp < CU) {
cat="Promising"
# Calculate the required sample size
n_star <- ceiling((z_total^2 * (p * (1 - p) * 2)) / (abs(delta))^2)*2
# Apply max N
n_new <- ifelse(n_star<=2*n, n_star, 2*n)
# add n
n_add <- ceiling((n_new-n)/2)*2
add_event_time1 <- rexp(n_add/2, rate = - log(1-p1) / 76)
add_event_time2 <- rexp(n_add/2, rate = - log(1-p2) / 76)
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_final_time <- add_rec_time + follow_up
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, final_time=add_final_time, subj=add_subj)
comb_data <- rbind(all_data, add_data)
final_data <- comb_data %>% subset(select=c(study_time, event_time, rec_time, group, interim_time, final_time, subj))
# event time compare with final time to get indicator: 1 = event, 0 = censored
final_data$cens <- ifelse(final_data$study_time <= final_data$final_time, 1, 0)
}
else {
n_new=n
final_data <- all_data %>% subset(select=c(study_time, event_time, rec_time, group, interim_time, final_time, subj))
# event time compare with final time: 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"}
}
count_simu <- count_simu + 1
# 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
# Standard error
SE_final <- sqrt(final_p * (1 - final_p) * (1/(n_new/2) + 1/(n_new/2)))
# Z statistic
z_final <- (final_p1 - final_p2) / SE_final
# P-value
p_value <- 2 * (1 - pnorm(abs(z_final)))
hypo <- ifelse(p_value < alpha, 1, 0)
final_results <- rbind(final_results, data.frame(p = p, p1 = p1, p2 = p2, cp = cp, cat = cat, n_new = n_new, p_value = p_value, hypo = hypo, count_simu = count_simu, n_int = nrow(surv_data)))
}
}
}
# interim analysis results printouts
est_n <- mean(final_results$n_new)
cp_mean <- mean(final_results$cp)
# Calculate the frequency and proportion of each category
cat_count <- table(final_results$cat)
prop_cat <- prop.table(cat_count)
# Print out some results
cat("New Sample Size / Original Sample Size:", est_n, "/", n, "(", est_n/n, ")", "\n")
cat("Mean Conditional Power:", cp_mean, "\n")
print(prop_cat)
# final analysis results printouts
sig <- sum(final_results$hypo)/max(final_results$count_simu)
nosig <- 1 - sig
cat("Probability of claim success:", sig, "\n")
cat("Proportion of not claim success:", nosig, "\n")
}
### under alternative hypothesis and type II error
simu(alpha=0.05, beta=0.2, rate1=0.7, rate2=0.3, t1=0.5, real_rate2=0.2, rec=1, CL=0.2, CU=0.8, num_simulations=100)
## New Sample Size / Original Sample Size: 56.77083 / 50 ( 1.135417 )
## Mean Conditional Power: 0.8397364
##
## Favorable Promising Unfavorable
## 0.7708333 0.1354167 0.0937500
## Probability of claim success: 0.9583333
## Proportion of not claim success: 0.04166667
### 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_rate2=0.3, rec=1, CL=0.2, CU=0.8, num_simulations=100)
## New Sample Size / Original Sample Size: 55.61224 / 50 ( 1.112245 )
## Mean Conditional Power: 0.6481831
##
## Favorable Promising Unfavorable
## 0.5816327 0.1122449 0.3061224
## Probability of claim success: 0.9081633
## Proportion of not claim success: 0.09183673
### under null hypothesis and type I error
simu(alpha=0.05, beta=0.2, rate1=0.7, rate2=0.3, t1=0.5, real_rate2=0.7, rec=1, CL=0.2, CU=0.8, num_simulations=100)
## New Sample Size / Original Sample Size: 50.60976 / 50 ( 1.012195 )
## Mean Conditional Power: 0.01472299
##
## Promising Unfavorable
## 0.01219512 0.98780488
## Probability of claim success: 0.06097561
## Proportion of not claim success: 0.9390244