Setting up the Function

#' @title Simulation of a Two-Group Parallel-Arm Trial With Interim Analyses
#' @docType Custom function for simulation from the rpact package
#' @author Andrew Althouse with edits by Zad Rafi
#' NOTE: If you want to confirm "type 1 error" under different stopping rules,
#' make death = in the two treatment arms (e.g. no treatment effect)
#' NOTE: I have set this one up to test the power for a treatment that would reduce mortality
#' from 40% in control group (1) to 30% in treatment group (2)
#' NOTE: Trial Design Parameters - Part 1
#' Here we will specify the basics: total N patients to enroll, and death rate for each treatment arm
#' NOTE: Trial Design Parameters - Part 2
#' Here we will define the interim analysis strategy and stopping rules
#' For this trial we will include provisions for efficacy stopping only (no pre-specified futility stopping)
#' We will use the rpact package to compute the stopping/success thresholds at the interim and final analysis
#' NOTE: Required packages: rpact and stringr
#' @param nSims # The number of simulations, the default is 1000
#' @param nPatients # here is where you specify the planned max number of patients you want included in each RCT
#' @param death1 # here is where you specify the event rate for patients receiving 'treatment 1' in these trial
#' @param death2 # here is where you specify the event rate for patients receiving 'treatment 2' in these trials
#' @param nLooks # here is where you put the number of looks that will take place (INCLUDING the final analysis)
#' @param analyses_scheduled # schedule of interim analyses
#' @param sided # Whether the test is 1-sided or 2-sided
#' @param alpha # Specified alpha level, the default is 0.05
#' @param informationRates #
#' @param trials # The total number of trials you wish to load in the table results. 
#' @param typeOfDesign # The type of design.
#' @param seed # Argument to set the seed for the simulations
#' @return list of dataframes and plots
#' @example See below this code block
interim_sim <- function(nPatients = 1000, death1 = 0.4, death2 = 0.3, 
                        nLooks = 4, analyses_scheduled = c(0.25, 0.50, 0.75, 1),
                        sided = 1, alpha = 0.05, 
                        informationRates = analyses_scheduled, 
                        typeOfDesign = "asOF", nSims = 1000, trials = 10,
                        seed = 1031) {

efficacy_thresholds <- numeric(nLooks)
design <- getDesignGroupSequential(
  sided = sided, alpha = alpha,
  informationRates = analyses_scheduled,
  typeOfDesign = typeOfDesign)
design_2 <- getDesignGroupSequential(typeOfDesign = "P") # Pocock
design_3 <- getDesignGroupSequential(typeOfDesign = "asP") # Alpha-spending Pocock
design_4 <- getDesignGroupSequential(typeOfDesign = "OF") # O'Brien-Fleming
designSet <- getDesignSet(designs = c(design, design_2,
                                      design_3, design_4),
                          variedParameters = "typeOfDesign")

RNGkind(kind = "L'Ecuyer-CMRG")
set.seed(seed)
for (j in 1:nLooks) {
  efficacy_thresholds[j] <- design$stageLevels[j]
}

analyses_nPatients <- analyses_scheduled * nPatients
efficacy_thresholds
pb <- txtProgressBar(min = 0, max = nSims, initial = 0, style = 3)
trialnum <- numeric(nSims)
or <- data.frame(matrix(ncol = nLooks, nrow = nSims))
lcl <- data.frame(matrix(ncol = nLooks, nrow = nSims))
ucl <- data.frame(matrix(ncol = nLooks, nrow = nSims))
pval <- data.frame(matrix(ncol = nLooks, nrow = nSims))
success <- data.frame(matrix(ncol = nLooks, nrow = nSims))

strings <- c("OR_%d", "LCL_%d", "UCL_%d",
             "Pval_%d", "Success_%d")

colnames(or) <- sprintf("OR_%d", (1:nLooks))
colnames(lcl) <- sprintf("LCL_%d", (1:nLooks))
colnames(ucl) <- sprintf("UCL_%d", (1:nLooks))
colnames(pval) <- sprintf("Pval_%d", (1:nLooks))
colnames(success) <- sprintf("Success_%d", (1:nLooks))
overall_success <- numeric(nSims)
df <- data.frame(trialnum, or, lcl, ucl,
                 pval, success, overall_success)

RNGkind(kind = "L'Ecuyer-CMRG")
set.seed(seed)
time <- system.time(
for (i in 1:nSims) {
  trialnum[i] <- i
  pid <- seq(1, nPatients, by = 1)
  treatment <- rep(1:2, nPatients / 2)
  deathprob <- numeric(nPatients)
  deathprob[treatment == 1] <- death1
  deathprob[treatment == 2] <- death2
  death <- rbinom(nPatients, 1, deathprob)
  trialdata <- data.frame(cbind(pid, treatment, death))

  for (j in 1:nLooks) {
    analysisdata <- subset(trialdata, pid <= analyses_nPatients[j])
    model <- glm(death ~ treatment,
                 family = binomial(link = "logit"),
                 data = analysisdata)
    or[i, j] <- exp(summary(model)$coefficients[2])
    lcl[i, j] <- exp(confint.default((model))[2, 1])
    ucl[i, j] <- exp(confint.default((model))[2, 2])
    pval[i, j] <- summary(model)$coefficients[8]
    success[i, j] <- ifelse(or[i, j] < 1 & pval[i, j] < efficacy_thresholds[j], 1, 0)
  }
  overall_success[i] <- 0
  for (j in 1:nLooks) {
    if (success[i, j] == 1) {
      overall_success[i] <- 1
    }
  }
    setTxtProgressBar(pb, i)
})

df <- data.frame(trialnum, or, lcl, ucl, pval,
                 success, overall_success)
simulation_results <- data.frame(matrix(vector(),
  nrow = nPatients, ncol = (length(df))))
colnames(simulation_results) <- c("trialnum", (do.call(
  rbind,
  lapply(1:length(strings),
    FUN = function(j) {
      (do.call(rbind, lapply(1:nLooks,
        FUN = function(i) ((sprintf(strings, i)))
      )[]))[, j]
    }
  )
)), "overall_success")
simulation_results[intersect(names(df), names(simulation_results))] <- (df[intersect(
  names(df),
  names(simulation_results)
)])
simulation_results <- as.data.frame(simulation_results)

outputs <- simulation_results[, c(-1, -length(simulation_results))]

cols <- as.character(1:nLooks)
rows <- as.character(1:length(strings))
interim <- matrix(
  nrow = length(strings), ncol = nLooks,
  dimnames = list((rows), (cols)), 
  data = rep(0, (length(strings) * nLooks)))

for (i in cols) {
  interim[, i] <- str_subset(colnames(outputs), i)
}

colnames(interim) <- sprintf("Interim_Look_%d", (1:nLooks))
colnames(simulation_results)[1] <- c("TrialNum")
colnames(simulation_results)[length(simulation_results)] <- c("Overall_Success")
simulation_results_trials <- head(simulation_results, trials)
results <- list(
  summary(design),
  summary(time),
  plot(design, 1) +
    theme_less(),
  plot(designSet, type = 1) +
    theme_less(),
  table(overall_success),
  simulation_results_trials)

names(results) <- c(
  "Design Summary", "Time to complete simulation", 
  "Main Design Plot", "Plot of Various Designs", 
  "Table of Overall Success", "Simulation Results")

return(results)
}

A Simulated Example

results <- interim_sim(nPatients = 1000, death1 = 0.4, death2 = 0.3, nLooks = 4,
                       analyses_scheduled = c(0.25, 0.50, 0.75, 1),
                       sided = 1, alpha = 0.025, typeOfDesign = "asOF",
                       informationRates =  c(0.25, 0.50, 0.75, 1),
                       nSims = 1000, trials = 20, seed = 1031)

Examining the Results

results[1:5]
#> $`Design Summary`
#> 
#> -----
#> 
#> *Sequential analysis with a maximum of 4 looks (group sequential design)*
#> 
#> O'Brien & Fleming type alpha spending design, one-sided overall 
#> significance level 2.5%, power 80%, undefined endpoint, inflation factor 1.0196, 
#> ASN H1 0.8387, ASN H01 0.9798, ASN H0 1.0168.
#> 
#> | Stage                             |       1 |       2 |       3 |       4 |
#> | ----- | ----- | ----- | ----- | ----- |
#> | Planned information rate          |     25% |     50% |     75% |    100% |
#> | Efficacy boundary (z-value scale) |   4.333 |   2.963 |   2.359 |   2.014 |
#> | Stage levels (one-sided)          | <0.0001 |  0.0015 |  0.0092 |  0.0220 |
#> | Cumulative alpha spent            | <0.0001 |  0.0015 |  0.0096 |  0.0250 |
#> | Cumulative power                  |  0.0018 |  0.1679 |  0.5400 |  0.8000 | 
#> 
#> $`Time to complete simulation`
#>    user  system elapsed 
#>   6.824   0.491   7.526 
#> 
#> $`Main Design Plot`

#> 
#> $`Plot of Various Designs`

#> 
#> $`Table of Overall Success`
#> overall_success
#>   0   1 
#> 134 866

To examine the results, I used the kableExtra package, though this is not necessary, and simply using the following script will suffice

table_results <- results[[6]]
View(table_results) 

Statistical Environment

The analyses were run on:

si <- sessionInfo()
print(si, RNG = TRUE, locale = TRUE)
#> R version 4.4.1 (2024-06-14)
#> Platform: aarch64-apple-darwin20
#> Running under: macOS Sonoma 14.6.1
#> 
#> Matrix products: default
#> BLAS:   /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/lib/libRblas.0.dylib 
#> LAPACK: /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/lib/libRlapack.dylib;  LAPACK version 3.12.0
#> 
#> Random number generation:
#>  RNG:     L'Ecuyer-CMRG 
#>  Normal:  Inversion 
#>  Sample:  Rejection 
#>  
#> locale:
#> [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
#> 
#> time zone: America/New_York
#> tzcode source: internal
#> 
#> attached base packages:
#>  [1] splines   grid      stats4    parallel  stats     graphics  grDevices utils     datasets  methods   base     
#> 
#> other attached packages:
#>  [1] cli_3.6.3             texPreview_2.1.0      tinytex_0.53          rmarkdown_2.28        brms_2.21.0          
#>  [6] bootImpute_1.2.1      knitr_1.48.1          boot_1.3-31           gtsummary_2.0.2       reshape2_1.4.4       
#> [11] ProfileLikelihood_1.3 ImputeRobust_1.3-1    gamlss_5.4-22         gamlss.dist_6.1-1     gamlss.data_6.0-6    
#> [16] mvtnorm_1.3-1         performance_0.12.3    summarytools_1.0.1    rpact_4.0.0           tidybayes_3.0.7      
#> [21] htmltools_0.5.8.1     Statamarkdown_0.9.2   car_3.1-2             carData_3.0-5         qqplotr_0.0.6        
#> [26] ggcorrplot_0.1.4.1    mitml_0.4-5           pbmcapply_1.5.1       Amelia_1.8.2          Rcpp_1.0.13          
#> [31] blogdown_1.19.1       doParallel_1.0.17     iterators_1.0.14      foreach_1.5.2         lattice_0.22-6       
#> [36] bayesplot_1.11.1      wesanderson_0.3.7     VIM_6.2.2             colorspace_2.1-1      here_1.0.1           
#> [41] progress_1.2.3        loo_2.8.0             mi_1.1                Matrix_1.7-0          broom_1.0.6          
#> [46] yardstick_1.3.1       svglite_2.1.3         Cairo_1.6-2           cowplot_1.1.3         mgcv_1.9-1           
#> [51] nlme_3.1-166          xfun_0.47             broom.mixed_0.2.9.5   reticulate_1.39.0     kableExtra_1.4.0     
#> [56] posterior_1.6.0       checkmate_2.3.2       parallelly_1.38.0     miceFast_0.8.2        randomForest_4.7-1.1 
#> [61] missForest_1.5        miceadds_3.17-44      quantreg_5.98         SparseM_1.84-2        MCMCpack_1.7-1       
#> [66] MASS_7.3-61           coda_0.19-4.1         latex2exp_0.9.6       rstan_2.32.6          StanHeaders_2.32.10  
#> [71] lubridate_1.9.3       forcats_1.0.0         stringr_1.5.1         dplyr_1.1.4           purrr_1.0.2          
#> [76] readr_2.1.5           tibble_3.2.1          ggplot2_3.5.1         tidyverse_2.0.0       ggtext_0.1.2         
#> [81] concurve_2.8.0        showtext_0.9-7        showtextdb_3.0        sysfonts_0.8.9        future.apply_1.11.2  
#> [86] future_1.34.0         tidyr_1.3.1           magrittr_2.0.3        mice_3.16.0           rms_6.8-2            
#> [91] Hmisc_5.1-3          
#> 
#> loaded via a namespace (and not attached):
#>   [1] nnet_7.3-19             TH.data_1.1-2           vctrs_0.6.5             digest_0.6.37           png_0.1-8              
#>   [6] shape_1.4.6.1           proxy_0.4-27            magick_2.8.4            fontLiberation_0.1.0    gWidgets2tcltk_1.0-8   
#>  [11] withr_3.0.1             ggpubr_0.6.0            survival_3.7-0          doRNG_1.8.6             memoise_2.0.1          
#>  [16] emmeans_1.10.4          MatrixModels_0.5-3      systemfonts_1.1.0       ragg_1.3.3              zoo_1.8-12             
#>  [21] V8_5.0.0                ggdist_3.3.2            DEoptimR_1.1-3          Formula_1.2-5           prettyunits_1.2.0      
#>  [26] rematch2_2.1.2          httr_1.4.7              rstatix_0.7.2           globals_0.16.3          rstudioapi_0.16.0      
#>  [31] extremevalues_2.3.4     pan_1.9                 generics_0.1.3          base64enc_0.1-3         curl_5.2.2             
#>  [36] mitools_2.4             desc_1.4.3              xtable_1.8-4            svUnit_1.0.6            pracma_2.4.4           
#>  [41] evaluate_0.24.0         hms_1.1.3               glmnet_4.1-8            bookdown_0.40           lmtest_0.9-40          
#>  [46] robustbase_0.99-4       matrixStats_1.4.1       svgPanZoom_0.3.4        class_7.3-22            gWidgets2_1.0-9        
#>  [51] pillar_1.9.0            caTools_1.18.3          compiler_4.4.1          stringi_1.8.4           jomo_2.7-6             
#>  [56] minqa_1.2.8             plyr_1.8.9              crayon_1.5.3            abind_1.4-8             metadat_1.2-0          
#>  [61] sp_2.1-4                mathjaxr_1.6-0          rapportools_1.1         twosamples_2.0.1        sandwich_3.1-1         
#>  [66] whisker_0.4.1           codetools_0.2-20        multcomp_1.4-26         textshaping_0.4.0       bcaboot_0.2-3          
#>  [71] openssl_2.2.1           flextable_0.9.6         bslib_0.8.0             QuickJSR_1.3.1          e1071_1.7-16           
#>  [76] gridtext_0.1.5          utf8_1.2.4              lme4_1.1-35.5           fs_1.6.4                itertools_0.1-3        
#>  [81] listenv_0.9.1           pkgbuild_1.4.4          estimability_1.5.1      ggsignif_0.6.4          tzdb_0.4.0             
#>  [86] pkgconfig_2.0.3         tools_4.4.1             cachem_1.1.0            viridisLite_0.4.2       DBI_1.2.3              
#>  [91] numDeriv_2016.8-1.1     fastmap_1.2.0           scales_1.3.0            sass_0.4.9              officer_0.6.6          
#>  [96] opdisDownsampling_1.0.1 insight_0.20.4          rpart_4.1.23            farver_2.1.2            survminer_0.4.9        
#>  [ reached getOption("max.print") -- omitted 58 entries ]

Article Citation