In mjuraska/seqDesign: Simulation and Group-Sequential Monitoring of Randomized Treatment Efficacy Trials with Time-to-Event Endpoints

\clearpage

rm(list=ls(all=TRUE))

knitr::opts_chunk$set(echo=FALSE, message=FALSE, warning=FALSE, cache=TRUE, fig.align='center', fig.width=8, fig.asp=0.5, out.width='100%')

# setting 'scipen=999' suppresses scientific format of large numbers
options(scipen=999)

library(CoRpower)
library(xtable)
library(ggplot2)
library(ldbounds)

# source seqDesign from your local master branch (this has several advantages at this stage of code development over sourcing from the remote master branch)
if (Sys.getenv("USERNAME")=="mjuraska"){
  source("h:/SCHARP/seqDesign/code/source_seqDesign.R")
} else {
  source("C:/Users/brian/Documents/R/source_seqDesign.R")  
}

# THE BELOW USER-SPECIFIED INPUT PARAMETERS ARE AUTOMATICALLY APPLIED TO THE ENTIRE REPORT ------------------

# a directory storing output from simTrial, monitorTrial, and censTrial
srcDir <- paste0(ifelse(Sys.getenv("USERNAME")=="mjuraska", "u:", "v:"), "/bsimpkin/seqDesign/Routput/NIAID/1to1VtoPratio/v5")

N.pla <- 13500          
N.vax <- 13500          
N.vax.arms <- 1      

# annual placebo incidence rate
infRate <- 0.028 / 0.9
infRate.fmt <- format(infRate, digits=3, nsmall=3)

# interim monitoring for benefit/efficacy
eventTarget <- 197
infFractionEff <- c(0.35, 0.7, 1)
atEventCount <- ceiling(infFractionEff * eventTarget)

# interim monitoring for lack of benefit/non-efficacy
eventTargetNonEff <- 130
infFractionNonEff <- c(0.2, 0.6, 1)
atEventCountNonEff <- ceiling(infFractionNonEff * eventTargetNonEff)

estimand <- "cox"

### efficacy monitoring ----------------------
# one-sided Wald test of H0: VE <= effNullVE
effNullVE <- 0.3

# VE design alternative for powering the study
effDesignAltVE <- 0.6

# boundary type: options are "OBF" or "Pocock"
effBound <- "OBF"

# overall one-sided alpha
effAlpha <- 0.025

# monitorTrial's output file name for the VE design alternative scenario used for computing cumulative probability of crossing efficacy boundary
cumProbCrossEff.monitorTrialFile <- paste0("monitorTrial_nPlac=", N.pla, "_nVacc=", N.vax, "_aveVE=", effDesignAltVE, "_infRate=", infRate.fmt, "_", estimand, ".RData")

### non-efficacy monitoring ------------------
# one-sided Wald test of H0: VE >= nonEffNullVE
nonEffNullVE <- 0.5

# boundary type: options are "FKG", "OBF" or "Pocock"
nonEffBound <- "FKG"

# overall one-sided alpha
nonEffAlpha <- 0.025

# monitorTrial's output file name for the VE = 0 scenario used for computing cumulative probability of crossing non-efficacy boundary
cumProbCrossNonEff.monitorTrialFile <- paste0("monitorTrial_nPlac=", N.pla, "_nVacc=", N.vax, "_aveVE=0_infRate=", infRate.fmt, "_", estimand, ".RData")

# monitorTrial's output file name for the VE design alternative scenario used for computing cumulative probability of crossing non-efficacy boundary
cumProbCrossNonEffUnderH1.monitorTrialFile <- paste0("monitorTrial_nPlac=", N.pla, "_nVacc=", N.vax, "_aveVE=", effDesignAltVE, "_infRate=", infRate.fmt, "_", estimand, ".RData")

# END OF USER-SPECIFIED VALUES --------------------------------------------------------------------

### efficacy boundary characterization -----------------------
# *one-sided* nominal significance levels for efficacy
tmp <- bounds(t=infFractionEff, iuse=ifelse(effBound=="OBF", 1, 2), alpha=effAlpha)
nominalAlphaEff <- 1 - pnorm(tmp$upper.bounds)

# previously, RCTdesign was used
# load(file.path(srcDir, paste0("monitoringBoundariesAt", paste(atEventCount, collapse="_"), "events.RData")))
# HRboundEff <- bounds$eff
# HRboundNonEff <- bounds$noneff

HRboundEff <- estHRbound(boundType="eff", nullHR=1 - effNullVE, alpha=2 * nominalAlphaEff, nEvents=atEventCount, randFrac=N.vax / (N.vax + N.pla))

cumProbCrossEff <- crossBoundCumProb(boundType="eff", nAnalyses=length(infFractionEff),
                                     monitorTrialFile=cumProbCrossEff.monitorTrialFile,
                                     monitorTrialDir=srcDir)


### non-efficacy boundary characterization -------------------
# *one-sided* nominal significance levels for non-efficacy
if (nonEffBound=="FKG"){
  nominalAlphaNonEff <- rep(nonEffAlpha, 3)  
} else {
  tmp <- bounds(t=infFractionNonEff, iuse=ifelse(nonEffBound=="OBF", 1, 2), alpha=nonEffAlpha)
  nominalAlphaNonEff <- 1 - pnorm(tmp$upper.bounds)  
}

HRboundNonEff <- estHRbound(boundType="nonEff", nullHR=1 - nonEffNullVE, alpha=2 * nominalAlphaNonEff, nEvents=atEventCountNonEff, randFrac=N.vax / (N.vax + N.pla))

cumProbCrossNonEff <- crossBoundCumProb(boundType="nonEff", nAnalyses=length(infFractionNonEff),
                                        monitorTrialFile=cumProbCrossNonEff.monitorTrialFile,
                                        monitorTrialDir=srcDir)
cumProbCrossNonEffUnderH1 <- crossBoundCumProb(boundType="nonEff", nAnalyses=length(infFractionNonEff),
                                               monitorTrialFile=cumProbCrossNonEffUnderH1.monitorTrialFile,
                                               monitorTrialDir=srcDir)


# formatting numeric output
N <- N.pla + N.vax
infFractionEff <- infFractionEff * 100
nominalAlphaEff <- format(round(nominalAlphaEff, digits=6), nsmall=6)
effVE <- format((1 - HRboundEff) * 100, digits=1, nsmall=1)
# HRboundEff <- format(HRboundEff, digits=4, nsmall=4)
cumProbCrossEff <- format(cumProbCrossEff * 100, digits=1, nsmall=1)

infFractionNonEff <- infFractionNonEff * 100
nominalAlphaNonEff <- format(round(nominalAlphaNonEff, digits=3), nsmall=3)
nonEffVE <- format((1 - HRboundNonEff) * 100, digits=1, nsmall=1)

# HRboundNonEff <- ifelse(is.na(HRboundNonEff), "n/a", format(HRboundNonEff, digits=4, nsmall=4))
cumProbCrossNonEff <- format(cumProbCrossNonEff * 100, digits=1, nsmall=1)
cumProbCrossNonEffUnderH1 <- ifelse(cumProbCrossNonEffUnderH1<0.001, "$<0.1$", format(cumProbCrossNonEffUnderH1 * 100, digits=1, nsmall=1))

Simulation Set-up

This report envisages a design with a \emph{single} vaccine arm; all operating characteristics extend to each of multiple vaccine arms assuming independent evaluation of efficacy in each vaccine arm vs. the shared placebo arm
$N = r N$ SARS-CoV-2-negative participants at baseline, randomly allocated in the 1:1 ratio to vaccine vs. placebo (if 10\% of the randomized participants are SARS-CoV-2-positive at baseline, then the total sample size is $r N / 0.9 = r N/0.9$)
All operating characteristics are for the r N baseline negative participants (r N.vax in the vaccine arm, r N.pla in the placebo arm)
r format(infRate * 100, digits=1, nsmall=1)\% annual COV-DIS incidence rate in the placebo arm
- Chosen as the minimal incidence rate satisfying the requirement that the target numbers of FAS-15d+ endpoints under both 2:1 and 1:1 allocation ratios be accrued with 80\% probability by 6.5 months since trial initiation in the scenario with $\overline{VE}$(2--32 weeks)$\,=60\%$ (see below for description of the VE model)
r atEventCount[length(atEventCount)] primary endpoints occurring $\geq 15$ days after the first dose in the FAS cohort (henceforth referred to as the FAS-15d+ endpoints) is the target endpoint total for the event-driven primary analysis
Efficacy monitoring at r infFractionEff[1]\%, r infFractionEff[2]\%, and r infFractionEff[3]\% of the target r atEventCount[length(atEventCount)] FAS-15d+ endpoints in the FAS cohort
Non-efficacy monitoring at r infFractionNonEff[1]\%, r infFractionNonEff[2]\%, and r infFractionNonEff[3]\% of approximately r atEventCountNonEff[length(atEventCountNonEff)] endpoints occurring $\geq 15$ days after the second dose (henceforthh referred to as the PP-15d+ endpoints) in the PP cohort
- The PP cohort includes all FAS participants who received the second dose
3-month uniform accrual
25 months of participant follow-up (i.e., 2 years of follow-up post-vaccination)
Stepwise VE model (all tables and figures are labeled by the time-averaged VE level over the time interval of (2, 32] weeks of follow-up denoted by $\overline{VE}$(2--32)):

\begin{center} \begin{tabular}{cc} \hline \hline Time Period & \ (Weeks) & VE Level\ \hline (0, 2]\T & $0.1 \times (45/43) \times \overline{VE}$(2--32)\ (2, 6] & $(2/3) \times (45/43) \times \overline{VE}$(2--32)\ (6, 32] & $(45/43) \times \overline{VE}$(2--32)\ (32, 108] & $\overline{VE}$(2--32)\ \hline \hline \end{tabular} \end{center}

2\% annual dropout in each arm
2\% probability of a missed second dose (used in definition of the PP cohort)
Group-sequential monitoring for benefit/efficacy: a group-sequential one-sided Wald test of $H_0$: $\mathrm{VE} \leq r effNullVE * 100\%$ using the Cox model and the O'Brien-Fleming boundary
- Analyses performed when r atEventCount[1] (r infFractionEff[1]\% information fraction), r atEventCount[2] (r infFractionEff[2]\% information fraction), and r atEventCount[length(atEventCount)] FAS-15d+ endpoints have been accrued (r atEventCount[length(atEventCount)] is the target number for powering the trial for $\overline{VE}\textrm{(2--32)} = 60\%$)
- At the time of each analysis, testing of $H_0$ is performed in the FAS cohort including only FAS-15d+ endpoints
- Reject $H_0$ if the lower bound of the monitoring-adjusted 95\% Wald CI for VE is $> r effNullVE * 100\%$
- 0.025 overall type 1 error rate for 3 one-sided tests (2 interim and a primary analysis)
Interim monitoring for lack of benefit/non-efficacy: a one-sided Wald test of $H_0$: $\mathrm{VE} \geq r nonEffNullVE * 100\%$ using the Cox model and nominal 95\% confidence intervals following Freidlin, Korn, and Gray (2010, Clin Trials) (referred to as the FKG approach)
- Analyses performed when approximately r atEventCountNonEff[1] (r infFractionNonEff[1]\% information fraction), r atEventCountNonEff[2] (r infFractionNonEff[2]\% information fraction), and r atEventCountNonEff[length(atEventCountNonEff)] PP-15d+ endpoints have been accrued
- Timing of analyses harmonized with sequential monitoring for benefit/efficacy
- At the time of each analysis, testing of $H_0$ is performed in the PP cohort including only PP-15d+ endpoints
- Reject $H_0$ if the upper bound of the nominal 95\% Wald CI for VE is $< r nonEffNullVE * 100\%$
- Nominal 0.025 type 1 error rate of each of the 3 one-sided tests

\clearpage

Summary of efficacy monitoring boundaries (using the Wald test)

tab <- data.frame(atEventCount, infFractionEff=paste0(infFractionEff, "\\%"), nominalAlphaEff, 
                  estVE=paste0("$\\mathrm{VE} \\geq ", effVE[, 1], "\\%$"),
                  ciVE=paste0("$(", effVE[, 3], "\\%, ", effVE[, 2], "\\%)$"),
                  cumProbCross=paste0(cumProbCrossEff, "\\%"))

print(xtable(tab,
             align=rep("c", NCOL(tab) + 1),
             digits=c(0, 0, 1, 6, 0, 0, 1),
             caption=NULL),
      sanitize.text.function=function(x){ x },
      math.style.negative=TRUE,
      include.rownames=FALSE,
      include.colnames=FALSE,
      hline.after=NULL,
      # 'add.to.row must be a list with 2 components 'pos' and 'command'
      add.to.row=list(pos=list(0, NROW(tab)),
                      command=c(paste0("\\hline \\hline  
                                        & & \\multicolumn{4}{c}{Efficacy Monitoring (O'Brien-Fleming)} \\\\\n
                                        &  & \\multicolumn{4}{c}{(Rejecting $H_0$: $\\mathrm{VE} \\leq ", effNullVE * 100, "\\%$)} \\\\\n
                                       \\cline{3-6}
                                       FAS-15d+ & Inf & Nominal & Est. VE & Adj. ", (1 - 2 * effAlpha) * 100, "\\% CI* & Cum Prob of Crossing \\\\\n
                                       Endpoints & Frac & Signif Level & at Bndary & at Bndary & Bndary if $\\overline{VE}\\textrm{(2--32)}=", effDesignAltVE * 100, "\\%$ \\\\\n
                                       \\hline "),
                                paste0("\\hline \\hline
                                       \\multicolumn{6}{l}{\\footnotesize * Group-sequential monitoring-adjusted ", (1 - 2 * effAlpha) * 100, "\\% CI}"))))

Summary of non-efficacy monitoring boundaries (using the Wald test)

tab <- data.frame(atEventCountNonEff, infFractionNonEff=paste0(infFractionNonEff, "\\%"), nominalAlphaNonEff, 
                  estVE=paste0("$\\mathrm{VE} \\leq ", nonEffVE[, 1], "\\%$"),
                  ciVE=paste0("$(", nonEffVE[, 3], "\\%, ", nonEffVE[, 2], "\\%)$"),
                  cumProbCross=paste0(cumProbCrossNonEff, "\\%"))

print(xtable(tab,
             align=rep("c", NCOL(tab) + 1),
             digits=c(0, 0, 1, 6, 0, 0, 1),
             caption=NULL),
      sanitize.text.function=function(x){ x },
      math.style.negative=TRUE,
      include.rownames=FALSE,
      include.colnames=FALSE,
      hline.after=NULL,
      # 'add.to.row must be a list with 2 components 'pos' and 'command'
      add.to.row=list(pos=list(0, NROW(tab)),
                      command=c(paste0("\\hline \\hline  
                                        & & \\multicolumn{4}{c}{Non-Efficacy Monitoring (FKG)} \\\\\n
                                        &  & \\multicolumn{4}{c}{(Rejecting $H_0$: $\\mathrm{VE} \\geq ", nonEffNullVE * 100, "\\%$)} \\\\\n
                                       \\cline{3-6}
                                       PP-15d+ & Inf & Nominal & Est. VE & ", ifelse(nonEffBound=="FKG", "", "Adj. "), (1 - 2 * effAlpha) * 100,"\\% CI* & Cum Prob of Crossing \\\\\n
                                       Endpoints & Frac & Signif Level & at Bndary & at Bndary & Bndary if $\\overline{VE}\\textrm{(2--32)}=0\\%$ \\\\\n
                                       \\hline "),
                                paste0("\\hline \\hline
                                       \\multicolumn{6}{l}{\\footnotesize * ", ifelse(nonEffBound=="FKG", "Nominal ", "Group-sequential monitoring-adjusted "), (1 - 2 * nonEffAlpha) * 100, "\\% CI}"))))

Cumulative probabilities of crossing the non-efficacy boundary at r paste0(atEventCountNonEff[1], ", ", atEventCountNonEff[2], ", and ", atEventCountNonEff[3]) PP-15d+ endpoints if $\overline{VE}\textrm{(2--32)}=r effDesignAltVE * 100\%$ are r cumProbCrossNonEffUnderH1[1]\%, r cumProbCrossNonEffUnderH1[2]\%, and r cumProbCrossNonEffUnderH1[3]\%
Interim monitoring for potential vaccine harm/disease enhancement: monitoring-adjusted one-sided conditional exact binomial test of $H_0$: $p_v \geq N_v/N$, where $p_v$ is the binomial probability of assignment to the vaccine arm conditional on the observed number of endpoints, and $N_v/N$ is the proportion of participants randomized to the vaccine arm
- Starting at the 12$^{\textrm{th}}$ FAS-15d+ endpoint and then continuously (i.e., after each FAS-15d+ endpoint) until the first interim analysis for lack of benefit (at r atEventCount[1] FAS-15d+ endpoints)
- 0.05 overall type 1 error rate and a constant nominal significance level of each test
1000 Monte-Carlo iterations
Software: R package seqDesign available at https://github.com/mjuraska/seqDesign

\clearpage

Power Calculations

# 'aveVE.vector' determines what results are reported in each row
aveVE.vector <- c(-2, -1.5, -1, -0.5, 0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9)

# generate a shell for the results
res <- as.data.frame(matrix(0, ncol=6, nrow=length(aveVE.vector)))
colnames(res) <- c("aveVE", "aveHR", "harm", "noneffInterim", "noneffFinal", "higheff")
res[,1] <- paste0(aveVE.vector*100,"%")
res[,2] <- 1 - aveVE.vector
if (!is.null(effNullVE)){ designPower <- as.data.frame(matrix(0, nrow=length(aveVE.vector), ncol=length(effNullVE))) }
colnames(designPower) <- paste0("design",1:length(effNullVE))
powerPP <- numeric(length(aveVE.vector))

# extract the estimated probabilities of each trial outcome and, if desired, unconditional power to reject H0 defined by 'effNullVE' from the output of 'monitorTrial'
row <- 0
for (aveVE in aveVE.vector){
  row <- row + 1
  filename <- paste0("monitorTrial_nPlac=", N.pla, "_nVacc=", N.vax, "_aveVE=", aveVE, "_infRate=",format(infRate, digits=3, nsmall=3),"_",estimand,".RData")
  load(file.path(srcDir, filename))

  bounds <- sapply(out, function(trial) trial[[1]]$boundHit)
  bounds.freq <- table(bounds, useNA="ifany")/length(out)
  res[row,3:6] <- bounds.freq[c("Harm","NonEffInterim","NonEffFinal","HighEff")] 
  res[row,] <- ifelse(is.na(res[row,]), 0, res[row,])

  if (!is.null(effNullVE)){ 
    altDetectedMatrix <- do.call("rbind",lapply(out, function(oneTrial){ oneTrial[[1]]$altDetected }))
    if (is.null(altDetectedMatrix)){
      designPower[row,] <- numeric(length(effNullVE))
    } else {
      designPower[row,] <- colSums(altDetectedMatrix, na.rm=TRUE)/length(out) 
    }    
  }

  filename <- paste0("monitorTrial_effCohort=pp_nPlac=", N.pla, "_nVacc=", N.vax, "_aveVE=", aveVE, "_infRate=",format(infRate, digits=3, nsmall=3),"_",estimand,".RData")
  load(file.path(srcDir, filename))

  # 'altDetected' is a logical vector of whether H0 was rejected in the PP cohort (TRUE) or not (FALSE)
  altDetected <- sapply(out, function(oneTrial){ oneTrial[[1]]$altDetected })
  altDetected[is.na(altDetected)] <- FALSE
  powerPP[row] <- mean(altDetected)
}

res$noneff <- res$noneffInterim + res$noneffFinal
res <- subset(res, select=c("aveVE","aveHR","harm","noneffInterim","noneffFinal","noneff"))
if (!is.null(effNullVE)){ res <- data.frame(res, designPower) }

# add unconditional power to reject H0: VE<=20% in the PP cohort (counting PP events starting at 28 + 15 days after enrollment)
#res$powerPP <- c(0, 0, 0, 0, 0.002, 0.007, 0.031, 0.070, 0.232, 0.532, 0.889, 0.984, 1)
res$powerPP <- powerPP

res[,3:NCOL(res)] <- res[,3:NCOL(res)]*100
res1 <- res[, -(5:6)]

# prints a table with probabilities of reaching each trial monitoring outcome and unconditional power to reject H0
print(xtable(res1, 
             align=rep(c("c","c","c"), c(5, ifelse(!is.null(effNullVE), length(effNullVE), 0), 1)), 
             digits=c(1,1,1,rep(1, NCOL(res1) - 2)), 
             caption=paste0("Probabilities ($\\times 100$) of reaching each possible trial monitoring outcome",ifelse(!is.null(effNullVE), paste0(" and unconditional power ($\\times 100$) to reject $H_0$: $VE \\leq ", effNullVE * 100, "\\%$ in the FAS and PP cohorts counting only FAS-15d+ and PP-15d+ endpoints, respectively,"),"")," for a 2-arm study design with ",N.pla," placebo and ",N.vax," vaccine recipients"),
             comment=FALSE
             ), 
      #scalebox=ifelse(!is.null(effNullVE),0.85,1),
      math.style.negative=TRUE,
      caption.placement="top", 
      include.rownames=FALSE, 
      include.colnames=FALSE,
      comment=FALSE,
      hline.after=4, 
      add.to.row=list(pos=list(0, NROW(res1)), 
                      command=c(paste0("\\hline \\hline  &  & Harm &  Non-eff ", ifelse(!is.null(effNullVE), paste0("& ", paste0("\\multicolumn{", length(effNullVE) + 1, "}{c}{Uncond Power} ")), ""), "\\\\\n \\cline{5-6} 
                                       $\\overline{VE}$(2--32) & $\\overline{HR}$(2--32) & FAS$^{\\ddag}$ VE$<0\\%$ & PP$^{\\dagger}$ VE$<$", nonEffNullVE * 100, "\\% & FAS$^{\\ddag}$ VE$>$", effNullVE * 100, "\\% & PP$^{\\dagger}$ VE$>$", effNullVE * 100, "\\% \\\\\n \\hline"), 
                                paste0(
                                  "\\hline \\hline 
                                  \\multicolumn{", 5+length(effNullVE), "}{l}{\\footnotesize $^{\\dagger}$ PP cohort counting only PP-15d+ endpoints} \\\\\n
                                  \\multicolumn{", 5+length(effNullVE), "}{l}{\\footnotesize $^{\\ddag}$ FAS cohort counting only FAS-15d+ endpoints} \\\\\n
                                  \\multicolumn{", 5+length(effNullVE), "}{l}{\\footnotesize N=", N.pla, ":", N.vax, " placebo:vaccine group} \\\\\n 
                                  \\multicolumn{", 5+length(effNullVE),"}{l}{\\footnotesize ", format(infRate*100, digits=2, nsmall=2), "\\% annual incidence rate in placebo group} \\\\\n 
                                  \\multicolumn{", 5+length(effNullVE),"}{l}{\\footnotesize 2\\% annual dropout rate in each group} \\\\\n
                                  \\multicolumn{", 5+length(effNullVE),"}{l}{\\footnotesize Cox model-based efficacy monitoring with O'Brien-Fleming boundaries (FAS-15d+ endpoints)} \\\\\n
                                  \\multicolumn{", 5+length(effNullVE), "}{l}{\\footnotesize ", ifelse(estimand=="combined", "Cox \\& cumulative incidence-based non-efficacy monitoring}", ifelse(estimand=="cuminc", "Cumulative incidence-based non-efficacy monitoring incl. post-1.5 months PP monitoring", "Cox model-based non-efficacy monitoring with nominal CIs (PP-15d+ endpoints)}")), " \\\\\n 
                                  \\multicolumn{", 5+length(effNullVE), "}{l}{\\footnotesize ", ifelse(estimand=="cox", "Cox model-based unconditional power", "Cumulative incidence-based FAS unconditional power"),"}"))))

\clearpage

res <- subset(res, aveHR<=1)

par(mar=c(4.5,5,3,1.5), las=1, cex.axis=1.2, cex.lab=1.3, cex.main=1.2)
plot(-1, 0, xlim=c(0, max(aveVE.vector)), ylim=0:1, xaxt="n", yaxt="n", xlab=bquote(paste(bar(VE), "(2-32)")), ylab="Probability")
axis(side=1, at=seq(0, max(aveVE.vector), by=0.1), labels=paste0(seq(0, max(aveVE.vector) * 100, by=10),"%"))
axis(side=2, at=c(0,0.05,0.9,seq(0.2,1,by=0.2)))
axis(side=4, at=c(0,0.05,0.9,seq(0.2,1,by=0.2)), labels=FALSE)
abline(h=c(0.05, 0.9), lwd=2, col="gray60")

with(res, lines(1-aveHR, harm/100, type="b", lwd=3.5, lty="solid", col="darkorange"))
#with(res, lines(1-aveHR, higheff/100, type="b", lwd=3.5, lty="longdash", col="red3"))
with(res, lines(1-aveHR, noneffInterim/100, type="b", lwd=3.5, lty="dotdash", col="red3"))
if (!is.null(effNullVE)){
  for (i in 1:length(effNullVE)){ lines(1-res$aveHR, res[,paste0("design",i)]/100, type="b", lwd=3.5, lty="solid", col="black") }
}  
with(res, lines(1-aveHR, powerPP/100, type="b", lwd=3.5, lty="longdash", col="magenta3"))

position <- ifelse(res$noneffInterim[res$aveHR == 0.9] > res$noneffInterim[res$aveHR == 0.8] & res$noneffInterim[res$aveHR == 1.0] > res$noneffInterim[res$aveHR == 0.8], c(1.2, 1.2), c(-0.1, 0.4))

with(res, text(1-aveHR[aveHR==0.8], noneffInterim[aveHR==0.8] / 100, "Non-efficacy\n at interim analysis", cex=1.2, adj=position))

with(res, text(1-aveHR[aveHR==1], max(harm[aveHR==1] / 100, 0.055), "Potential harm", adj=c(-0.1, -0.1), cex=1.2))

#with(res, text(1-aveHR[aveHR==0.4], higheff[aveHR==0.4] / 100, "High efficacy", adj=c(-0.1, 0.5), cex=1.2))
#with(res, text(1-aveHR[aveHR==0.5], powerPP[aveHR==0.5] / 100, "FAS-15d+\nunconditional power\n[VE>20%]", cex=1.2, adj=c(0.95, -0.4)))
#with(res, text(1-aveHR[aveHR==0.4], design1[aveHR==0.4] / 100, "PP\nunconditional power\n[VE>20%]", cex=1.2, adj=c(0, 1.3)))
#text(0.45, 0.6, "ITT unconditional power\n[VE(0-6)>25%]", cex=1.1, pos=4)

position <- min(res$powerPP[res$aveHR<0.5], res[res$aveHR<0.5, grep('design', names(res))], 90)/100

legend(0.55, position - 0.09, lty=c("solid", "longdash"), lwd=3.5, col=c("black", "magenta3"), title=paste0("Unconditional Power\n[VE < ", effNullVE * 100, "%]"), legend=c("FAS (FAS-15d+ endpoints)", "PP (PP-15d+ endpoints)"), cex=1.2, bty="n")

text(0.6, max(res$harm[res$aveHR<0.5]/100, res$noneffInterim[res$aveHR<0.5]/100, 0.05) + 0.15,
     paste0(
       "N=",N.pla,":",N.vax," pla:vax\n3-month uniform accrual\n",format(infRate*100, digits=2, nsmall=2),"% annual pla incidence\n2% annual dropout\nCox model O'Brien-Fleming eff\n    monitoring (FAS-15d+ endpoints)\n",ifelse(estimand=="combined","Cox & cum inc-based ",ifelse(estimand=="cuminc","Cum inc ","Cox model ")), "FKG non-eff\n    monitoring (PP-15d+ endpoints)\n", ifelse(estimand=="cox","Cox model uncond power","Cum inc uncond power")), cex=1, pos=4)

\clearpage

Interim Monitoring

df <- data.frame(aveVE=paste0(aveVE.vector * 100, "%"), pEffA1=0, pEffA2=0, pEffA3=0)

for (i in 1:length(aveVE.vector)){
  load(file.path(srcDir, paste0("monitorTrial_nPlac=", N.pla, "_nVacc=", N.vax, "_aveVE=", aveVE.vector[i], "_infRate=", format(infRate, digits=3, nsmall=3), "_cox.RData")))

  effCross <- sapply(out, function(x){ ifelse(x[[1]]$boundHit=="Eff", NROW(x[[1]]$summEff), NA) })
  effCross <- factor(effCross, levels=1:3)
  df[i, -1] <- 100 * as.vector(table(effCross)) / length(out)
}

print(xtable(df,
             align=rep("r", 5),
             digits=c(0, 0, 1, 1, 1),
             caption=paste0("Probabilities ($\\times 100$) of rejecting $H_0$: $\\mathrm{VE} \\leq ", effNullVE * 100, "\\%$ at the first interim analysis, second interim analysis, and the primary analysis")
             ),
      math.style.negative=TRUE,
      include.rownames=FALSE,
      include.colnames=FALSE,
      caption.placement="top",
      hline.after=NULL,
      add.to.row=list(pos=list(0, NROW(df)), 
                      command=c(paste0("\\hline \\hline  & \\multicolumn{3}{c}{Rejecting $\\mathrm{VE} \\leq ", effNullVE * 100, "\\%$} \\\\\n 
                                \\cline{2-4} $\\overline{VE}$(2--32) & 1$^{\\textrm{st}}$ IA & 2$^{\\textrm{nd}}$ IA & Prim A \\\\\n \\hline"),
                                "\\hline \\hline")),
      comment=FALSE
      )

#fullVE <- function(aveVE, vePeriods){ aveVE * (vePeriods[3]-1)/((vePeriods[3]-1)-(vePeriods[2]-1)/2) }

pwrPP.filename <- paste0("VEpwPP_nPlac=",N.pla,"_nVacc=",N.vax,"_infRate=",format(infRate, digits=3, nsmall=3),".RData")
load(file.path(srcDir, pwrPP.filename))

par(mar=c(4.5,5,3,1.5), las=1, cex.axis=1.2, cex.lab=1.3, cex.main=1.2)
plot(-1, 0, xlim=c(0,0.8), ylim=0:1, xaxt="n", yaxt="n", xlab="VE(1.5-6 months)", ylab="")
mtext("Per-Protocol Unconditional Power [VE(1.5-6) > 20%]", side=2, line=3.5, cex=1.3, las=0)
axis(side=1, at=seq(0,0.8,by=0.1), labels=paste0(seq(0,80,by=10),"%"))
axis(side=2, at=c(0.9, seq(0,1,by=0.2)))
axis(side=4, at=c(seq(0,1,by=0.2),0.9), labels=FALSE)
abline(h=0.9, lty="dotted", lwd=2)

colNames <- c("black", "darkgoldenrod2", "magenta3")
ltyNames <- c("solid", "dotdash", "longdash")
for (i in 1:3){
  # the x-coordinates are the true values of VE(1.5-6)
  x <- seq(0,0.8,by=0.1)
  # the y-coordinates are powers to reject H0: VE(1.5-6)<=20% in favor of H1: VE(1.5-6)>20%
  y <- sapply(pwList, function(oneAveVEData){ oneAveVEData$VEpwPP[i] })
  lines(x, y, type="b", lty=ltyNames[i], lwd=3.5, col=colNames[i])
}
text(0, 0.7,
     paste0("N=",N.pla,":",N.vax," placebo:vaccine group\n",infRate*100,"% annual placebo group incidence rate\n5% annual dropout rate\nVE halved in first 1.5 months\n", "3-month uniform accrual"),
     cex=1, pos=4)
legend(0.55, 0.4, lty=ltyNames, col=colNames, lwd=3, title="Probability of >=1\nmissed vaccination",
       legend=paste0(c(0,5,10),"%"), bty="n", cex=1.2)

\clearpage

# a list named 'decisionTimes'
monitorTrialFile <- paste0("monitorTrial_nPlac=", N.pla, "_nVacc=", N.vax, "_aveVE=", aveVE.vector, "_infRate=", infRate.fmt, "_cox.RData")
times <- decisionTimes(monitorTrialFile, monitorTrialDir=srcDir)

dframe <- data.frame(aveVE=factor(rep(aveVE.vector, each=1000)), times=do.call(c, times))

ggplot(dframe, aes(x=aveVE, y=times)) +
  geom_boxplot(outlier.shape=NA, size=0.8) +
  stat_summary(fun.y="mean", geom="point", shape=23, size=2, colour="red3", stroke=1.1) +
  xlab(bquote(paste(bar(VE), "(2-32)"))) +
  ylab("Weeks since Trial Start to Crossing\nHarm/Nonefficacy/Efficacy Boundary or\nReaching the Target Number of Endpoints\nIf No Boundary Crossed") +
  scale_x_discrete(breaks=aveVE.vector, labels=aveVE.vector * 100) +
  scale_y_continuous(breaks=seq(4, 32, by=2), expand=c(0.02, 0.02), sec.axis=dup_axis(name=NULL)) +
  coord_cartesian(ylim=c(3, 32)) +
  theme_bw()

\clearpage

# a list named 'eventSplits with 1 component for each VE level
load(file.path(srcDir, "eventSplits.RData"))

medEventCounts <- sapply(eventSplits, function(df){ quantile(rowSums(df), probs=0.5, names=FALSE) })
medEventCountsV <- sapply(eventSplits, function(df){ quantile(df$V, probs=0.5, names=FALSE) })
medEventCountsP <- sapply(eventSplits, function(df){ quantile(df$P, probs=0.5, names=FALSE) })

LB95 <- 1 - exp(log(1 - aveVE.vector) + qnorm(0.975) * sqrt(4 / medEventCounts))
UB95 <- 1 - exp(log(1 - aveVE.vector) - qnorm(0.975) * sqrt(4 / medEventCounts))

tableEventSplitAndCI <- data.frame(aveVE=round(aveVE.vector * 100, digits=0), splitMITT=paste0(round(medEventCountsV, digits=0), ":", round(medEventCountsP, digits=0)), 
                                         LB95=round(LB95 * 100, digits=0), UB95=round(UB95 * 100, digits=0))
tableEventSplitAndCI <- tableEventSplitAndCI[-(1:3), ]

print(xtable(tableEventSplitAndCI,
             align=rep("r", 5),
             digits=rep(0, 5),
             caption="Median endpoints in the FAS cohort occurring $\\geq 28 + 15$ days after the first vaccination (matching the primary analysis cohort) and nominal 95\\% confidence bounds for VE at the time of reporting of the primary analysis results"
             ),
      math.style.negative=TRUE,
      include.rownames=FALSE,
      include.colnames=FALSE,
      caption.placement="top",
      hline.after=NULL,
      add.to.row=list(pos=list(0, NROW(tableEventSplitAndCI)), 
                      command=c("\\hline \\hline Ave & & \\multicolumn{2}{c}{VE} \\\\\n \\cline{3-4} VE & Vaccine: & 95\\% LB & 95\\% UB \\\\\n (\\%) & Placebo* & (\\%) & (\\%) \\\\\n \\hline", 
                                "\\hline \\hline \\multicolumn{4}{l}{\\footnotesize $^{\\ast}$ Median endpoint counts}")),
      comment=FALSE
      )

\clearpage

Trial Duration

stage2 <- 108
cumProbTrialDuration1VaxArm <- function(simTrial.filename, monitorTrial.filename, dirname, stage2=stage2){
  load(file.path(dirname, simTrial.filename))
  trialData <- trialObj[["trialData"]]

  load(file.path(dirname, monitorTrial.filename))

  bounds <- sapply(out, function(trial){ trial[[1]]$boundHit })
  stopTimes <- sapply(out, function(trial){ trial[[1]]$stopTime })
  for (i in 1:length(bounds)){
    if (bounds[i] %in% c("Eff","HighEff")){
      data.i <- trialData[[i]]
      endStage2 <- max(data.i$entry + stage2)
      trialStop <- min(max(data.i$exit), endStage2)
      stopTimes[i] <- trialStop
    }
  }
  return(quantile(stopTimes, prob=seq(0,1,by=0.01)))
}

aveVE <- c(-0.5, seq(0, 0.8, by=0.2))
# store the quantiles of the trial duration for every 'aveVE'
quantiles <- matrix(0, nrow=length(aveVE), ncol=length(seq(0,1,by=0.01)))

for (i in 1:length(aveVE)){
  simTrial.filename <- paste0("simTrial_nPlac=",N.pla,"_nVacc=",N.vax,"_aveVE=",aveVE[i],"_infRate=",format(infRate, digits=3, nsmall=3),".RData")
  monitorTrial.filename <- paste0("monitorTrial_nPlac=",N.pla,"_nVacc=",N.vax,"_aveVE=",aveVE[i],"_infRate=",format(infRate, digits=3, nsmall=3),"_",estimand,".RData")
  quantiles[i,] <- cumProbTrialDuration1VaxArm(simTrial.filename, monitorTrial.filename, srcDir, stage2=stage2)
}

# convert into months
quantiles <- quantiles/(52/12)
quantiles <- cbind(0, quantiles, ifelse(quantiles[, NCOL(quantiles)]<28, 28, NA))

par(mar=c(4.5,5,3,1.5), las=1, cex.axis=1.2, cex.lab=1.3, cex.main=1.2)
plot(-10, 0, xlim=c(0, 28), ylim=0:1, xaxt="n", yaxt="n", 
     xlab="t = Trial Duration (Months)", ylab="Probability (Trial Duration <= t)")
axis(side=1, at=seq(0, 28, by=2))
axis(side=1, at=3)
axis(side=2, at=seq(0,1,by=0.2))
axis(side=4, at=seq(0,1,by=0.2), labels=FALSE)
segments(x0=3, y0=0, y1=0.5, lwd=2, col="gray30")
segments(x0=4, y0=0, y1=1, lwd=2, col="gray30")
#segments(x0=14, y0=0, y1=1, lwd=2, col="gray30")

linesCol <- c("black", "blue", "red", "green", "darkorange", "magenta3")
for (i in 1:length(aveVE)){
  lines(quantiles[i,], c(0, seq(0,1,by=0.01), 1), lwd=3.5, col=linesCol[i])
}

legend(x=12, y=0.55, legend=paste(aveVE*100,"%",sep=""), col=linesCol, lwd=3, 
       cex=1.2, title="Average\nVE(2-32)", y.intersp=0.7, bty="n")
text(18, 0.4, 
     paste0("N=", N.pla, ":", N.vax, " pla:vax\n", format(infRate * 100, digits=2, nsmall=2), "% pla incidence\n 2% dropout rate\n 3 mo accrual"), 
     adj=0, cex=1.2)
text(x=3, y=0.4, labels="Accrual\ncompleted", cex=1.2, pos=2)
text(x=4, y=0.95, labels="Vaccinations\ncompleted", cex=1.2, pos=2)  
#text(x=14, y=0.95, labels="Vaccinations through\nmonth 2 completed", cex=1.2, pos=2)

cumProbTrialDuration <- function(trialDataCens.filename, dirname, stage2=156){
  load(file.path(dirname, trialDataCens.filename))

  trialStopTime <- numeric(length(trialListCensor))
  for (i in 1:length(trialListCensor)){
    dataI <- trialListCensor[[i]]
    trialStopTime[i] <- max(pmin(dataI$entry + stage2, dataI$exit))
  }
  return(quantile(trialStopTime, prob=seq(0,1,by=0.01)))
}

if (N.vax.arms>1){
  aveVE <- matrix(c(0,0,0.4,0,0.4,0.4), nrow=3, ncol=2, byrow=TRUE)

  # store the quantiles of the trial duration for every 'aveVE'
  quantiles <- matrix(0, nrow=NROW(aveVE), ncol=length(seq(0,1,by=0.01)))

  for (j in 1:NROW(aveVE)){
    trialDataCens.filename <- paste0("trialDataCens_nPlac=",N.pla,"_nVacc=",paste(rep(N.vax,N.vax.arms), collapse="_"),"_aveVE=",paste(aveVE[j,], collapse="_"),"_infRate=",infRate,"_",estimand,".RData")
    quantiles[j,] <- cumProbTrialDuration(trialDataCens.filename, srcDir)
  }

  # convert into months
  quantiles <- quantiles/(52/12)

  par(mar=c(4.5,5,3,1.5), las=1, cex.axis=1.2, cex.lab=1.3, cex.main=1.2)
  plot(-10, 0, xlim=c(0,54), ylim=0:1, xaxt="n", yaxt="n", xlab=paste(N.vax.arms,"-Vaccine-Arm Trial Duration (Months)",sep=""), ylab="Cumulative Probability")
  axis(side=1, at=c(0,seq(10,50,by=10),54))
  axis(side=2, at=seq(0,1,by=0.2))
  axis(side=4, at=seq(0,1,by=0.2), labels=FALSE)
  segments(x0=18, y0=0, y1=0.5, lwd=2, col="gray30")
  segments(x0=24, y0=0, y1=0.65, lwd=2, col="gray30")
  segments(x0=30, y0=0, y1=0.8, lwd=2, col="gray30")

  linesCol <- c("black", "blue", "red", "green", "darkorange")[1:NROW(aveVE)]
  for (i in 1:NROW(aveVE)){ lines(quantiles[i,], seq(0,1,by=0.01), lwd=2, col=linesCol[i]) }

  aveVEcomb <- NULL
  for (k in 1:NROW(aveVE)){ aveVEcomb <- c(aveVEcomb,paste("(",paste(paste(aveVE[k,]*100, "%", sep=""), collapse=", "),")",sep="")) }

  legend(x=0.5, y=0.975, legend=aveVEcomb, col=linesCol, lwd=2, cex=1.2, title="Average\nVE(0-18)", bty="n")
  text(0, 0.2, paste0("N=",N.pla,":",paste(rep(N.vax,N.vax.arms),collapse=":")," pla:",paste(rep("vac",N.vax.arms),collapse=":"),"\n","18-month accrual; average 309 per month\n",infRate*100,"% annual placebo incidence\n5% annual dropout\nVE halved in first 6 months"), adj=0, cex=1.2)
  text(x=18, y=0.45, labels="Accrual\ncompleted", cex=1.2, pos=2)
  text(x=24, y=0.6, labels="Vaccinations through\nmonth 6 completed", cex=1.2, pos=2)
  text(x=30, y=0.75, labels="Vaccinations through\nmonth 12 completed", cex=1.2, pos=2)
}

\clearpage

Endpoint Accrual

aveVE <- 0.6
pMissVax <- 0.02
tableEventAccrual <- tabEventAccrual(file.path(srcDir, paste0("simTrial_nPlac=",N.pla,"_nVacc=",N.vax,"_aveVE=",aveVE,"_infRate=",format(infRate, digits=3, nsmall=3),".RData")),
                                     atEvents=sort(c(seq(30, 200, by=10), atEventCount)), lagTimeMITT=15 / 7, lagTimePP=(28 + 15) / 7, prob=0.5, na.ub=0.3)
# convert time variables to month
tableEventAccrual[, 2] <- tableEventAccrual[, 2] / (52 / 12)
tableEventAccrual[, 3] <- tableEventAccrual[, 3] / (52 / 12)

print(xtable(tableEventAccrual,
             digits=c(0, 0, 1, 1),
             align=rep("r", NCOL(tableEventAccrual) + 1),
             caption=paste0("Median time since the first vaccination to accrue given COV-DIS endpoint totals occurring more that 14 days after the first vaccination in the FAS cohort and those occurring more than 14 days after the last vaccination in the PP cohort, under $\\overline{VE}\\textrm{(2--32)} = ", aveVE*100, "\\%$ with the VE model described in the simulation setup ($N=",N.pla,"$ in the placebo group, $N=",N.vax,"$ in the vaccine group, 25 months of participant follow-up, ", format(infRate * 100, digits=2, nsmall=2), "\\% annual placebo incidence rate, 2\\% annual dropout rate, and ",pMissVax*100,"\\% probability of a missed second dose).")),
      NA.string="n/a",
      table.placement="H",
      caption.placement="top",
      include.rownames=FALSE,
      include.colnames=FALSE,
      #scalebox=0.9,
      hline.after=NULL,
      add.to.row=list(pos=list(0, NROW(tableEventAccrual)),
                      command=c("\\hline \\hline Endpoint & Med Months* & Med Months* \\\\\n Total & FAS$^{\\dagger}$ & PP$^{\\ddag}$\\\\\n \\hline",
                                "\\hline \\hline 
                                \\multicolumn{3}{l}{\\footnotesize * Median time (in months)} \\\\\n
                                \\multicolumn{3}{l}{\\footnotesize \\quad since first enrollment} \\\\\n
                                \\multicolumn{3}{l}{\\footnotesize $^{\\dagger}$ Endpoints in the FAS cohort starting} \\\\\n
                                \\multicolumn{3}{l}{\\footnotesize \\quad $>14$ days after enrollment} \\\\\n 
                                \\multicolumn{3}{l}{\\footnotesize $^{\\ddag}$ Endpoints in the PP cohort starting} \\\\\n
                                \\multicolumn{3}{l}{\\footnotesize \\quad $>14$ days after completed immunization}")), comment=FALSE)

\clearpage

aveVE <- 0.6
pMissVax <- 0.02
tableEventAccrual <- tabEventAccrual(file.path(srcDir, paste0("simTrial_nPlac=",N.pla,"_nVacc=",N.vax,"_aveVE=",aveVE,"_infRate=",format(infRate, digits=3, nsmall=3),".RData")),
                                     atEvents=sort(c(seq(30, 200, by=10), atEventCount)), lagTimeMITT=15 / 7, lagTimePP=(28 + 15) / 7, prob=0.8, na.ub=0.3)
# convert time variables to month
tableEventAccrual[, 2] <- tableEventAccrual[, 2] / (52 / 12)
tableEventAccrual[, 3] <- tableEventAccrual[, 3] / (52 / 12)

print(xtable(tableEventAccrual,
             digits=c(0, 0, 1, 1),
             align=rep("r", NCOL(tableEventAccrual) + 1),
             caption=paste0("80$^{\\textrm{th}}$ percentile of the distribution of time since the trial start to accrue given COV-DIS endpoint totals occurring more that 14 days after the first vaccination in the FAS cohort and those occurring more than 14 days after the last vaccination in the PP cohort, under $\\overline{VE}\\textrm{(2--32)} = ", aveVE*100, "\\%$ with the VE model described in the simulation setup  ($N=",N.pla,"$ in the placebo group, $N=",N.vax,"$ in the vaccine group, 25 months of participant follow-up, ", format(infRate * 100, digits=2, nsmall=2), "\\% annual placebo incidence rate, 2\\% annual dropout rate, and ",pMissVax*100,"\\% probability of a missed second dose).")),
      NA.string="n/a",
      table.placement="H",
      caption.placement="top",
      include.rownames=FALSE,
      include.colnames=FALSE,
      #scalebox=0.9,
      hline.after=NULL,
      add.to.row=list(pos=list(0, NROW(tableEventAccrual)),
                      command=c("\\hline \\hline Endpoint & 80$^{\\textrm{th}}$ \\%ile Months* & 80$^{\\textrm{th}}$ \\%ile Months* \\\\\n Total & FAS$^{\\dagger}$ & PP$^{\\ddag}$\\\\\n \\hline",
                                "\\hline \\hline 
                                \\multicolumn{3}{l}{\\footnotesize * 80$^{\\textrm{th}}$ percentile of time (in months)} \\\\\n
                                \\multicolumn{3}{l}{\\footnotesize \\quad since first enrollment} \\\\\n
                                \\multicolumn{3}{l}{\\footnotesize $^{\\dagger}$ Endpoints in the FAS cohort starting} \\\\\n
                                \\multicolumn{3}{l}{\\footnotesize \\quad $>14$ days after enrollment} \\\\\n 
                                \\multicolumn{3}{l}{\\footnotesize $^{\\ddag}$ Endpoints in the PP cohort starting} \\\\\n
                                \\multicolumn{3}{l}{\\footnotesize \\quad $>14$ days after completed immunization}")), comment=FALSE)

\clearpage

aveVE <- 0.8
pMissVax <- 0.02
tableEventAccrual <- tabEventAccrual(file.path(srcDir, paste0("simTrial_nPlac=",N.pla,"_nVacc=",N.vax,"_aveVE=",aveVE,"_infRate=",format(infRate, digits=3, nsmall=3),".RData")),
                                     atEvents=sort(c(seq(30, 200, by=10), atEventCount)), lagTimeMITT=15 / 7, lagTimePP=(28 + 15) / 7, prob=0.5, na.ub=0.3)
# convert time variables to month
tableEventAccrual[, 2] <- tableEventAccrual[, 2] / (52 / 12)
tableEventAccrual[, 3] <- tableEventAccrual[, 3] / (52 / 12)

print(xtable(tableEventAccrual,
             digits=c(0, 0, 1, 1),
             align=rep("r", NCOL(tableEventAccrual) + 1),
             caption=paste0("Median time since the first vaccination to accrue given COV-DIS endpoint totals occurring more that 14 days after the first vaccination in the FAS cohort and those occurring more than 14 days after the last vaccination in the PP cohort, under $\\overline{VE}\\textrm{(2--32)} = ", aveVE*100, "\\%$ with the VE model described in the simulation setup ($N=",N.pla,"$ in the placebo group, $N=",N.vax,"$ in the vaccine group, 25 months of participant follow-up, ", format(infRate * 100, digits=2, nsmall=2), "\\% annual placebo incidence rate, 2\\% annual dropout rate, and ",pMissVax*100,"\\% probability of a missed second dose).")),
      NA.string="n/a",
      table.placement="H",
      caption.placement="top",
      include.rownames=FALSE,
      include.colnames=FALSE,
      #scalebox=0.9,
      hline.after=NULL,
      add.to.row=list(pos=list(0, NROW(tableEventAccrual)),
                      command=c("\\hline \\hline Endpoint & Med Months* & Med Months* \\\\\n Total & FAS$^{\\dagger}$ & PP$^{\\ddag}$\\\\\n \\hline",
                                "\\hline \\hline 
                                \\multicolumn{3}{l}{\\footnotesize * Median time (in months)} \\\\\n
                                \\multicolumn{3}{l}{\\footnotesize \\quad since first enrollment} \\\\\n
                                \\multicolumn{3}{l}{\\footnotesize $^{\\dagger}$ Endpoints in the FAS cohort starting} \\\\\n
                                \\multicolumn{3}{l}{\\footnotesize \\quad $>14$ days after enrollment} \\\\\n 
                                \\multicolumn{3}{l}{\\footnotesize $^{\\ddag}$ Endpoints in the PP cohort starting} \\\\\n
                                \\multicolumn{3}{l}{\\footnotesize \\quad $>14$ days after completed immunization}")), comment=FALSE)

\clearpage

aveVE <- 0.6
pMissVax <- 0.02
tableEventAccrual <- tabEventAccrual(file.path(srcDir, paste0("simTrial_nPlac=",N.pla,"_nVacc=",N.vax,"_aveVE=",aveVE,"_infRate=",format(infRate, digits=3, nsmall=3),".RData")),
                                     atWeeks=1:12 * 52 / 12, lagTimeMITT=15 / 7, lagTimePP=(28 + 15) / 7, prob=0.5, na.ub=0.3)
# convert time variables to month
tableEventAccrual[, 1] <- tableEventAccrual[, 1] / (52 / 12)

print(xtable(tableEventAccrual,
             digits=c(0, 0, 0, 0),
             align=rep("r", NCOL(tableEventAccrual) + 1),
             caption=paste0("Median COV-DIS endpoint totals occurring more that 14 days after the first vaccination in the FAS cohort and those occurring more than 14 days after the last vaccination in the PP cohort observed by a given time since the first vaccination, under $\\overline{VE}\\textrm{(2--32)} = ", aveVE*100, "\\%$ with the VE model described in the simulation setup  ($N=",N.pla,"$ in the placebo group, $N=",N.vax,"$ in the vaccine group, 25 months of participant follow-up, ", format(infRate * 100, digits=2, nsmall=2), "\\% annual placebo incidence rate, 2\\% annual dropout rate, and ",pMissVax*100,"\\% probability of a missed second dose).")),
      NA.string="n/a",
      table.placement="H",
      caption.placement="top",
      include.rownames=FALSE,
      include.colnames=FALSE,
      #scalebox=0.9,
      hline.after=NULL,
      add.to.row=list(pos=list(0, NROW(tableEventAccrual)),
                      command=c("\\hline \\hline  & Med Endpoints & Med Endpoints \\\\\n Months* & FAS$^{\\dagger}$ & PP$^{\\ddag}$\\\\\n \\hline",
                                "\\hline \\hline 
                                \\multicolumn{3}{l}{\\footnotesize * Months since first enrollment} \\\\\n
                                \\multicolumn{3}{l}{\\footnotesize $^{\\dagger}$ Median endpoints in the FAS cohort starting} \\\\\n
                                \\multicolumn{3}{l}{\\footnotesize \\quad $>14$ days after enrollment} \\\\\n 
                                \\multicolumn{3}{l}{\\footnotesize $^{\\ddag}$ Median endpoints in the PP cohort starting} \\\\\n
                                \\multicolumn{3}{l}{\\footnotesize \\quad $>14$ days after completed immunization}")), comment=FALSE)

\clearpage

aveVE <- 0.6
pMissVax <- 0.02
tableEventAccrual <- tabEventAccrual(file.path(srcDir, paste0("simTrial_nPlac=",N.pla,"_nVacc=",N.vax,"_aveVE=",aveVE,"_infRate=",format(infRate, digits=3, nsmall=3),".RData")),
                                     atWeeks=1:12 * 52 / 12, lagTimeMITT=15 / 7, lagTimePP=(28 + 15) / 7, prob=0.2, na.ub=0.3)
# convert time variables to month
tableEventAccrual[, 1] <- tableEventAccrual[, 1] / (52 / 12)

print(xtable(tableEventAccrual,
             digits=c(0, 0, 0, 0),
             align=rep("r", NCOL(tableEventAccrual) + 1),
             caption=paste0("20$^{\\textrm{th}}$ percentile of the distribution of COV-DIS endpoint totals occurring more that 14 days after the first vaccination in the FAS cohort and those occurring more than 14 days after the last vaccination in the PP cohort observed by a given time since the first vaccination, under $\\overline{VE}\\textrm{(2--32)} = ", aveVE*100, "\\%$ with the VE model described in the simulation setup ($N=",N.pla,"$ in the placebo group, $N=",N.vax,"$ in the vaccine group, 25 months of participant follow-up, ", format(infRate * 100, digits=2, nsmall=2), "\\% annual placebo incidence rate, 2\\% annual dropout rate, and ",pMissVax*100,"\\% probability of a missed second dose).")),
      NA.string="n/a",
      table.placement="H",
      caption.placement="top",
      include.rownames=FALSE,
      include.colnames=FALSE,
      #scalebox=0.9,
      hline.after=NULL,
      add.to.row=list(pos=list(0, NROW(tableEventAccrual)),
                      command=c("\\hline \\hline  & 20$^{\\textrm{th}}$ \\%ile Endpoints & 20$^{\\textrm{th}}$ \\%ile Endpoints \\\\\n Months* & FAS$^{\\dagger}$ & PP$^{\\ddag}$\\\\\n \\hline",
                                "\\hline \\hline 
                                \\multicolumn{3}{l}{\\footnotesize * Months since first enrollment} \\\\\n
                                \\multicolumn{3}{l}{\\footnotesize $^{\\dagger}$ 20$^{\\textrm{th}}$ percentile endpoints in the FAS cohort starting} \\\\\n
                                \\multicolumn{3}{l}{\\footnotesize \\quad $>14$ days after enrollment} \\\\\n 
                                \\multicolumn{3}{l}{\\footnotesize $^{\\ddag}$ 20$^{\\textrm{th}}$ percentile endpoints in the PP cohort starting} \\\\\n
                                \\multicolumn{3}{l}{\\footnotesize \\quad $>14$ days after completed immunization}")), comment=FALSE)

\clearpage

aveVE <- 0.6
pMissVax <- 0.02
infRate1 <- 0.5 * infRate
tableEventAccrual <- tabEventAccrual(file.path(srcDir, "moreEventAccrualScenarios", 
                                               paste0("simTrial_nPlac=",N.pla,"_nVacc=",N.vax,"_aveVE=",aveVE,"_infRate=",format(infRate1, digits=3, nsmall=3),".RData")),
                                     atWeeks=1:12 * 52 / 12, lagTimeMITT=15 / 7, lagTimePP=(28 + 15) / 7, prob=0.5, na.ub=0.3)
# convert time variables to month
tableEventAccrual[, 1] <- tableEventAccrual[, 1] / (52 / 12)

print(xtable(tableEventAccrual,
             digits=c(0, 0, 0, 0),
             align=rep("r", NCOL(tableEventAccrual) + 1),
             caption=paste0("Median COV-DIS endpoint totals occurring more that 14 days after the first vaccination in the FAS cohort and those occurring more than 14 days after the last vaccination in the PP cohort observed by a given time since the first vaccination, under $\\overline{VE}\\textrm{(2--32)} = ", aveVE*100, "\\%$ with the VE model described in the simulation setup  ($N=",N.pla,"$ in the placebo group, $N=",N.vax,"$ in the vaccine group, 25 months of participant follow-up, ", format(infRate1 * 100, digits=2, nsmall=2), "\\% annual placebo incidence rate, 2\\% annual dropout rate, and ",pMissVax*100,"\\% probability of a missed second dose).")),
      NA.string="n/a",
      table.placement="H",
      caption.placement="top",
      include.rownames=FALSE,
      include.colnames=FALSE,
      #scalebox=0.9,
      hline.after=NULL,
      add.to.row=list(pos=list(0, NROW(tableEventAccrual)),
                      command=c("\\hline \\hline  & Med Endpoints & Med Endpoints \\\\\n Months* & FAS$^{\\dagger}$ & PP$^{\\ddag}$\\\\\n \\hline",
                                "\\hline \\hline 
                                \\multicolumn{3}{l}{\\footnotesize * Months since first enrollment} \\\\\n
                                \\multicolumn{3}{l}{\\footnotesize $^{\\dagger}$ Median endpoints in the FAS cohort starting} \\\\\n
                                \\multicolumn{3}{l}{\\footnotesize \\quad $>14$ days after enrollment} \\\\\n 
                                \\multicolumn{3}{l}{\\footnotesize $^{\\ddag}$ Median endpoints in the PP cohort starting} \\\\\n
                                \\multicolumn{3}{l}{\\footnotesize \\quad $>14$ days after completed immunization}")), comment=FALSE)

\clearpage

aveVE <- 0.6
pMissVax <- 0.02
infRate1 <- 0.75 * infRate
tableEventAccrual <- tabEventAccrual(file.path(srcDir, "moreEventAccrualScenarios",
                                               paste0("simTrial_nPlac=",N.pla,"_nVacc=",N.vax,"_aveVE=",aveVE,"_infRate=",format(infRate1, digits=3, nsmall=3),".RData")),
                                     atWeeks=1:12 * 52 / 12, lagTimeMITT=15 / 7, lagTimePP=(28 + 15) / 7, prob=0.5, na.ub=0.3)
# convert time variables to month
tableEventAccrual[, 1] <- tableEventAccrual[, 1] / (52 / 12)

print(xtable(tableEventAccrual,
             digits=c(0, 0, 0, 0),
             align=rep("r", NCOL(tableEventAccrual) + 1),
             caption=paste0("Median COV-DIS endpoint totals occurring more that 14 days after the first vaccination in the FAS cohort and those occurring more than 14 days after the last vaccination in the PP cohort observed by a given time since the first vaccination, under $\\overline{VE}\\textrm{(2--32)} = ", aveVE*100, "\\%$ with the VE model described in the simulation setup  ($N=",N.pla,"$ in the placebo group, $N=",N.vax,"$ in the vaccine group, 25 months of participant follow-up, ", format(infRate1 * 100, digits=2, nsmall=2), "\\% annual placebo incidence rate, 2\\% annual dropout rate, and ",pMissVax*100,"\\% probability of a missed second dose).")),
      NA.string="n/a",
      table.placement="H",
      caption.placement="top",
      include.rownames=FALSE,
      include.colnames=FALSE,
      #scalebox=0.9,
      hline.after=NULL,
      add.to.row=list(pos=list(0, NROW(tableEventAccrual)),
                      command=c("\\hline \\hline  & Med Endpoints & Med Endpoints \\\\\n Months* & FAS$^{\\dagger}$ & PP$^{\\ddag}$\\\\\n \\hline",
                                "\\hline \\hline 
                                \\multicolumn{3}{l}{\\footnotesize * Months since first enrollment} \\\\\n
                                \\multicolumn{3}{l}{\\footnotesize $^{\\dagger}$ Median endpoints in the FAS cohort starting} \\\\\n
                                \\multicolumn{3}{l}{\\footnotesize \\quad $>14$ days after enrollment} \\\\\n 
                                \\multicolumn{3}{l}{\\footnotesize $^{\\ddag}$ Median endpoints in the PP cohort starting} \\\\\n
                                \\multicolumn{3}{l}{\\footnotesize \\quad $>14$ days after completed immunization}")), comment=FALSE)

\clearpage

Endpoints in Vaccine Arm for Correlates Evaluation

infThroughWeekVec <- c(7, 13) * 52 / 12
pMissVax <- 0.02

for (infThroughWeek in infThroughWeekVec){
  aveVE <- 0.9
  p <- c(0.01, 0.05, 0.25, 0.5, 0.75, 0.95, 0.99)
  ppname <- "pp1"

  res <- as.data.frame(matrix(0, nrow=2, ncol=length(p)+1))

  trialDataCens.filename <- paste0("trialDataCens_nPlac=",N.pla,"_nVacc=",N.vax,"_aveVE=",aveVE,"_infRate=",format(infRate, digits=3, nsmall=3),"_",estimand,".RData")
  load(file.path(srcDir, trialDataCens.filename))

  # for each trial, create a vector of Stage 1 infection counts by treatment
  infecListMITT <- infecListPP <- as.list(NULL)
  for (i in 1:length(trialListCensor)){
    dataI <- trialListCensor[[i]]
    dataI$futime <- dataI$exit - dataI$entry
    infecListMITT[[i]] <- as.vector( table( as.factor(dataI$trt)[dataI$event & dataI$futime>=15 / 7 & dataI$futime<=infThroughWeek] ) )
    infecListPP[[i]] <- as.vector( table( as.factor(dataI$trt)[dataI[[ppname]]==1 & dataI$event & dataI$futime>=(28 + 15) / 7 & dataI$futime<=infThroughWeek] ) )
  }

  # for each trial, calculate the number of infections diagnosed by infThroughWeek weeks among vaccine recipients
  NinfMITT <- sapply(infecListMITT, function(vectInfbyTx){ sum(vectInfbyTx[-1], na.rm=TRUE) })
  NinfPP <- sapply(infecListPP, function(vectInfbyTx){ sum(vectInfbyTx[-1], na.rm=TRUE) })

  res[1,1] <- mean(NinfMITT)
  res[1,-1] <- quantile(NinfMITT, prob=p)    
  res[2,1] <- mean(NinfPP)
  res[2,-1] <- quantile(NinfPP, prob=p)  

  res <- round(res, 0)
  colnames(res) <- c("Mean", paste(p*100,"%",sep=""))
  print(xtable(res,
               digits=0,
               align=rep("c",length(p)+2),
               caption=paste0("Distribution of the number of month 0.5--",infThroughWeek*12/52," FAS and month 1.5--",infThroughWeek*12/52," PP COV-DIS endpoints in the vaccine group, for scenarios with $\\overline{VE}\\textrm{(2--32)}=",aveVE*100,"\\%$ and the VE model described in the simulation setup ($N=",N.pla,"$ in the placebo group, $N=",N.vax,"$ in the vaccine group, and per-protocol cohort is identified assuming ",pMissVax*100,"\\% probability of $\\geq 1$ missed vaccination).")
               ),
        table.placement="H",
        caption.placement="top",
        include.rownames=FALSE,
        include.colnames=FALSE,
        #scalebox=0.9,
        hline.after=NULL,
        add.to.row=list(pos=list(-1,0,1,2), 
                        command=c(paste0("\\hline \\hline & \\multicolumn{7}{c}{Percentiles of distribution of \\# of} \\\\\n & \\multicolumn{7}{c}{COV-DIS endpoints in vax arm} \\\\\n \\cline{2-8}"),
                                  paste0("Mean & 1\\% & 5\\% & 25\\% & 50\\% & 75\\% & 95\\% & 99\\% \\\\\n  \\hline \\multicolumn{8}{c}{Month 0.5--",infThroughWeek*12/52," COV-DIS endpoints in the FAS cohort} \\\\\n"),
                                  paste0("\\hline \\multicolumn{8}{c}{Month 1.5--",infThroughWeek*12/52," COV-DIS endpoints in PP* cohort} \\\\\n"),
                                  paste0("\\hline \\hline \\multicolumn{8}{l}{\\footnotesize N=",N.pla,":",N.vax," Placebo:vaccine arm} \\\\\n \\multicolumn{8}{l}{\\footnotesize ", format(infRate * 100, digits=2, nsmall=2),"\\% annual incidence rate in placebo arm} \\\\\n \\multicolumn{8}{l}{\\footnotesize $\\overline{VE}\\textrm{(2--32)} = ", aveVE*100, "\\%$} \\\\\n \\multicolumn{8}{l}{\\footnotesize 2\\% dropout rate in each arm} \\\\\n \\multicolumn{8}{l}{\\footnotesize *",pMissVax*100,"\\% probability of a missed second dose}"))),
        comment=FALSE
        )
}

\clearpage

infThroughWeek <- 7 * 52 / 12
pMissVax <- 0.02
p <- 0.5
ppname <- "pp1"

aveVEvector <- c(0, seq(0.4, 0.8, by=0.1))

medNfailuresMITT <- medNvaxFailuresMITT <- medNplaFailuresMITT <- NULL
for (aveVE in aveVEvector){

  trialDataCens.filename <- paste0("trialDataCens_nPlac=",N.pla,"_nVacc=",N.vax,"_aveVE=",aveVE,"_infRate=",format(infRate, digits=3, nsmall=3),"_",estimand,".RData")
  load(file.path(srcDir, trialDataCens.filename))

  # for each trial, create a vector of Stage 1 infection counts by treatment
  infecListMITT <- as.list(NULL)
  for (i in 1:length(trialListCensor)){
    dataI <- trialListCensor[[i]]
    dataI$futime <- dataI$exit - dataI$entry
    infecListMITT[[i]] <- as.vector( table( as.factor(dataI$trt)[dataI$event & dataI$futime>=15 / 7 & dataI$futime<=infThroughWeek] ) )
  }

  # for each trial, calculate the number of infections diagnosed by infThroughWeek weeks among vaccine recipients
  nVaxFailuresMITT <- sapply(infecListMITT, function(vectInfbyTx){ sum(vectInfbyTx[-1], na.rm=TRUE) })

  # for each trial, calculate the number of infections diagnosed by infThroughWeek weeks among placebo recipients
  nPlaFailuresMITT <- sapply(infecListMITT, function(vectInfbyTx){ sum(vectInfbyTx[1], na.rm=TRUE) })

  # for each trial, calculate the total number of infections diagnosed by infThroughWeek weeks
  nFailuresMITT <- sapply(infecListMITT, function(vectInfbyTx){ sum(vectInfbyTx, na.rm=TRUE) })

  medNvaxFailuresMITT <- c(medNvaxFailuresMITT, quantile(nVaxFailuresMITT, prob=p))
  medNplaFailuresMITT <- c(medNplaFailuresMITT, quantile(nPlaFailuresMITT, prob=p))
  medNfailuresMITT <- c(medNfailuresMITT, quantile(nFailuresMITT, prob=p))
}

LB95 <- 1 - exp(log(1 - aveVEvector) + qnorm(0.975) * sqrt(9 / (2 * medNfailuresMITT)))
UB95 <- 1 - exp(log(1 - aveVEvector) - qnorm(0.975) * sqrt(9 / (2 * medNfailuresMITT)))

tableMITTfailureSplitAndCI <- data.frame(aveVE=round(aveVEvector * 100, digits=0), splitMITT=paste0(round(medNvaxFailuresMITT, digits=0), ":", round(medNplaFailuresMITT, digits=0)), 
                                         LB95=round(LB95 * 100, digits=0), UB95=round(UB95 * 100, digits=0))

print(xtable(tableMITTfailureSplitAndCI,
             align=rep("r", 5),
             digits=rep(0, 5),
             caption="Month 0.5-7 COV-DIS primary endpoint splits in the FAS cohort and 95\\% confidence bounds for VE at the time of stopping the trial"
             ),
      math.style.negative=TRUE,
      include.rownames=FALSE,
      include.colnames=FALSE,
      caption.placement="top",
      hline.after=NULL,
      add.to.row=list(pos=list(0, NROW(tableMITTfailureSplitAndCI)), 
                      command=c("\\hline \\hline Ave & & \\multicolumn{2}{c}{VE} \\\\\n \\cline{3-4} VE(0--6) & Vaccine: & 95\\% LB$^{\\dagger}$ & 95\\% UB$^{\\dagger}$ \\\\\n (\\%) & Placebo* & (\\%) & (\\%) \\\\\n \\hline", 
                                "\\hline \\hline \\multicolumn{4}{l}{\\footnotesize $^{\\ast}$ Median endpoint counts} \\\\\n \\multicolumn{4}{l}{\\footnotesize $^{\\dagger}$ Under proportional hazards}")),
      comment=FALSE
      )

\clearpage

aveVE <- c(0, 0.6)
p <- c(0.01, 0.025, 0.05, seq(0.1,0.9,by=0.1), 0.95, 0.975, 0.99)

res <- as.data.frame(matrix(0, nrow=length(aveVE), ncol=length(p)))

for (i in 1:length(aveVE)){
  simTrial.filename <- paste0("simTrial_nPlac=",N.pla,"_nVacc=",N.vax,"_aveVE=",aveVE[i],"_infRate=",format(infRate, digits=3, nsmall=3),".RData")
  load(file.path(srcDir, simTrial.filename))
  Ninf <- sapply(trialObj$NinfStage1, function(vectInfbyTx){ vectInfbyTx[1] + max(vectInfbyTx[-1], na.rm=TRUE) })
  res[i,] <- quantile(Ninf, prob=p)
}

res <- round(res, 0)
res <- cbind(c("0%","60%"), res)
colnames(res) <- c("VE (0-6)", paste(p*100,"%",sep=""))
print(xtable(res,
             digits=0,
             align=c("c","l",rep("c",15)),
             caption=paste("Distribution of the number of Month 0--6 endpoints pooled over the control vaccine group and the experimental vaccine group, ignoring sequential monitoring for non-efficacy (n=",N.pla," in the control vaccine group and n=",N.vax," in the experimental vaccine group)",sep="")
             ),
      table.placement="H",
      caption.placement="top",
      include.rownames=FALSE,
      include.colnames=FALSE,
      scalebox=0.85,
      hline.after=0,
      add.to.row=list(pos=list(-1,0,NROW(res)), command=c("\\hline \\hline & \\multicolumn{15}{c}{Percentiles of distribution of number of Month 0--24 endpoints} \\\\\n","VE(0-6) & 1\\% & 2.5\\% & 5\\% & 10\\% & 20\\% & 30\\% & 40\\% & 50\\% & 60\\% & 70\\% & 80\\% & 90\\% & 95\\% & 97.5\\% & 99\\% \\\\\n",paste0("\\hline \\hline \\multicolumn{16}{l}{\\footnotesize N=",N.pla,":",paste(rep(N.vax,N.vax.arms),collapse=":")," control:",paste(rep("experimental vaccine",N.vax.arms),collapse=":")," group} \\\\\n \\multicolumn{16}{l}{\\footnotesize 3-month accrual} \\\\\n \\multicolumn{16}{l}{\\footnotesize ",infRate*100,"\\% annual incidence in the control vaccine group} \\\\\n \\multicolumn{16}{l}{\\footnotesize 5\\% annual dropout}"))),
      comment=FALSE
      )


  aveVE <- c(0, 0.6)
  p <- c(0.01, 0.025, 0.05, seq(0.1,0.9,by=0.1), 0.95, 0.975, 0.99)

  res <- as.data.frame(matrix(0, nrow=length(aveVE), ncol=length(p)))

  for (j in 1:length(aveVE)){
  #read in the censored data
    trialDataCens.filename <- paste0("trialDataCens_nPlac=",N.pla,"_nVacc=",N.vax,"_aveVE=",aveVE[j],"_infRate=",format(infRate, digits=3, nsmall=3),"_",estimand,".RData")
    load(file.path(srcDir, trialDataCens.filename))

    # for each trial, create a vector of Stage 1 infection counts by treatment
    infecList <- as.list(NULL)
    for (i in 1:length(trialListCensor)){
      dataI <- trialListCensor[[i]]
      dataI$futime <- dataI$exit - dataI$entry
      stage1ind <- dataI$event & dataI$futime<=26
      infecList[[i]] <- as.vector( table( as.factor(dataI$trt)[stage1ind] ) )
    }

    # for each trial, calculate the number of Stage 1 endpoints pooled over all treatment arms
    Ninf.all <- sapply(infecList, function(vectInfbyTx){ sum(vectInfbyTx, na.rm=TRUE) })

    # for each trial, calculate the number of Stage 1 endpoints pooled over the placebo arm and the vaccine arm with
    # the maximum number of endpoints
    Ninf <- sapply(infecList, function(vectInfbyTx){ vectInfbyTx[1] + max(vectInfbyTx[-1], na.rm=TRUE) })

    res[j,] <- quantile(Ninf, prob=p)
  }

  res <- round(res, 0)
  res <- cbind(c("0%","60%"), res)
  colnames(res) <- c("VE (0-6)", paste(p*100,"%",sep=""))
  print(xtable(res,
               digits=0,
               align=c("c","l",rep("c",15)),
               caption=paste("Distribution of the number of Month 0--6 endpoints pooled over the control vaccine group and the experimental vaccine group, accounting for sequential monitoring for non-efficacy (n=",N.pla," in the control vaccine group and n=",N.vax," in the experimental vaccine group)",sep="")
               ),
        table.placement="H",
        caption.placement="top",
        include.rownames=FALSE,
        include.colnames=FALSE,
        scalebox=0.85,
        hline.after=NULL,
        add.to.row=list(pos=list(-1,0,NROW(res)), command=c("\\hline \\hline & \\multicolumn{15}{c}{Percentiles of distribution of number of Month 0--6 endpoints} \\\\\n","VE(0-6) & 1\\% & 2.5\\% & 5\\% & 10\\% & 20\\% & 30\\% & 40\\% & 50\\% & 60\\% & 70\\% & 80\\% & 90\\% & 95\\% & 97.5\\% & 99\\% \\\\\n  \\hline", paste0("\\hline \\hline \\multicolumn{16}{l}{\\footnotesize N=",N.pla,":",paste(rep(N.vax,N.vax.arms),collapse=":")," control:",paste(rep("experimental  vaccine",N.vax.arms),collapse=":")," group} \\\\\n \\multicolumn{16}{l}{\\footnotesize 3-month accrual} \\\\\n \\multicolumn{16}{l}{\\footnotesize ",infRate*100,"\\% annual incidence in the control vaccine group} \\\\\n \\multicolumn{16}{l}{\\footnotesize 5\\% dropout rate in each group} \\\\\n \\multicolumn{16}{l}{\\footnotesize Cumulative incidence-based non-efficacy monitoring}"))),
        comment=FALSE
        )

# 'VEest' solves the equation 'log HRest - 1.96 * sqrt(4/d) = log(1 - UB)' for log HRest and then outputs the transformation 1 - exp(log HRest)
# 1:1 randomization ratio
VEest <- function(n, UB){ 1 - exp(log(1 - UB) + qnorm(0.975) * sqrt(9 / (2 * n))) }

# 'nPlaAtBoundary' computes the number of endpoints in the control arm at the boundary for a given total number of endpoints, a VE point estimate at the boundary,
# and the number at risk in each arm
# a binomial model is used, i.e., VE = 1 - ((n-n0)/N1) / (n0/N0)
nPlaAtBoundary <- function(n, VEpointEst, Npla, Nvax){ round((Npla * n) / (Nvax * (1 - VEpointEst) + Npla), digits=0) }

nPP <- c(66, 99, 131)
#nPP <- nITT - (nPlaCasesByTaumax + nVaxCasesByTaumax)

# number of placebo recipients who we expect to be in follow-up and at risk for the endpoint at 1.5 months
nPlaControlsByTaumax <- computeN(Nrand=N.pla, tau=0, taumax=1.5, VEtauToTaumax=0, VE0toTau=0, risk0=(0.021 / 0.9) * 1.5 / 12, dropoutRisk=0.02 * 1.5 / 12, propCasesWithS=1)$nControls

# number of vaccine recipients who we expect to be in follow-up and at risk for the endpoint at 1.5 months under VE(0-1.5) = 0%
nVaxControlsByTaumax <- computeN(Nrand=N.vax, tau=0, taumax=1.5, VEtauToTaumax=0, VE0toTau=0, risk0=(0.021 / 0.9) * 1.5 / 12, dropoutRisk=0.02 * 1.5 / 12, propCasesWithS=1)$nControls

VEestPP <- VEest(nPP, 0.5)

nPlaAtBoundaryPP <- nPlaAtBoundary(nPP, VEestPP, nPlaControlsByTaumax, nVaxControlsByTaumax)

# nPlaCasesByTaumax <- computeN(Nrand=N.pla, tau=0.5, taumax=1.5, VEtauToTaumax=0, VE0toTau=0, risk0=0.015 * 1 / 12, dropoutRisk=0.02 * 1 / 12, propCasesWithS=1)$nCases
# # number of ITT endpoints by taxmax in the vaccine arm under VE(0-1.5) = 0%
# nVaxCasesByTaumax <- computeN(Nrand=N.vax, tau=0.5, taumax=1.5, VEtauToTaumax=0, VE0toTau=0, risk0=0.015 * 1 / 12, dropoutRisk=0.02 * 1 / 12, propCasesWithS=1)$nCases
# 
# nITT <- nPP + nPlaCasesByTaumax + nVaxCasesByTaumax
# #nITT <- c(52, 77, 103)
# 
# VEestITT <- VEest(nITT, 0.5)
# 
# nPlaAtBoundaryITT <- nPlaAtBoundary(nITT, VEestITT, N.pla, N.vax)

tableNoneffBoundary <- data.frame(nPP=nPP, splitPP=paste0(nPP - nPlaAtBoundaryPP, ":", nPlaAtBoundaryPP), VEestPP=round(VEestPP * 100, digits=0))

# # non-eficacy monitoring starts at the fixed number of 43 PP endpoints (71 expected ITT endpoints)
# 
# # number of Rotarix recipients who we expect to be in follow-up and at risk for the endpoint at 2.5 months (2983)
# computeNval <- computeN(Nrand=N.pla, tau=0, taumax=2.5, VEtauToTaumax=0, VE0toTau=0, risk0=pexp(2.5, rate=-log(1 - 0.017) / 20), 
#                         dropoutRisk=pexp(2.5, rate=-log(1 - 0.05) / 22.5), propCasesWithS=1)
# nPlaTau <- computeNval$nCases + computeNval$nControls
# 
# # number of subunit vaccine recipients who we expect to be in follow-up and at risk for the endpoint at 2.5 months under RVE=-0.25% (5966)
# computeNval <- computeN(Nrand=N.vax, tau=0, taumax=2.5, VEtauToTaumax=-0.25, VE0toTau=-0.25, risk0=pexp(2.5, rate=-log(1 - 0.017) / 20), 
#                         dropoutRisk=pexp(2.5, rate=-log(1 - 0.05) / 22.5), propCasesWithS=1)
# nVaxTau <- computeNval$nCases + computeNval$nControls
# 
# # number of Rotarix PP endpoints required for binomial RVE=-70% given the total number of PP endpoints is 43 (10)
# n0 <- round((nPlaTau * 43) / (nVaxTau * (1 + 0.7) + nPlaTau), digits=0)
# 1 - ((43-n0)/nVaxTau)/(n0/nPlaTau)
# 
# # get RVE point estimate for a given number of endpoints and given 95% upper confidence bound (RVE=-47%)
# RVEest <- 1 - exp(log(1 - 0.1) + qnorm(0.975) * sqrt(9 / (2 * 71)))
# 
# # number of Rotarix ITT endpoints required for binomial RVE=-45% given the total number of ITT endpoints is 71 (10)
# n0 <- round((N.pla * 71) / (N.vax * (1 + abs(RVEest)) + N.pla), digits=0)
# 1 - ((71-n0)/N.vax)/(n0/N.pla)
# 
# # expected number of endpoints in the Rotarix arm within the last year among those who remain in follow-up (29)
# nPlaCasesWithinLastYear <- computeN(Nrand=N.pla, tau=10.5, taumax=22.5, VEtauToTaumax=0, VE0toTau=0, risk0=pexp(12, rate=-log(1 - 0.017) / 20),
#                                     dropoutRisk=pexp(12, rate=-log(1 - 0.05) / 22.5), propCasesWithS=1)$nCases
# 
# # expected number of endpoints in the subunit vaccine arm within the last year among those who remain in follow-up under RVE=-70% (99)
# nVaxCasesWithinLastYear <- computeN(Nrand=N.vax, tau=10.5, taumax=22.5, VEtauToTaumax=-0.7, VE0toTau=-0.7, risk0=pexp(12, rate=-log(1 - 0.017) / 20),
#                                     dropoutRisk=pexp(12, rate=-log(1 - 0.05) / 22.5), propCasesWithS=1)$nCases
# 
# # i.e., under RVE=-70%, we expect 43 + 128 = 171 total endpoints 1 year after the first non-efficacy analysis
# # get RVE point estimate for a given number of endpoints and given 95% upper confidence bound (RVE=-24%)
# RVEest <- 1 - exp(log(1 - 0.1) + qnorm(0.975) * sqrt(9 / (2 * 171)))
# 
# # number of Rotarix PP endpoints required for binomial RVE=-24% given the total number of PP endpoints is 171 (49)
# n0 <- round((nPlaTau * 171) / (nVaxTau * (1 + abs(RVEest)) + nPlaTau), digits=0)
# 1 - ((171-n0)/nVaxTau)/(n0/nPlaTau)
# 
# # get RVE point estimate for a given number of endpoints and given 95% upper confidence bound (RVE=-47%)
# RVEest <- 1 - exp(log(1 - 0.1) + qnorm(0.975) * sqrt(9 / (2 * 199)))
# 
# # number of Rotarix ITT endpoints required for binomial RVE=-45% given the total number of ITT endpoints is 71 (10)
# n0 <- round((N.pla * 71) / (N.vax * (1 + abs(RVEest)) + N.pla), digits=0)
# 1 - ((71-n0)/N.vax)/(n0/N.pla)

\clearpage

aveVEsets=matrix(c(0.4,0,
                   0.4,0.2,
                   0.4,0.3,
                   0.5,0.3,
                   0.5,0.4,
                   0.6,0.3,
                   0.6,0.4), nrow=7, ncol=2, byrow=TRUE)

typeVec <- "head"
for (type in typeVec){
  rankSelectPwVector <- headPwVector <- NULL
  for (i in 1:NROW(aveVEsets)){
    rankSelectPwr.filename <- paste0("rankSelectPwr_nPlac=",N.pla,"_nVacc=",paste(rep(N.vax,N.vax.arms), collapse="_"),"_aveVE=",paste(aveVEsets[i,], collapse="_"),"_infRate=",infRate,"_",estimand,".RData")
    load(file.path(srcDir, rankSelectPwr.filename))

    if (type=="head"){
      rankSelectPwVector <- c(rankSelectPwVector, out$rankSelectPw*100)  # power to correctly identify the winning efficacious vaccine
      headPwVector <- c(headPwVector, out$headHeadPw[1,1]*100)  # power to detect positive relative VE(0-stage1) of Vx1 vs Vx2
    }
    if (type=="pool"){
      rankSelectPwVector <- c(rankSelectPwVector, out$poolRankSelectPw*100)  # power to correctly identify the best pool
      headPwVector <- c(headPwVector, out$poolHeadPw[1,1]*100)  # power to detect positive relative VE(0-stage1) of, e.g., Ve3-Vx4 vs Vx1-Vx2
    }
  }
  aveVEchars <- apply(aveVEsets, 1, function(aveVE){ paste0("(",paste(aveVE*100,collapse=", "),")") })
  nRows <- length(aveVEchars)
  res <- data.frame(aveVEchars=aveVEchars, headPwVector=headPwVector, rankSelectPwVector=rankSelectPwVector)
  print(xtable(res,
               digits=1,
               align=rep("c",4),
               caption=ifelse(type=="head",
                                     "Power ($\\times 100$) to detect that relative VE(0--18) $> 0\\%$ comparing head-to-head vaccine regimens 1 vs. 2 and VE(0--18) $> 0\\%$ for vaccine regimen 1, and probability ($\\times 100$) of correct ranking and selection of the winning most efficacious vaccine regimen",
                                     "Power ($\\times 100$) to detect that relative VE(0--18) $> 0\\%$ comparing head-to-head pooled vaccine regimens 3--4 vs. 1--2 and VE(0--18) $> 0\\%$ for the pooled vaccine regimen 3--4, and probability ($\\times 100$) of correct ranking and selection among the pooled pairs of the winning most efficacious regimen"
                              )
               ),
        table.placement="H",
        caption.placement="top",
        include.rownames=FALSE,
        include.colnames=FALSE,
        #scalebox=0.9,        
        hline.after=0,
        add.to.row=list(pos=list(-1,0,nRows), command=c(paste0("\\hline \\hline & Power ($\\times 100$) & \\\\\n True average VE (\\%)$^1$ & ",ifelse(type=="head","Vx1 vs. Vx2","Vx3-4 vs. Vx1-2")," \\& & Probability ($\\times 100$) \\\\\n"),paste0("(Vx1, Vx2) & ",ifelse(type=="head","Vx1 vs. Placebo","Vx3-4 vs. Placebo"),"$^2$ & select best ",ifelse(type=="head","vaccine","pooled Vx"),"$^3$ \\\\\n"),paste0("\\hline\\hline \\multicolumn{3}{l}{\\footnotesize $^1$ VE halved in the first 6 months} \\\\\n \\multicolumn{3}{l}{\\footnotesize $^2$ Cumulative incidence-based Wald tests of both ",ifelse(type=="head","Vx1/Vx2","Vx3-4/Vx1-2")," and} \\\\\n \\multicolumn{3}{l}{\\footnotesize \\; \\,",ifelse(type=="head","Vx1/Placebo","Vx3-4/Placebo")," VE(0--18) with 1-sided $\\alpha=0.025$} \\\\\n \\multicolumn{3}{l}{\\footnotesize $^3$ Correct selection $=$ ",ifelse(type=="head","Vx1","pooled Vx3-4")," has highest estimated VE(0--36) and} \\\\\n \\multicolumn{3}{l}{\\footnotesize \\; \\, VE(0--18) significantly $>0\\%$} \\\\\n \\multicolumn{3}{l}{\\footnotesize N=",N.pla,":",paste(rep(N.vax,N.vax.arms),collapse=":")," pla:",paste(rep("vac",N.vax.arms),collapse=":")," group} \\\\\n \\multicolumn{3}{l}{\\footnotesize 18-month accrual; average 309 per month} \\\\\n \\multicolumn{3}{l}{\\footnotesize ",infRate*100,"\\% annual incidence in the placebo group} \\\\\n \\multicolumn{3}{l}{\\footnotesize 5\\% annual dropout} \\\\\n \\multicolumn{3}{l}{\\footnotesize Cumulative incidence-based monitoring}")))
        )
}

\clearpage

harmData <- read.csv(file.path(srcDir, paste0("harmBounds_N=60_alphaPerTest=0.0140_pVacc=",round(N.vax/(N.pla+N.vax),2),".csv")))
harmData <- subset(harmData, N>=10)
par(mar=c(5.5,5,1,1), las=1, cex.axis=1.1, cex.lab=1.2)
{plot(-10,0, xlim=c(10,60), ylim=c(9,40), xaxt="n", yaxt="n", xlab="", ylab="", type="n")
mtext("Total Number of Infections\n(Placebo + 1 Vaccine Arm)", side=1, line=3.5, cex=1.2)
mtext("Number of Vaccine Infections\n(1 Vaccine Arm)", side=2, line=2.5, cex=1.2, las=0)
axis(side=1, at=seq(10,60,by=5))
axis(side=2, at=seq(10,40,by=5))
axis(side=3, at=seq(10,60,by=5), labels=FALSE)
axis(side=4, at=seq(10,40,by=5), labels=FALSE)
abline(0,1, col="gray30")
text(28,32,"Y=X", col="gray30")
points(harmData$N, harmData$V, pch=19, col="darkred", cex=0.5)}
for (i in c(20,40,60)){
  segments(x0=i, y0=8, y1=harmData$V[harmData$N==i], lwd=2, lty="dashed", col="darkorange")
  segments(x0=10, x1=i, y0=harmData$V[harmData$N==i], lwd=2, lty="dashed", col="darkorange")
}
text(42, 21, "Potential Harm\nStopping\nBoundaries", adj=0, cex=1.1, font=2, col="blue")

VEbinomModel <- function(p, r){ 1 - r*p/(1-p) }

UL95 <- function(pHat, n, r, bound){ VEbinomModel(pHat, r) + qnorm(0.975)*sqrt(pHat*(r^2)/(n*((1-pHat)^3))) - bound }

LL95 <- function(pHat, n, r, bound){ VEbinomModel(pHat, r) - qnorm(0.975)*sqrt(pHat*(r^2)/(n*((1-pHat)^3))) - bound }

getpHatNonEff <- function(nVec, randRatioPV, upperBound){
  return(sapply(nVec, function(n){ uniroot(UL95, interval=c(0.01,0.9), n=n, r=randRatioPV, bound=upperBound)$root }))
}

getpHatHighEff <- function(nVec, randRatioPV, lowerBound){
  return(sapply(nVec, function(n){ uniroot(LL95, interval=c(0.01,0.9), n=n, r=randRatioPV, bound=lowerBound)$root }))
}

randRatio <- N.pla/N.vax

# calculate the MC average number of infections (pooled over placebo + 1 vaccine arm) triggering the first non-efficacy IA for trials with average VE=20%
load(file.path(srcDir,paste0("monitorTrial_nPlac=",N.pla,"_nVacc=",N.vax,"_aveVE=0.2_infRate=",infRate.fmt,"_",estimand,".RData")))
firstNoneffCnt <- round(mean(do.call("c", sapply(out, function(trial){ trial[[1]]$firstNonEffCnt }))),0)
nInf <- seq(firstNoneffCnt,240,by=20)

pHat <- getpHatNonEff(nVec=nInf, randRatioPV=randRatio, upperBound=0.4)
nVax <- round(pHat*nInf,0)
nSplitVaxPla <- paste0(nVax,":",nInf-nVax)
VEhat <- VEbinomModel(pHat, randRatio)*100
tableInfSplits <- data.frame(nInf, nSplitVaxPla, round(VEhat,0))

\clearpage

par(mar=c(5.5,5,1,1), las=1, cex.axis=1.1, cex.lab=1.2)
plot(-10,0, xlim=c(0,240), ylim=c(-100,20), xaxt="n", xlab="", ylab="", type="n")
mtext("Total Number of Infections by Month 18\n(Placebo + 1 Vaccine Arm)", side=1, line=3.5, cex=1.2)
mtext("Vaccine Efficacy (%)", side=2, line=2.5, cex=1.2, las=0)
axis(side=1, at=seq(0,240,by=20))
axis(side=2, at=seq(-25,20,by=5))
axis(side=3, at=seq(0,240,by=20), labels=FALSE)
axis(side=4, at=seq(-25,20,by=5), labels=FALSE)
for (i in 1:length(Ninf)){
  segments(x0=Ninf[i], y0=-30, y1=VEhat[i], lwd=2, col="darkred")
}
points(Ninf, VEhat, pch=19, col="darkred", cex=0.6)
lines(Ninf, VEhat, col="darkred")
text(5, 15, "Non-Efficacy\nStopping Boundaries", adj=0, cex=1.1, font=2, col="blue")

print(xtable(tableNoneffBoundary,
             align=rep("r", 4),
             digits=rep(0, 4),
             caption="COV-DIS endpoint splits occurring $\\geq 15$ days after the last vaccination in the PP cohort, and point estimates of VE at the non-efficacy boundary using Freidlin et al.'s nominal CI approach"
             ),
      math.style.negative=TRUE,
      include.rownames=FALSE,
      include.colnames=FALSE,
      caption.placement="top",
      comment=FALSE,
      hline.after=NULL,
      add.to.row=list(pos=list(0, NROW(tableNoneffBoundary)), 
                      command=c("\\hline \\hline Total & Vaccine: & Est. VE \\\\\n 
                                Endpoints$^{\\dagger}$ & Placebo$^{\\ast}$ & (\\%) \\\\\n 
                                \\hline ", 
                                "\\hline \\hline 
                                \\multicolumn{3}{l}{\\footnotesize $^{\\ast}$ Under $\\mathrm{VE} := 1 - \\textrm{vax-arm attack rate} \\big/$} \\\\\n
                                \\multicolumn{3}{l}{\\footnotesize \\quad $ \\big/ \\textrm{pla-arm attack rate}$} \\\\\n 
                                \\multicolumn{3}{l}{\\footnotesize $^{\\dagger}$ Scenario with $\\textrm{VE} = 0\\%$}"))
      )

mjuraska/seqDesign documentation built on Dec. 14, 2022, 4:26 p.m.

rdrr.io home R language documentation Run R code online

CRAN packages Bioconductor packages R-Forge packages GitHub packages

Note that we can't provide technical support on individual packages. You should contact the package authors for that.

mjuraska/seqDesign
Simulation and Group-Sequential Monitoring of Randomized Treatment Efficacy Trials with Time-to-Event Endpoints

In mjuraska/seqDesign: Simulation and Group-Sequential Monitoring of Randomized Treatment Efficacy Trials with Time-to-Event Endpoints

Simulation Set-up

Power Calculations

Interim Monitoring

Trial Duration

Endpoint Accrual

Endpoints in Vaccine Arm for Correlates Evaluation

R Package Documentation

Browse R Packages

We want your feedback!

mjuraska/seqDesign Simulation and Group-Sequential Monitoring of Randomized Treatment Efficacy Trials with Time-to-Event Endpoints

In mjuraska/seqDesign: Simulation and Group-Sequential Monitoring of Randomized Treatment Efficacy Trials with Time-to-Event Endpoints

Simulation Set-up

Power Calculations

Interim Monitoring

Trial Duration

Endpoint Accrual

Endpoints in Vaccine Arm for Correlates Evaluation

R Package Documentation

Browse R Packages

We want your feedback!

mjuraska/seqDesign
Simulation and Group-Sequential Monitoring of Randomized Treatment Efficacy Trials with Time-to-Event Endpoints