R/monitorTrial.R
In mjuraska/seqDesign: Simulation and Group-Sequential Monitoring of Randomized Treatment Efficacy Trials with Time-to-Event Endpoints
################# Begin code for monitorTrial ################
## Some arguments described just below the argument list
#' Group Sequential Monitoring of Simulated Efficacy Trials for the Event of Potential Harm, Non-Efficacy, and High Efficacy
#'
#' \code{monitorTrial} applies a group sequential monitoring procedure to data-sets generated by \code{simTrial}, which may result in modification or termination of each simulated trial.
#'
#' @param dataFile if \eqn{\code{saveDir} = \code{NULL}}, a list returned by \code{simTrial}; otherwise a name (character string) of an \code{.RData} file created by \code{simTrial}
#' @param stage1 the final week of stage 1 in a two-stage fixed-follow-up trial
#' @param stage2 the final week of stage 2 in a two-stage fixed-follow-up trial, i.e., the maximum total follow-up time
#' @param harmMonitorRange a 2-component numeric vector specifying the range of the pooled number of events (pooled over the placebo and vaccine arm accruing events the fastest) over which the type I error rate, specified in \code{harmMonitorAlpha}, shall be spent (per vaccine arm). Note that \code{harmMonitorRange} does not specify a range over which potential-harm stopping boundaries will be computed; instead, it specifies when potential-harm monitoring will start, and the range over which the type I error rate \code{harmMonitorAlpha} will be spent. If \eqn{\code{nonEffStartMethod} = \code{"FKG"}} or \code{"fixed"}, then the second value is ignored and can be replaced with \code{NA}.
#' @param harmMonitorAlpha a numeric value (0.05 by default) specifying the overall type I error rate for potential-harm monitoring (per vaccine arm). To turn off potential-harm monitoring, set \code{harmMonitorAlpha} equal to 0.00001.
#' @param alphaPerTest a per-test nominal significance level for potential-harm monitoring. If \code{NULL} (default), a per-test significance level is calculated that yields the overall type I error rate of \code{harmMonitorAlpha} at the end of \code{harmMonitorRange}.
#' @param nonEffStartMethod a character string specifying the method used for determining when non-efficacy monitoring is to start. The default method of Freidlin, Korn, and Gray (2010) ("\code{FKG}") calculates the minimal pooled event count (pooled over the placebo and vaccine arm accruing events the fastest) such that a hazard-ratio-based VE point estimate of 0\% would result in declaring non-efficacy, i.e., the upper bound of the two-sided \eqn{(1-\code{alphaNoneff}) 100\%} confidence interval for VE based on the asymptotic variance of the log-rank statistic equals the non-efficacy threshold specified as component \code{upperVEnonEff} in the list \code{nonEffStartParams}. If this list component is left unspecified, the argument \code{upperVEnonEff} is used as the non-efficacy threshold. The alternative method ("\code{fixed}") starts non-efficacy monitoring at a fixed pooled event count (pooled over the placebo and vaccine arm accruing events the fastest) specified by component \code{N1} in the list \code{nonEffStartParams}.
#' @param nonEffStartParams a list with named components specifying parameters required by \code{nonEffStartMethod} (\code{NULL} by default)
#' @param nonEffIntervalUnit a character string specifying whether intervals between two adjacent non-efficacy interim analyses should be event-driven (default option "\code{counts}") or calendar time-driven (option "\code{time}")
#' @param nonEffInterval a numeric vector (a number of events or a number of weeks) specifying the timing of non-efficacy interim analyses. If a single numeric value is specified, then all interim looks are equidistant (in terms of the number of events or weeks), and the value specifies the constant increment of information or time between two adjacent interim looks. If a numeric vector with at least two components is specified, then, following the initial interim look, the timing of subsequent interim looks is determined by (potentially differential) increments of information or time specified by this vector.
#' @param nonEffCohorts a named list specifying all cohorts (for both timing and analysis) used in non-efficacy monitoring. The required list components characterize the 'timing cohort,' i.e., the cohort events in which determine the analysis timepoints in an event-driven design (components \code{times}, \code{timeUnit}, and \code{timingCohort}), and the analysis cohort(s), i.e., the cohort(s) in which inference is made about the null hypothesis of non-efficacy (component \code{cohort1} is required, and \code{cohort2}, etc. are optional for additional analysis cohorts if, e.g., non-efficacy monitoring is conducted in both the ITT and per-protocol cohorts). As for the timing cohort, \code{times} specifies the analysis timepoints in terms of event counts (\eqn{\code{timeUnit} = \code{"counts"}} is the only implemented option). \code{timingCohort} is a list characterizing the timing cohort by components \code{lagTime}, a lag time (in weeks) that controls event inclusion for timing (if no lag is desired, it can be set to \code{NULL} or 0), and \code{cohortInd}, a character string naming an indicator variable included in the data frames outputted by \code{\link{simTrial}}, which also controls event inclusion for the purpose of determining analysis timepoints. The other top-level list components \code{cohort1}, \code{cohort2}, etc. are each a list that must contain a component named \code{estimand}, which can take on values \code{"cox"}, \code{"cuminc"}, or \code{"combined"}. Optional list components in \code{cohort1}, \code{cohort2}, etc. are \code{lagTime}, \code{cohortInd}, \code{nullVE}, and \code{nominalAlphas}. \code{lagTime} specifies a lag time (in weeks) that controls event inclusion for the analysis cohort. If no lag time is desired, then this component can be ignored, set to \code{NULL}, or set to 0. \code{cohortInd} is a character string naming an indicator variable included in the data frames outputted by \code{\link{simTrial}} to be used to subset participants into the analysis cohort. This component allows inclusion of, e.g., a "per-protocol" variable (e.g., by setting \code{cohortInd} to \code{"pp1"}). \code{nullVE} specifies the one-sided null hypothesis as \eqn{H_0: VE \geq \code{nullVE}}. The rejection of \eqn{H_0} constitutes evidence for non-efficacy. \code{nominalAlphas} specifies nominal significance levels in a two-arm event-driven trial design.
#' @param effCohort a named list specifying the timing and analysis cohorts used in group-sequential efficacy monitoring (if part of the monitoring plan). The required list components characterize the 'timing cohort,' i.e., the cohort events in which determine the analysis timepoints in an event-driven design (components \code{times}, \code{timeUnit}, and \code{timingCohort}), and a single analysis cohort, i.e., the cohort in which inference is made about the null hypothesis of efficacy (components \code{estimand}, \code{lagTime}, \code{cohortInd}, \code{nullVE}, and \code{nominalAlphas}). As for the timing cohort, \code{times} specifies the analysis timepoints in terms of event counts (\eqn{\code{timeUnit} = \code{"counts"}} is the only implemented option). \code{timingCohort} is a list characterizing the timing cohort by components \code{lagTime}, a lag time (in weeks) that controls event inclusion for timing (if no lag is desired, it can be set to \code{NULL} or 0), and \code{cohortInd}, a character string naming an indicator variable included in the data frames outputted by \code{\link{simTrial}}, which also controls event inclusion for the purpose of determining analysis timepoints. As for the analysis cohort, \code{estimand} can take on values \code{"cox"}, \code{"cuminc"}, or \code{"combined"}. \code{lagTime} specifies a lag time (in weeks) that controls event inclusion for the analysis cohort. If no lag time is desired, then this component can be ignored, set to \code{NULL}, or set to 0. \code{cohortInd} is a character string naming an indicator variable included in the data frames outputted by \code{\link{simTrial}} to be used to subset participants into the analysis cohort. This component allows inclusion of, e.g., a "per-protocol" variable (e.g., by setting \code{cohortInd} to \code{"pp1"}). \code{nullVE} specifies the one-sided primary null hypothesis as \eqn{H_0: VE \leq \code{nullVE}}. The rejection of \eqn{H_0} constitutes evidence for clinically relevant efficacy. \code{nominalAlphas} specifies nominal significance levels in a two-arm event-driven trial design.
#' @param lowerVEnoneff specifies an additional criterion for declaring non-efficacy if the hypothesis test is based on Wald confidence interval(s). It requires that the lower bound of the two-sided Wald CI(s) for the VE estimand(s), at the confidence level determined by \code{nonEffCohorts$cohort1$nominalAlphas}, etc., lie(s) below \code{lowerVEnoneff} (typically set equal to 0). If \code{NULL} (default), this criterion is ignored.
#' @param highVE specifies a criterion for declaring high-efficacy: the lower bound of the two-sided \eqn{(1-\code{alphaHigh}) 100\%} confidence interval for the VE estimand lies above \code{highVE} (typically a number in the 0.5--1 range). To turn off high efficacy monitoring, set \code{highVE} equal to 1.
#' @param stage1VE specifies a criterion for advancement of a treatment's evaluation into Stage 2: the lower bound of the two-sided \eqn{(1-\code{alphaStage1}) 100\%} confidence interval for the VE estimand lies above \code{stage1VE} (typically set equal to 0)
#' @param lowerVEuncPower a numeric vector with each component specifying a one-sided null hypothesis \eqn{H_0: VE(0-\code{stage1}) \leq \code{lowerVEuncPower} \times 100\%}. Unconditional power (i.e., accounting for sequential monitoring) to reject each \eqn{H_0} is calculated, where the rejection region is defined by the lower bound of the two-sided \eqn{(1-\code{alphaUncPower}) 100\%} confidence interval for the VE estimand being above the respective component of \code{lowerVEuncPower} (typically values in the 0--0.5 range).
#' @param alphaHigh one minus the nominal confidence level of the two-sided confidence interval used for high efficacy monitoring
#' @param alphaStage1 one minus the nominal confidence level of the two-sided confidence interval used for determining whether a treatment's evaluation advances into Stage 2
#' @param alphaUncPower one minus the nominal confidence level of the two-sided confidence interval used to test one-sided null hypotheses \eqn{H_0: VE(0-\code{stage1}) \leq \code{lowerVEuncPower} \times 100\%} against alternative hypotheses \eqn{H_1: VE(0-\code{stage1}) > \code{lowerVEuncPower} \times 100\%}. The same nominal confidence level is applied for each component of \code{lowerVEuncPower}.
#' @param saveFile a character string specifying the name of the output \code{.RData} file. If \code{NULL} (default), a default file name will be used.
#' @param saveDir a character string specifying a path for \code{dataFile}. If supplied, the output is also saved as an \code{.RData} file in this directory; otherwise the output is returned as a list.
#' @param verbose a logical value indicating whether information on the output directory, file name, and monitoring outcomes should be printed out (default is \code{TRUE})
#'
#' @details All time variables use week as the unit of time. Month is defined as 52/12 weeks.
#'
#' Potential harm monitoring starts at the \code{harmMonitorRange[1]}-th infection pooled over the placebo group and the vaccine regimen that accrues infections the fastest. The potential harm analyses continue at each additional infection until the first interim analysis for non-efficacy. The monitoring is implemented with exact one-sided binomial tests of H0: \eqn{p \le p0} versus H1: \eqn{p > p0}, where \eqn{p} is the probability that an infected participant was assigned to the vaccine group, and \eqn{p0} is a fixed constant that represents the null hypothesis that an infection is equally likely to be assigned vaccine or placebo. Each test is performed at the same prespecified nominal/unadjusted alpha-level (\code{alphaPerTest}), chosen based on simulations such that, for each vaccine regimen, the overall type I error rate by the \code{harmMonitorRange[2]}-th arm-pooled infection (i.e., the probability that the potential harm boundary is reached when the vaccine is actually safe, \eqn{p = p0}) equals \code{harmMonitorAlpha}.
#'
#' Non-efficacy is defined as evidence that it is highly unlikely that the vaccine has a beneficial effect measured as VE(0--\code{stage1}) of \code{upperVEnoneff} x 100\% or more. The non-efficacy analyses for each vaccine regimen will start at the first infection (pooled over the vaccine and placebo arm) determined by \code{nonEffStartMethod}. Stopping for non-efficacy will lead to a reported two-sided (1-\code{alphaNoneff}) x 100\% CI for VE(0--\code{stage1}) with, optionally, the lower confidence bound below \code{lowerVEnoneff} and the upper confidence bound below \code{upperVEnoneff}, where \code{estimand} determines the choice of the VE(0--\code{stage1}) estimand. This approach is similar to the inefficacy monitoring approach of Freidlin, Korn, and Gray (2010). If \code{estimand = "combined"}, stopping for non-efficacy will lead to reported (1-\code{alphaNoneff}) x 100\% CIs for both VE parameters with, optionally, lower confidence bounds below \code{lowerVEnoneff} and upper confidence bounds below \code{upperVEnoneff}. If \code{laggedMonitoring = TRUE}, stopping for non-efficacy will lead to reported (1-\code{alphaNoneff}) x 100\% CIs for both VE(0--\code{stage1}) and VE(\code{lagTime}--\code{stage1}) with, optionally, lower confidence bounds below \code{lowerVEnoneff} and upper confidence bounds below \code{upperVEnoneff}.
#'
#' High efficacy monitoring allows early detection of a highly protective vaccine if there is evidence that VE(0--\code{stage2}) \eqn{>} \code{highVE} x 100\%. It is synchronized with non-efficacy monitoring during Stage 1, and a single high-efficacy interim analysis during Stage 2 is conducted halfway between the end of Stage 1 and the end of the trial. While monitoring for potential harm and non-efficacy restricts to \code{stage1} infections, monitoring for high efficacy counts all infections during \code{stage1} or \code{stage2}, given that early stopping for high efficacy would only be warranted under evidence for durability of the efficacy.
#'
#' The following principles and rules are applied in the monitoring procedure:
#' \itemize{
#' \item Exclude all follow-up data from the analysis post-unblinding (and include all data pre-unblinding).
#' \item The monitoring is based on modified ITT analysis, i.e., all subjects documented to be free of the study endpoint at baseline are included and analyzed according to the treatment assigned by randomization, ignoring how many vaccinations they received (only pre-unblinding follow-up included).
#' \item If a vaccine hits the harm boundary, immediately discontinue vaccinations and accrual into this vaccine arm, and unblind this vaccine arm (continue post-unblinded follow-up until the end of Stage 1 for this vaccine arm).
#' \item If a vaccine hits the non-efficacy boundary, immediately discontinue vaccinations and accrual into this vaccine arm, keep blinded and continue follow-up until the end of Stage 1 for this vaccine arm.
#' \item If and when the last vaccine arm hits the non-efficacy (or harm) boundary, discontinue vaccinations and accrual into this vaccine arm, and unblind (the trial is over, completed in Stage 1).
#' \item Stage 1 for the whole trial is over on the earliest date of the two events: (1) all vaccine arms have hit the harm or non-efficacy boundary; and (2) the last enrolled subject in the trial reaches the final \code{stage1} visit.
#' \item Continue blinded follow-up until the end of Stage 2 for each vaccine arm that reaches the end of \code{stage1} with a positive efficacy (as defined by \code{stage1VE}) or high efficacy (as defined by \code{highVE}) result.
#' \item If at least one vaccine arm reaches the end of \code{stage1} with a positive efficacy or high efficacy result, continue blinded follow-up in the placebo arm until the end of Stage 2.
#' \item Stage 2 for the whole trial is over on the earliest date of the two events: (1) all subjects in the placebo arm and each vaccine arm that registered efficacy or high efficacy in \code{stage1} have failed or been censored; and (2) all subjects in the placebo arm and each vaccine arm that registered efficacy or high efficacy in \code{stage1} have completed the final \code{stage2} visit.
#' }
#'
#' The above rules have the following implications:
#' \itemize{
#' \item If a vaccine hits the non-efficacy boundary but Stage 1 for the whole trial is not over, then one includes in the analysis all follow-up through the final \code{stage1} visit for that vaccine regimen, including all individuals accrued up through the date of hitting the non-efficacy boundary (which will be the total number accrued to this vaccine arm).
#' \item If a vaccine hits the harm boundary, all follow-up information through the date of hitting the harm boundary is included for this vaccine; no follow-up data are included after this date.
#' \item If and when the last vaccine arm hits the non-efficacy (or harm) boundary, all follow-up information through the date of hitting the non-efficacy (or harm) boundary is included for this vaccine; no follow-up data are included after this date.
#' }
#'
#' @return If \code{saveDir} (and, optionally \code{saveFile}) is specified, the output list (named \code{out}) is saved as an \code{.RData} file in \code{saveDir} (the path to \code{saveDir} is printed); otherwise it is returned. The output object is a list of length equal to the number of simulated trials, each of which is a list of length equal to the number of treatment arms, each of which is a list with (at least) the following components:
#' \itemize{
#' \item \code{boundHit}: a character string stating the monitoring outcome in this treatment arm, i.e., one of \code{"Harm"}, \code{"NonEffInterim"}, \code{"NonEffFinal"}, \code{"Eff"}, or \code{"HighEff"}. The first four outcomes can occur in Stage 1, whereas the last outcome can combine data over Stage 1 and Stage 2.
#' \item \code{stopTime}: the time of hitting a stopping boundary since the first subject enrolled in the trial
#' \item \code{stopInfectCnt}: the pooled number of infections at \code{stopTime}
#' \item \code{summObj}: a \code{data.frame} containing summary information from each non-/high efficacy interim analysis
#' \item \code{finalHRci}: the final CI for the hazard ratio, available if \code{estimand!="cuminc"} and there is at least 1 infection in each arm
#' \item \code{firstNonEffCnt}: the number of infections that triggered non-efficacy monitoring (if available)
#' \item \code{totInfecCnt}: the total number of \code{stage1} (\code{stage2} if \code{boundHit = "HighEff"}) infections
#' \item \code{totInfecSplit}: a table with the numbers of \code{stage1} (\code{stage2} if \code{boundHit = "HighEff"}) infections in the treatment and control arm
#' \item \code{lastExitTime}: the time between the first subject's enrollment and the last subject's exiting from the trial
#' }
#'
#' @references Freidlin B., Korn E. L., and Gray R. (2010), A general inefficacy interim monitoring rule for randomized clinical trials. \emph{Clinical Trials} 7(3):197-208.
#'
#' @examples
#' simData <- simTrial(N=c(1000, 1000), aveVE=c(0, 0.4),
#' VEmodel="half", vePeriods=c(1, 27, 79), enrollPeriod=78,
#' enrollPartial=13, enrollPartialRelRate=0.5, dropoutRate=0.05,
#' infecRate=0.06, fuTime=156, visitSchedule=seq(0, 156, by=4),
#' missVaccProb=0.05, VEcutoffWeek=26, nTrials=5,
#' stage1=78, randomSeed=300)
#'
#' ### trial design: fixed follow-up; cumulative incidence-based VE estimand;
#' no efficacy monitoring; non-efficacy monitoring using
#' the Freidlin et al. method; hypothesis tests based on Wald
#' confidence intervals
#' ### RIGHT NOW, THE BELOW CALL IS BROKEN
#' monitorData <- monitorTrial(dataFile=simData, stage1=78, stage2=156,
#' harmMonitorRange=c(10, NA), harmMonitorAlpha=0.05,
#' nonEffStartMethod="FKG", nonEffInterval=20,
#' nonEffCohorts=list(timeUnit="counts",
#' # right now 'timingCohort' is required
#' timingCohort=list(lagTime=0),
#' cohort1=list(estimand="cuminc",
#' nullVE=0.4,
#' nominalAlphas=0.025)),
#' # it appears that right now 'effCohort' must be specified
#' # even when there is no efficacy monitoring;
#' # this call of monitorTrial() still doesn't run
#' # because it requires event counts for timing
#' effCohort=list(timeUnit="counts",
#' timingCohort=list(lagTime=0)),
#' lowerVEnoneff=0, highVE=0.7, lowerVEuncPower=0,
#' alphaHigh=0.05, alphaUncPower=0.05)
#'
#' ### trial design: event-driven; hazard-based VE estimand;
#' harmonized efficacy and non-efficacy monitoring;
#' hypothesis tests using the score test in the Cox model
#' ### THE SCORE TEST OPTION NEEDS TO BE ADDED
#' monitorData <- monitorTrial(dataFile=simData, stage1=78, stage2=156,
#' harmMonitorRange=c(10, 50), harmMonitorAlpha=0.05,
#' nonEffCohorts=list(
#' times=c(50, 100, 150),
#' timeUnit="counts",
#' timingCohort=list(lagTime=0),
#' cohort1=list(estimand="cox",
#' nullVE=0.4,
#' nominalAlphas=c(0.0001, 0.0060, 0.0231))),
#' effCohort=list(times=c(50, 100, 150),
#' timeUnit="counts",
#' timingCohort=list(lagTime=0),
#' estimand="cox",
#' nullVE=0,
#' nominalAlphas=c(0.0001, 0.0060, 0.0231)),
#' highVE=1, lowerVEuncPower=0,
#' alphaHigh=0.05, alphaUncPower=0.05)
#'
#' ### alternatively, to save the .RData output file (no '<-' needed):
#' ###
#' ### simTrial(N=c(1000, 1000), aveVE=c(0, 0.4),
#' ### VEmodel="half", vePeriods=c(1, 27, 79), enrollPeriod=78,
#' ### enrollPartial=13, enrollPartialRelRate=0.5, dropoutRate=0.05,
#' ### infecRate=0.06, fuTime=156, visitSchedule=seq(0, 156, by=4),
#' ### missVaccProb=0.05, VEcutoffWeek=26, nTrials=5,
#' ### stage1=78, saveDir="./", randomSeed=300)
#' ###
#' ### THIS CALL NEEDS TO BE REVISED TO A SIMPLE BUT WORKING CODE
#' ### THE INTENT HERE IS ONLY TO ILLUSTRATE HOW OUTPUT FILES CAN BE READ IN
#' ### THE PURPOSE IS NOT TO SHOW ANY ADDITIONAL TRIAL DESIGNS
#' ### monitorTrial(dataFile=
#' ### "simTrial_nPlac=1400_nVacc=1000_1000_aveVE=0.2_0.4_infRate=0.04.RData",
#' ### stage1=78, stage2=156, harmMonitorRange=c(10,100), alphaPerTest=NULL,
#' ### nonEffStartMethod="FKG", nonEffInterval=20, lowerVEnoneff=0,
#' ### upperVEnoneff=0.4, highVE=0.7, stage1VE=0, lowerVEuncPower=0,
#' ### alphaNoneff=0.05, alphaHigh=0.05, alphaStage1=0.05, alphaUncPower=0.05,
#' ### estimand="cuminc", lagTime=26, saveDir="./")
#'
#' @seealso \code{\link{simTrial}}, \code{\link{censTrial}}, \code{\link{rankTrial}}, \code{\link{estHRbound}}, \code{\link{crossBoundCumProb}}, and \code{\link{decisionTimes}}
#'
#' @export
monitorTrial <- function (dataFile,
stage1,
stage2,
## range over which to "spend" the type-I error specified in argument
## 'harmMonitorRange'. It should be a vector of length 2 giving the
## (start, stop) infection-count range over which the type-I error
## will be spent. Please note that this argument does *NOT* dictate
## when harm monitoring will END. The range dictates when it will
## START, and over what range the type-I error will be spread. If
## you plan to use nonEffStartMethod "FKG" or "fixed" without lagged
## Monitoring (i.e. laggedMonitoring=FALSE) then the 'stop' value of
## the range, if provided, will not be used and need not be specified
## or may be replaced with an NA.
harmMonitorRange,
## Total Type I error for potential harm monitoring (per vacc. arm)
harmMonitorAlpha=0.05,
## 'harmMonitorContol' is a list containing other parameters controlling
## the harm monitoring. Currently the only one used in 'maxCnt'. This
## *MUST* be specified and represesents the set of infection counts for
## which harm bounds will be constructed. Harm monitoring ends when
## 'maxCnt' has been exceeded, whether or not you've reached the
## criteria for initiaton of non-efficacy monitoring, so you should:
## (a) choose a nonEff start method that guarantees starting
## at/before 'maxCnt', or
## (b) choose maxCnt to be a value larger than you're non-efficacy
## starting count will ever be (maybe use 2 or 3 times the
## upper range for harm monitoring (harmMonitorRange[2]).
## If maxCnt is left unspecified, the default value of
## 3 x harmMonitorRange[2] will be used for it.
#harmMonitorControl=list( maxCnt=NULL),
## if you have determined the constant alpha value to do the binomial
## 'potential-harm-monitoring' tests at, they you can pass it through
## this argument, in which case argument 'harmMonitorAlpha' need not
## be specified
alphaPerTest=NULL,
## character vector of methods to use to determine when to start non-efficacy
## monitoring. Each method requires different input information, and that
## information should be input via the 'nonEffStartParams' argument - which
## should be a list that the specifies the inputs for the method. Each
## method's input requirements should be specified in the documentation.
##
## Methods:
## FKG - the starting method suggested in Freidlin, Korn and Gray's 2010
## 'Clinical Trials' paper: "A general inefficacy interim
## monitoring rule for randomized clinical trials"
## It boils down to start monitoring at the earlier infection count
## such that an estimated effect <= 0 would cause the trial to
## stop. In our code that translates to a 95% CI (based on the
## asymptotic variance of the log-rank statistic) around an
## estimated VE of 0% would exclude the VE specified by parameter
## 'upperVEnonEff'. The expectation is that this will be harmonized
## with the non-efficacy monitoring, and so will use the arguments
## 'upperVEnonEff' 'alphaNoneff' from argument list. However, should
## you want to use different values, you may do say by passing them
## via the 'nonEffStartParams' argument list. The same naming should
## be used within the list.
## Parameters: upperVEnonEff alphaNoneff
##
## fixed - Starts at the infection count specified by paramter 'N1'.
## If 'N1' is set to 75, then non-efficacy always begins at the
## 75th infection
## Parameter: N1
##
nonEffStartMethod=c("FKG", "fixed"),
## A *list* (not a vector) of parameters needed by the method specified
## in argument 'nonEffStartMethod'. Some methods have defaults in place,
## for those you do not need to use 'nonEffStartParams' unless you want
## values other than the defaults.
nonEffStartParams=NULL,
nonEffIntervalUnit=c("counts","time"),
# a single numeric value (for constant increments) or a numeric vector
# (for variable increments)
nonEffInterval,
# if noneff monitoring is being used (which it always is), then
# nonEffCohorts must provide info on each cohort being monitored.
# the argument is a list of lists. The internal lists should be
# named "cohort1", "cohort2", ..., "cohort<N>".
# Each cohort's list must contain a component named 'estimand'
# specifying either 'cox','cuminc', or 'combined' (for both).
# Optional list components are:
# 'lagTime' - indicates the lag time that controls event inclusion
# for the cohort. If no lag time is desired, then this
# component can be excluded, set to NULL, or set to 0
# 'cohortInd'- Character string naming an indicator variable to be
# used to subset participants to the desired cohort.
# The variable must already exist in the simTrials
# output. This allows inclusion of e.g a "per-protocol"
# variable
#
# Example specifications:
# - 2 cohorts, nominal CI monitoring
# nonEffCohorts <- list(
# times=NULL, timeUnit="counts",
# timingCohort=list(lagTime=NULL, cohortInd=NULL),
# cohort1 = list(
# estimand="cox", lagTime=NULL, nullVE=0.5,
# nominalAlphas=0.05
# ),
# cohort2 = list(
# estimand="cox", lagTime=6, cohortInd="pp1",
# nullVE=0.5, nominalAlphas=0.05)
# )
# )
#
# - 1 cohort, sequential control of type-I error
# nonEffCohorts <- list(
# times=c(50,75,100), timeUnit="counts",
# timingCohort=list(lagTime=NULL, cohortInd=NULL),
# cohort1 = list(
# estimand="cox", lagTime=NULL, cohortInd="pp1",
# nullVE=0.6, nominalAlphas=c(0.0030, 0.0183, 0.0440)
# ),
# )
#
## The default set-up (below) monitors one cohort, with 'cox' estimand
## and lagTime defaulting to NULL
nonEffCohorts = list(
times=NULL, timeUnit="counts",
timingCohort=list(lagTime=NULL, cohortInd=NULL),
cohort1 = list(
estimand="cox", lagTime=NULL, cohortInd=NULL,
nullVE=NULL, nominalAlphas=NULL
)
),
# timeLag = lag to use in converting 'times' counts into time
# (not needed if timeUnit = "time" or if no time lag is used)
#
# -- These args are only needed if using sequential monitoring bounds
# + totalAlpha = total (two-sided) alpha for efficacy testing
# + totalEvents = number of events expected (or planned for) at the
# last analysis time. Needed to convert to 'information time'
# for use of spending function approach
# + spendingFunction = Tells which kind of bound to create
# Specification is either: (1) the text string "OBF", indicating
# use of O'Brien-Fleming bounds, or (2) an alpha spending function.
# The alpha spending function (ASF) 'f' must satisfy the
# following conditions:
# (a) takes a single argument with values in [0,1],
# (b) returns values in [0,1]
# (c) must be strictly increasing
# (d) f(0)=0 and f(1)=1
# We provide a function asf() which allows creation of such
# functions for some common families of bounds.
# asf("Pocock") creates a pocock ASF,
# asf("power", phi=x) creates an ASF from the power family
# with parameter phi set to 'x'.
# asf("HSD", phi=x) creates an ASF from the Hwang-Shih-DeCani
# Family with parameter phi set to 'x'.
#
# More optional components are:
# 'lagTime' - indicates the lag time that controls event inclusion
# for the cohort. If no lag time is desired, then this
# component can be excluded, set to NULL, or set to 0
# 'cohortInd'- Character string naming an indicator variable to be
# used to subset participants to the desired cohort.
# The variable must already exist in the simTrials
# output. This allows inclusion of e.g a "per-protocol"
# variable
# NOT ALL IMPLEMENTED YET
#
# Minimum required:
# times, timeUnit, estimand, lagTime, totalAlpha, nullVE,
# totalEvents
# effCohort = list( times, timeUnit, timeLag, timeCohortInd,
# estimand, lagTime, cohortInd,
# totalAlpha, nullVE, totalEvents,
# spendingFunction),
effCohort = list( times=NULL, timeUnit="counts",
timingCohort=list(lagTime=NULL, cohortInd=NULL),
nullVE = NULL,
estimand="cox", lagTime=NULL, cohortInd=NULL,
nominalAlphas=NULL),
## lowerVEnoneff is not required. Specify only if you want this
## condition as part of your monitoring.
lowerVEnoneff=NULL,
#upperVEnoneff,
highVE,
stage1VE,
lowerVEuncPower=NULL,
#alphaNoneff,
alphaHigh,
alphaStage1,
alphaUncPower=NULL,
saveFile= NULL,
saveDir = NULL,
verbose = TRUE ) {
## selected Arguments
##
## lowerVEnoneff:
## lower confidence bound to be below 'lowerVEnoneff' to meet criterion 1
## for declaring non-efficacy
##
## upperVEnoneff:
## upper confidence bound to be below 'upperVEnoneff' to meet criterion 2
## for declaring non-efficacy
##
## highVE:
## lower confidence bound to be above 'highVE' for declaring high efficacy
##
## stage1VE:
## lower confidence bound to be above 'stage1VE' for advancing into Stage 2
##
## lowerVEuncPower:
## vector of VE value for which the user wishes to determine unconditional
## power to reject H0: VE <= (lowerVEuncPower * 100)%
## at 2-sided alpha level 'alphaUncPower'
##
## alphaNoneff:
## One minus confidence level of 2-sided CI for non-efficacy monitoring
##
## alphaHigh:
## One minus confidence level of 2-sided CI for high efficacy monitoring
##
## alphaStage1:
## One minus confidence level of 2-sided CI for testing whether an arm
## advances into Stage 2
##
## alphaUncPower:
## One minus confidence level of 2-sided CI for unconditional power to
## reject H0: VE <= (lowerVEuncPower * 100)%
nonEffIntervalUnit <- match.arg(nonEffIntervalUnit)
#nonEffTimeUnit <- match.arg(nonEffTimeUnit)
nonEffStartMethod <- match.arg(nonEffStartMethod)
if (!is.null(nonEffCohorts$times) ) {
## disable start method, since not needed
nonEffStartMethod <- ""
}
if ( is.list(dataFile) ) {
trialObj <- dataFile
} else {
if ( !is.null(saveDir) ){
## load in RData object (a list named 'trialObj' )
load(file.path(saveDir, dataFile))
} else {
load(dataFile)
}
}
## check contents of 'trialData'
if ( !all( c("entry","exit","event","trt") %in% names(trialObj$trialData[[1]]) ) )
stop("Trial data must contain columns: entry, exit, event, trt\n")
nTrtArms <- as.integer( trialObj$nArms - 1 )
nTrials <- length(trialObj[["trialData"]])
## set altVE based on lowerVEuncPower, etc. if those were specified instead
alphaAltVE <- altVE <- NULL
if ( is.null(altVE) && !is.null(lowerVEuncPower)) {
altVE <- lowerVEuncPower
}
if (is.null(alphaAltVE) && !is.null(altVE)) {
alphaAltVE <- ifelse(!is.null(alphaUncPower), alphaUncPower, alphaStage1)
}
## ----------------------------------------------------------------------
## The following code is needed to determine infection count at which to
## begin the nonEff monitoring. Also needed for harm-monitoring
## derive the prob. of assignment to each vaccine arm (relative to placebo
## only, not relative to other vaccine arms)
null.p <- numeric(length=nTrtArms)
for (i in 1:nTrtArms) {
arms.i <- c(1, i+1)
assn.probs <- trialObj$trtAssgnProbs[ arms.i ]
null.p[i] <- assn.probs[2] / sum(assn.probs)
}
## if we have multiple vaccine groups and they are of different sizes,
## then we must stop. The code is not set up to handle that yet.
if (nTrtArms > 1 && ( max(null.p)-min(null.p) > sqrt(.Machine$double.eps) ) )
{
stop("The code is not currently set up to handle harm monitoring when\n",
"there are multiple treatment arms and their sizes differ. \n\n")
} else {
## all values of null.p are equal, just choose the first
null.p <- null.p[1]
}
## This will also determine when nonEfficacy monitoring starts
## (and therefore how long harm monitoring will last).
## Two methods are found here and two later on. These two do not
## require any specific information about the trial to determine,
## whereas the later methods do.
## Method implemented here:
## Start nonEff monitoring at the smallest infection count for which we
## would be able to reject the null hypothesis that VE = upperVEnoneff
## if our point estimate of VE = 0. This computation requires us to
## assume that our test statistic is asymptotically equivalent to the
## log-rank statistic, which (under the null of no diff - not the setting
## we're in, right?) has a distribution that is
## Normal( 0, sqrt(1/(p*(1-p)*D))), where
## D = number of endpoints between the placebo and trt arm being monitored,
## p = prob. of assignment to the vaccine arm vs. the placebo arm (all
## other arms ignored for computing 'p')
##
## if 'timeUnit' is not specified but 'timingCohort' is, assume timeUnit is 'counts'
if ( is.null(nonEffCohorts$timeUnit) ) {
if ( !is.null(nonEffCohorts$timingCohort) ) {
nonEffCohorts$timeUnit <- "counts"
} else {
stop("You must specify 'timeUnit' within 'nonEffCohorts' as 'time' or specify",
" 'timingCohort' to provide info needed to translate counts into times\n\n")
}
}
## Create indicators of whether either a lagTime or cohortInd has been
## specified for timing of non-efficacy analyses
if ( nonEffCohorts$timeUnit != "time" ) {
if ( !is.null(nonEffCohorts$timingCohort$lagTime) &&
nonEffCohorts$timingCohort$lagTime > 0 ) {
laggedMonitoring <- TRUE
} else {
laggedMonitoring <- FALSE
}
neTimeCohInd <- ifelse( !is.null(nonEffCohorts$timingCohort$cohortInd),
TRUE, FALSE )
## if either indicator is set (is "ne" timing diff. from Harm timing
neTimingDiff <- ifelse(laggedMonitoring || neTimeCohInd, TRUE, FALSE )
}
if ( !is.null(nonEffCohorts$times) && nonEffCohorts$timeUnit == "counts") {
nonEffCnts <- nonEffCohorts$times
if ( !neTimingDiff ) {
N1 <- nonEffCnts[1]
N1.ne <- N1
} else {
N1.ne <- nonEffCnts[1]
}
}
if ( nonEffStartMethod == "FKG") {
if (!is.null(nonEffStartParams)) {
## values specified directly by user (maybe - we don't
## know what's actually in 'nonEffStartMethod' yet
upperVE <- nonEffStartParams$upperVEnoneff
alphaNE <- nonEffStartParams$alphaNoneff
if ( is.null(upperVE) && is.null(alphaNE) )
warning("The argument 'nonEffStartParams' has been specified but",
"it does not contain\n", "parameters used by method FKG\n\n",
immediate.=TRUE)
} else {
upperVE <- NULL
alphaNE <- NULL
}
## Fill in NULL values for either parameter using the monitorTrial
## arguments upperVEnoneff and alphaNoneff
if ( is.null(upperVE) )
upperVE <- nonEffCohorts$cohort1$nullVE
if ( is.null(alphaNE) )
alphaNE <- nonEffCohorts$cohort1$nominalAlphas
if ( is.null(upperVE) || is.null(alphaNE) )
stop( "One or both of the arguments needed by nonEffStartMethod 'FKG'",
"are missing\n.", "They were not specified via argument ",
"'nonEffStartParams' and their values also could not be obtained\n",
"from the arguments 'nonEffCohorts$cohort1$nullVE' and ",
" 'nonEffCohorts$cohort1$nominalAlphas'\n",
"Exiting...\n" )
N1.ne <- ceiling(qnorm(1 - alphaNE/2)^2 /( null.p*(1-null.p)*log(1-upperVE)^2 ))
if ( !neTimingDiff ) {
N1 <- N1.ne
}
} else if ( nonEffStartMethod == "fixed" ) {
if (is.null( nonEffStartParams$N1 ) )
stop("Argument 'nonEffStartMethod' was specified as 'fixed'. ",
"This method requires\n", "you to provide an argument named",
"'N1' via the list argument 'nonEffStartParams'\n.",
"e.g. nonEffStartParams <- list(N1=60). Please fix.\n" )
N1.ne <- nonEffStartParams$N1
if (!neTimingDiff ) {
N1 <- N1.ne
}
}
## N1 should only exist if laggedMonitoring is FALSE, while N1.ne exists
## otherwise (N1.ne also exists if N1 does).
## Harm monitoring will end at N1 infections, so when "N1" exists, set the upper
## limit of 'harmMonitorRange' to N1 so boundaries will be defined based on it
if ( exists("N1") ) {
harmMonitorRange <- c(harmMonitorRange[1], N1)
}
## -------------------------- get harm bounds ---------------------------
## do a check on 'harmMonitorRange' to be sure it's correctly specified
if ( length(harmMonitorRange) != 2 || is.na(harmMonitorRange[2]) ||
is.infinite(harmMonitorRange[2]) ) {
if (neTimingDiff) {
stop("When the non-efficacy timing cohort differs from the harm monitoring\n",
"cohort, then the harmMonitoringRange must be fully specified - i.e. it\n",
"must be length 2 vector with both values being integers.\n\n")
} else {
stop("When nonEffStartMethod is NOT one of 'FKG' or 'fixed', then ",
"harmMonitoringRange must be fully specified.\n", "It should be a ",
"length 2 vector with both values being integers.\n\n")
}
}
# calculate stopping boundaries for harm
## NOTE: harmMonitorRange dictates that range over which the type-I error is
## spent, it does not specify the time at which harm monitoring stops. That is
## dictated by when non-eff monitoring begins
if ( is.null(alphaPerTest) ) {
## choose the value of alphaPerTest to spend the type I error over harmMonitorRange
alphaPerTest <- getAlphaPerTest(harmMonitorRange, null.p,
totalAlpha = harmMonitorAlpha)
}
## For generating the bounds, we need both 'N' and the upper limit of
## 'harmMonitorRange' to be an upper limit on the number of tests to do
## (i.e. an upper limit on when the first non-efficacy will be done)
if (exists("N1")) {
maxCnt = N1
} else if ( !exists("maxCnt") ) {
maxCnt = 2 * harmMonitorRange[2]
}
## get harm bounds
harmBounds <- getHarmBound(
N = maxCnt,
per.test = alphaPerTest,
harmBoundRange = c(harmMonitorRange[1], maxCnt),
null.p = null.p,
dataDir = saveDir,
verbose = verbose)
## -------------------------- end get harm bounds ---------------------------
## creates a list of length 'nTrials' each element of which is a list of
## length 'nTrtArms'
out <- rep( list(vector("list",nTrtArms)), nTrials )
## censor all trials to stage 1 (fastest to do all at once). Still need to
## keep access to uncensored data, for highEff monitoring
if ( !missing(stage1) ) {
stg1dat <- censorTrial(trialObj$trialData, times=stage1,
timeScale="follow-up")
}
## Make note of whether 'N1' exists (it will for certain types of non-eff. monitoring
## timing). If it does, then N1 will remain constant throughout the trial.
## Otherwise we'll need to reset it at the beginning of each new trial/arm
constant.N1 <- ifelse( exists("N1"), TRUE, FALSE )
## Start looping over trials
for (i in 1:nTrials ) {
## extract data for the i-th trial
if ( !missing(stage1) ) {
datI <- stg1dat[[i]]
datIall <- trialObj[["trialData"]][[ i ]]
} else {
datI <- trialObj[["trialData"]][[ i ]]
}
## create separate set of only infections ('events'), then
## order the events by trial time at which they are observed
eventDF <- datI[datI$event == 1, ]
eventDF <- eventDF[order(eventDF$exit), ]
#cat("Trial :", i, " \n")
## Now start comparisons - each active arm vs. placebo arm, one at a time
for (j in 1:nTrtArms) {
## subset *all data* down to the two arms being compared
datI.j <- datI[datI$trt %in% c(0,j), ]
## convert 'trt' to indicator variable before passing it to
## 'applyStopRules' (i.e. convert the non-zero values to 1)
datI.j$trt <- as.integer(datI.j$trt > 0 )
## same for the 'all' data.frame
if ( exists("datIall") ){
datIall.j <- datIall[datIall$trt %in% c(0,j), ]
datIall.j$trt <- as.integer(datIall.j$trt > 0 )
}
## 1. do harm monitoring for j-th treatment arm vs. placebo
## if we don't have a constant value of N1 (across trials), then
## remove the existing value of N1
if ( !constant.N1 && exists("N1") )
rm("N1")
## subset events in relevant arms: j-th active trt and placebo(trt=0)
E.j <- eventDF[eventDF$trt %in% c(0,j), ]
nInfec <- nrow(E.j) # counts infections through 'stage1'
E.j$nInf <- 1:nInfec
## set N1.ne for the cases for which N1 was just defined (immed. above)
if ( !exists("N1.ne") && exists("N1") )
N1.ne <- N1
## ---------------- Process nonEff timing args -------------------
## We need to process the timing information for non-efficacy monitoring now,
## because it's used to determine when harm monitoring will end
## The code is LONG :(
## determine end of stage 1 time, needed if any specified monitoring
## times/counts are not attained in stage 1.
## [ endStg1 == calendar time at which stage 1 follow-up ends ]
endStg1 <- max( datI.j$exit )
if ( neTimingDiff ) {
if (neTimeCohInd) {
timeCoh <- ( datI.j[[nonEffCohorts$timingCohort$cohortInd]]==1 )
nedatI.j <- datI.j[ timeCoh, ]
} else {
nedatI.j <- datI.j
}
if (laggedMonitoring) {
nedatI.j <- censorTrial( nedatI.j, times=nonEffCohorts$timingCohort$lagTime,
timeScale="follow-up", type="left")
}
nInfec.ne <- sum( nedatI.j$event == 1 )
}
## create new variable, 'nInfec_nonEff' that will contain the total number of
## infections under the chosen non-efficacy monitoring cohort
nInfec_nonEff <- ifelse(neTimingDiff, nInfec.ne, nInfec)
## Determine how many non-eff monitoring cohorts were defined
n.neCohorts <- sum( substr(names(nonEffCohorts),1,6) == "cohort" )
## initialize 'alphaNoneff.ij', a temp list containing the nominal alpha
## values for each non-efficacy cohort. This list we be modified, as req.
## based on the current simulated trial, which is why we reinitialize it
## for trial
alphaNoneff.ij <- lapply(1:n.neCohorts, function(idx, neC) {
neC[[paste0("cohort",idx)]]$nominalAlphas
}, neC=nonEffCohorts)
## determine times of non-efficacy testing (times, not counts), then monitor
if (is.null(nonEffCohorts$times)) {
if (nonEffIntervalUnit == "counts") {
## makes sequence of counts at which analyses will be done
if (length(nonEffInterval)==1){
# then a constant increment in the endpoint count is assumed
nonEffCnts <- seq(from = N1.ne, to = nInfec_nonEff, by = nonEffInterval)
} else {
# then 'nonEffInterval' is a vector of potentially variable increments
# specifying endpoint counts for subsequent interim looks following the
# initial look
nonEffCnts <- cumsum(c(N1.ne, nonEffInterval))
}
## Convert counts into times.
if ( neTimingDiff ) {
## subsetted and lagged (if required) dataset 'nedatI.j' was already created earlier
nonEffTimes.ij <- getInfectionTimes(nedatI.j, cnts=nonEffCnts, lagTime=NULL )
} else {
nonEffTimes.ij <- getInfectionTimes(datI.j, cnts=nonEffCnts)
}
} else {
## User specified a time-sequence for monitoring. Figure out when the first
## time and last times will be, then create the sequence
firstLastnonEffTimes <-
getInfectionTimes(datI.j, cnts=c(N1.ne, nInfec_nonEff),
lagTime=nonEffTimeLag )
if (length(nonEffInterval)==1){
nonEffTimes.ij <-
seq(from = firstLastnonEffTimes[1],
to = firstLastnonEffTimes[2],
by = nonEffInterval)
} else {
# then 'nonEffInterval' is a vector of potentially variable time increments
nonEffTimes.ij <- cumsum( c(firstLastnonEffTimes[1], nonEffInterval) )
## truncate so that nonEffTimes don't go beyond the trial duration
nonEffTimes.ij <- nonEffCohorts$times[ nonEffCohorts$times <= firstLastnonEffTimes[2] ]
}
}
} else {
## if nonEffTimes was specified
if (nonEffCohorts$timeUnit == "time") {
nonEffTimes.ij <- sort(nonEffCohorts$times)
## check if any specified times are beyond the extent of stage 1
if ( any(nonEffTimes.ij > endStg1) ) {
wTooLrg <- which( nonEffTimes.ij > endStg1 )
## if just one is too large, then set that last time to 'endStg1',
## no need to modify 'alpha's
if ( length(wTooLrg) == 1 ) {
nonEffTimes.ij[wTooLrg] <- endStg1
} else if (length(wTooLrg)>1) {
## If 'nominalAlphas' is length 1 for all cohorts, then we don't need
## modify nominal alphas. We set the time of the first missed test to
## the end of stage 1 fu time and drop the remaining tests
if ( all( sapply(alphaNoneff.ij, length) == 1 )) {
nonEffTimes.ij[ wTooLrg[1] ] <- endStg1
nonEffTimes.ij <- nonEffTimes.ij[ 1:wTooLrg[1] ]
} else {
## If not all are length 1, Issue warning but continue anyway
cat("Warning: More than one of the specified non-efficacy test times was ",
"not reached.\n", "The code cannot adjust the nominal alphas properly",
" for this.\n", "A final test will be inserted at the end of stage1 ",
"follow-up and will be carried\n", "out using the nominal alpha ",
"level associated with the final scheduled test\n")
## replace the first missed test with a test at the end of stage 1
## using the nominalAlpha level corresponding to final analysis
nonEffTimes.ij[ wTooLrg[1] ] <- endStg1
nonEffTimes.ij <- nonEffTimes.ij[ 1:wTooLrg[1] ]
## modify alpha for any cohorts that have more than one specified
for (idx in 1:n.neCohorts) {
if ( (L <- length(alphaNoneff.ij[[idx]]) ) > 1 )
alphaNoneff.ij[[idx]] <- alphaNoneff.ij[[idx]][ c(1:(wTooLrg[1]-1), L) ]
}
}
}
}
} else {
if ( neTimingDiff ) {
## subsetted and lagged (if required) dataset 'nedatI.j' was already created earlier
nonEffTimes.ij <- getInfectionTimes(nedatI.j, cnts=nonEffCnts)
} else {
nonEffTimes.ij <- getInfectionTimes(datI.j, cnts=nonEffCnts)
}
## check if we "lost" any tests, due to the counts not being reached
if ( length(nonEffTimes.ij) < length(nonEffCnts) ) {
nLost <- length(nonEffCnts) - length(nonEffTimes.ij)
## if we only "lost" 1 test, or if only one alpha specifed for each noneff
## cohort, then we don't need to modify the alphaNoneff.ij values
if (nLost==1 || all( sapply( alphaNoneff.ij, length) == 1) ) {
nonEffTimes.ij <- c(nonEffTimes.ij, endStg1)
} else {
## if we "lost" >1 test, or if have > 1 alpha specifed for any non-eff
## cohort, then we *do* need to modify some alphaNoneff.ij values
## Issue warning but continue anyway
cat("Warning: More than one of the specified non-efficacy test times was",
" not reached.\n", "The code cannot adjust the nominal alphas properly",
" for this.\n", "A final test will be inserted at the end of stage1 ",
"follow-up and will be carried\n", "out using the nominal alpha level",
" associated with the final scheduled test\n\n")
## add a test at the end stage 1 follow-up. Perform the test using
## using the 'alphaLevelNoneff' value corresponding to final analysis
nonEffTimes.ij <- c(nonEffTimes.ij, endStg1)
## modify alpha for any cohorts that have more than one specified
for (idx in 1:n.neCohorts) {
if ( (L <- length(alphaNoneff.ij[[idx]]) ) > 1 ){
alphaNoneff.ij[[idx]] <-
alphaNoneff.ij[[idx]][ c(1:(length(nonEffTimes.ij)-1), L) ]
}
}
}
}
}
}
## ---------------- End Processing nonEff timing args -------------------
## ----- Process efficacy analysis inputs here so we know what test --------
## ----- to perform if we stop early for harm or non-eff --------
## Now get 'N1' for the cases where we don't have it, which should only
## be when laggedMonitoring was specified
if ( !exists("N1") && neTimingDiff ) {
## determine the total infection count (based on full set of data) that corresponds
## to the first non-eff test time
N1 <- sum( datI.j$event==1 & datI.j$exit <= nonEffTimes.ij[1] )
}
## Ensure we have harm bounds up to 'N1', if not, regenerate them
if (N1 > max(harmBounds$N)) {
maxCnt <- round(1.5 * N1 )
## get harm bounds
harmBounds <- getHarmBound(
N = maxCnt,
per.test = alphaPerTest,
harmBoundRange = c(harmMonitorRange[1], maxCnt),
null.p = null.p,
dataDir = saveDir,
verbose = verbose)
}
## subset to restrict harm monitoring to the first 'N1' infections
harmBounds.j <- harmBounds[harmBounds$N <= N1, c("N","V","P")]
## doEvaluate harm monitoring bounds. Just based on infected ppts, so
## only uses 'E.j'
harmRes <- do_harm_monitoring(E.j, bounds = harmBounds.j )
## if there is indication of "harm" then need to create store output
## object then move to next trial
if ( harmRes$isHarm ) {
if ( !is.null( effCohort$cohortInd ) ) {
datI.j <- datI.j[ datI.j[[ effCohort$cohortInd ]] == 1, ]
}
fst <- finalStage1Test(
datI.j,
analysisType = "stopTime",
lowerVE = altVE,
alphaLevel = alphaAltVE,
estimand = effCohort$estimand,
lagTime = effCohort$lagTime,
time = harmRes$stopTime,
randFraction = null.p )
fest <- fst$stage1summ
VEcir <- !is.null(fest$VE_CIR)
VE_CI <- if (VEcir) {
c(fest$VE_CIR_loCI, fest$VE_CIR_upCI)
} else {
c(fest$VE_Cox_loCI, fest$VE_Cox_upCI)
}
out[[i]][[j]] <-
list(
boundHit = "Harm",
stopTime = harmRes$stopTime,
stopInfectCnt = harmRes$stopInfectCnt,
stopInfectSplit = harmRes$stopInfectSplit,
stopVE = ifelse(VEcir, fest$VE_CIR, fest$VE_Cox),
stopVE_CI= VE_CI,
VE_estimand = effCohort$estimand,
stage1complete = FALSE,
stage2complete = FALSE,
altDetected = fst$altDetected )
next
}
## If harm bounds not hit, move to nonEff monitoring
## 2. begin non-eff monitoring prep
## --------------------------------
## I don't think I need to skip this situation. The later code that deals with
## test times that are not reached should be used to handle it
#if (N1.ne > nInfec_nonEff ) {
## deal with annoying case where the trial didn't have enough infec.s
## to reach the start point of the non-eff. interim analyses
## For now we do nothing, except skip all the interim analysis code
## and go directly to the finalStage1 test
# nonEffTimes.ij <- NA
## create an empty object futRes needed later to make code work
# futRes <- list()
# effRes <- list()
## Set 'N1' to NA ?? Need to think about it
## N1 <- NA
#} else { # This 'else' goes on FOREVER!
## ------- run non-eff ---------
## We need to loop over the defined cohorts, then concatenate results across them
## The details on the non-eff cohorts exist within sublists named 'cohort1',
## 'cohort2'...
#n.neCohorts <- sum( substr(names(nonEffCohorts),1,6) == "cohort" )
cohortFutRes <- vector( "list", n.neCohorts )
for(coh.idx in 1:n.neCohorts) {
## name of component
coh.i <- paste0("cohort", coh.idx)
ind <- nonEffCohorts[[ coh.i ]]$cohortInd
## if cohortInd is specified, use it to subset the data
if (!is.null( ind ) ) {
cohI.j <- datI.j[ datI.j[[ind]] == 1, ]
} else {
cohI.j <- datI.j
}
cohortFutRes[[ coh.idx ]] <-
applyStopRules (
cohI.j,
testTimes = nonEffTimes.ij,
boundType = "nonEff",
boundLabel = "NonEffInterim",
lowerVE = lowerVEnoneff,
upperVE = nonEffCohorts[[ coh.i ]]$nullVE,
alphaLevel = alphaNoneff.ij[[coh.idx]],
lagTime = nonEffCohorts[[ coh.i ]]$lagTime,
estimand = nonEffCohorts[[ coh.i ]]$estimand,
returnAll = ifelse(n.neCohorts>1, TRUE, FALSE),
randFraction = null.p )
}
if (n.neCohorts > 1) {
## extract the 'summObj' components
summObj <- lapply(cohortFutRes, function(y) y$summObj)
names(summObj) <- paste0("cohort", 1:n.neCohorts)
## if all cohorts report a bound hit, then we'll need to dig in deeper
## to see if they all hit at the same timepoint ever.
## If they didn't all hit, then we're done here!
if ( all( sapply(cohortFutRes, function(y) y$boundWasHit) ) ) {
# extract the stopTimes from each cohort
stopIndxList <- lapply( cohortFutRes, function(x) x$stopIndx )
## find common values in all components
commonStopIndx <- unlist( Reduce(intersect, stopIndxList) )
if ( length(commonStopIndx) > 0 ) {
stopIndx <- min( commonStopIndx )
## time of test corresponding to stopIndx
stopTime <- nonEffTimes.ij[ stopIndx ]
## subset each set of results in summObj to rows 1:stopIndx
summObj <- lapply(summObj, function(y) y[1:stopIndx,] )
boundWasHit <- TRUE
boundHit <- "NonEffInterim"
} else {
boundWasHit <- FALSE
boundHit <- NA
stopTime <- NA
}
} else {
boundWasHit <- FALSE
boundHit <- NA
stopTime <- NA
}
futRes <- list( boundWasHit = boundWasHit,
boundHit = boundHit,
boundType = "NonEffInterim",
stopTime = stopTime,
summObj = summObj)
} else {
## just one cohort, so set futRes to it
futRes <- cohortFutRes[[1]][ c("boundWasHit", "boundHit",
"boundType", "stopTime", "summObj") ]
}
## if noneff hit, then store results and go to next arm
## (note: output will be a list because 'futRes' is a list
if ( futRes$boundWasHit ) {
if ( !is.null( effCohort$cohortInd ) ) {
datI.j <- datI.j[ datI.j[[ effCohort$cohortInd ]] == 1, ]
}
fst <- finalStage1Test(
datI.j,
analysisType = "stopTime",
lowerVE = altVE,
alphaLevel = alphaAltVE,
estimand = effCohort$estimand,
lagTime = effCohort$lagTime,
time = futRes$stopTime,
randFraction = null.p )
fest <- fst$stage1summ
VEcir <- !is.null(fest$VE_CIR)
VE_CI <- if (VEcir) {
c(fest$VE_CIR_loCI, fest$VE_CIR_upCI)
} else {
c(fest$VE_Cox_loCI, fest$VE_Cox_upCI)
}
out[[i]][[j]] <-
c( futRes,
list( stopVE = ifelse(VEcir, fest$VE_CIR, fest$VE_Cox),
stopVE_CI= VE_CI,
VE_estimand = effCohort$estimand,
firstNonEffCnt = N1,
stage1complete = FALSE,
stage2complete = FALSE,
altDetected = fst$altDetected ) )
next
}
## Insert efficacy analyses here. We'll need to rework later to
## allow it to be executed earlier, so parts of it can be included
## if some test are done prior to stopping for harm or non-eff
## There's no chance that we would have stopped for efficacy
## but instead ended up stopping for harm or non-eff, the two
## results are two disparate
#effCohort = list( times=c(52,77,103), timeUnit="counts", timeLag=2,
# nullVE = 0.2, totalAlpha=0.05, totalEvents=103,
# estimand="cox", lagTime=2,
# nominalAlphas=c(0.0030, 0.0183, 0.0440) ),
## clear out variable, if it exists
if ( exists("nominalAlphas.ij") )
rm( nominalAlphas.ij )
## if
## Process to derive 'effTimes.ij', etc.
if (effCohort$timeUnit == "time") {
effTimes.ij <- sort(effCohort$times)
## check if any specified times are beyond the extent of stage 1
if ( any(effTimes.ij > endStg1) ) {
wTooLrg <- which( effTimes.ij > endStg1 )
## if just one is too large, then set that last time to 'endStg1'
if ( length(wTooLrg) == 1 ) {
effTimes.ij[wTooLrg] <- endStg1
} else if ( length(effCohort$nominalAlphas) == 1 ) {
effTimes.ij[ wTooLrg[1] ] <- endStg1
effTimes.ij <- effTimes.ij[ 1:wTooLrg[1] ]
} else {
## Issue warning but continue anyway
cat("Warning: More than one of the specified efficacy test times was ",
"not reached.\n", "The code cannot adjust the nominalAlphas properly",
" for this.\n", "A final test will be inserted at the end of stage1 ",
"follow-up and will be carried\n", "out using the nominal alpha ",
"level associated with the final scheduled test\n")
## replace the first missed test with a test at the end of stage 1.
## test using 'alphaLevelNoneff' corresponding to final analysis
effTimes.ij[ wTooLrg[1] ] <- endStg1
effTimes.ij <- effTimes.ij[ 1:wTooLrg[1] ]
nominalAlphas.ij <-
effCohort$nominalAlphas[ c( 1:(wTooLrg[1]-1),
length(effCohort$nominalAlphas) ) ]
}
}
} else {
## extract info on timingCohort
if ( is.null( effCohort$timingCohort) )
stop("Component 'timingCohort' must be specified for the efficacy analysis\n.",
"It is needed to translate counts to calendar time for analyses\n")
if ( is.null(effCohort$timingCohort$lagTime) )
stop("Component 'lagTime' is a required component of 'timingCohort' in the\n",
"efficacy testing specification. Set it's value to zero if not needed\n")
if (!is.null( effCohort$timingCohort$cohortInd ) ) {
## subset data using 'cohortInd' and then determine the times
timingI.j <- datI.j[ datI.j[[ effCohort$timingCohort$cohortInd ]] == 1, ]
effTimes.ij <- getInfectionTimes(
timingI.j, cnts=sort(effCohort$times),
lagTime= effCohort$timingCohort$lagTime )
} else {
## cohortInd not specified, just use lagTime
effTimes.ij <- getInfectionTimes(
datI.j, cnts=sort(effCohort$times),
lagTime= effCohort$timingCohort$lagTime )
}
## check if we "lost" any tests, due to the counts not being reached
if ( length(effTimes.ij) < length(effCohort$times) ) {
nLost <- length(effCohort$times) - length(effTimes.ij)
## if we only lost 1, or if we're using a fixed nominal alpha level
## then add a test at the end of stage 1 follow-up
if (nLost==1 || length(effCohort$nominalAlphas)==1 ) {
effTimes.ij <- c(effTimes.ij, endStg1)
} else {
## Issue warning but continue anyway
cat("Warning: More than one of the specified efficacy test times was",
" not reached.\n", "The code cannot adjust the nominalAlphas properly",
" for this.\n", "A final test will be inserted at the end of stage1 ",
"follow-up and will be carried\n", "out using the nominal alpha level",
" associated with the final scheduled test\n\n")
## add a test at the end stage 1 follow-up. Perform the test using
## using the 'alphaLevelNoneff' value corresponding to final analysis
effTimes.ij <- c(effTimes.ij, endStg1)
nominalAlphas.ij <-
effCohort$nominalAlphas[ c( 1:(length(effTimes.ij)-1),
length(effCohort$nominalAlphas) ) ]
}
}
}
## If we've not created a customized set of values 'nominalAlphas.ij' in the
## code above, create it now
if ( !exists("nominalAlphas.ij") )
nominalAlphas.ij <- effCohort$nominalAlphas
if ( !is.null( effCohort$cohortInd ) ) {
effdatI.j <- datI.j[ datI.j[[ effCohort$cohortInd ]] == 1, ]
} else {
effdatI.j <- datI.j
}
effRes <-
applyStopRules(
effdatI.j,
testTimes = effTimes.ij,
boundType = "eff",
boundLabel = "Eff",
lowerVE = effCohort$nullVE,
alphaLevel = nominalAlphas.ij,
estimand= effCohort$estimand,
lagTime = effCohort$lagTime,
randFraction = null.p )
## 3. if nonEff NOT hit, then move on to high-eff monitoring
## We first do *only* at the stage1 times (i.e. same times as for
## nonEff monitoring ), but using all data
highEffRes <- applyStopRules(
datIall.j,
testTimes = nonEffTimes.ij,
boundType = "highEff",
boundLabel = "HighEff", ## assumed by later function
lowerVE = highVE,
alphaLevel = alphaHigh,
estimand = effCohort$estimand,
randFraction = null.p )
if ( highEffRes$boundWasHit ) {
fst <- finalStage1Test(
effdatI.j,
analysisType = "stopTime",
lowerVE = altVE,
lagTime = effCohort$lagTime,
time = highEffRes$stopTime,
randFraction = null.p )
fest <- fst$stage1summ
VEcir <- !is.null(fest$VE_CIR)
VE_CI <- if (VEcir) {
c(fest$VE_CIR_loCI, fest$VE_CIR_upCI)
} else {
c(fest$VE_Cox_loCI, fest$VE_Cox_upCI)
}
## store output, then determine value for 'altDetected'
## (note: output will be a list because 'highEffRes' is a list
out[[i]][[j]] <-
c( highEffRes,
list( summNonEff = futRes$summObj,
stopVE = ifelse(VEcir, fest$VE_CIR, fest$VE_Cox),
stopVE_CI= VE_CI,
VE_estimand = effCohort$estimand,
firstNonEffCnt = N1,
stage1complete = FALSE,
stage2complete = FALSE,
altDetected = fst$altDetected ) )
next
}
#} ## end of (N1.ne > nInfec_nonEff) if/else clause
##-----------------------------------------------------------------------
## --- skip this for now ---
if (FALSE) {
## Reached the end of stage 1 w/o hitting any boundaries (yee-ha!)
## celebrate by doing 'end-of-stage-1' dance
fst <- finalStage1Test(
datI.j,
analysisType = "final",
stage1VE = stage1VE,
lowerVE = altVE,
alphaLevel = alphaStage1,
estimand = effCohort$estimand,
lagTime = effCohort$lagTime,
boundLabels=c("Eff", "NonEffFinal"),
randFraction = null.p )
fest <- fst$stage1summ
VEcir <- !is.null(fest$VE_CIR)
VE_CI <- if (VEcir) {
c(fest$VE_CIR_loCI, fest$VE_CIR_upCI)
} else {
c(fest$VE_Cox_loCI, fest$VE_Cox_upCI)
}
## not storing into out[[i]][[j]] yet, we'll add on if not final-non-eff
finalStg1List <-
c( fst,
list( stopVE = ifelse(VEcir, fest$VE_CIR, fest$VE_Cox),
stopVE_CI= VE_CI,
VE_estimand = effCohort$estimand,
firstNonEffCnt = N1,
summNonEff = futRes$summObj,
stage1complete = TRUE,
stage2complete = FALSE ) )
## if we hit the bound at stage 1 analysis, then we're done with this cohort
if (fst$boundHit == "NonEffFinal") {
out[[i]][[j]] <- finalStg1List
next
}
}
## --- End of skipping for now ---
## -- Temp replacement for skipped code --
finalStg1List <-
c( effRes[ c("boundWasHit", "boundHit", "boundType",
"stopIndx", "stopTime") ] ,
altDetected = effRes$boundWasHit,
list(
# stopVE = ifelse(VEcir, fest$VE_CIR, fest$VE_Cox),
# stopVE_CI= VE_CI,
# VE_estimand = effEst,
firstNonEffCnt = N1,
summNonEff = futRes$summObj,
summEff = effRes$summObj,
stage1complete = TRUE,
stage2complete = FALSE ) )
## -- End temp replacement for skipped code --
## 5. If no stage1 highEff, and we "passed" the end-of-stg1 analysis, then do
## the stage2 highEff analysis
## (a) timing of it is half way between end of stage 1 and end of trial
## compute as end-of-stg1 (max(datI.j$exit)) plus half the remaining time
## (stage2- stage1)/2 until the last person from stage1 completes stg2
stg2highEffTime <- ( max(datI.j$exit) + (stage2 - stage1)/2 )
stg2highEff <-
applyStopRules(
datIall.j,
testTimes = stg2highEffTime,
boundType = "highEff",
boundLabel = "HighEff", ## assumed by later function
lowerVE = highVE,
alphaLevel = alphaHigh,
estimand=effCohort$estimand,
randFraction = null.p )
## Note: 'altDetected' result is included in 'finalStg1List', which is why
## you will not see it added to any of the following output objects
if ( stg2highEff$boundWasHit ) {
## alter some values in present in 'finalStg1List'
finalStg1List$boundHit <- stg2highEff$boundHit
finalStg1List$boundType <- stg2highEff$boundType
finalStg1List$stopTime <- stg2highEff$stopTime
## remove components from stg2Eff already present in 'finalStg1List'
dropComps <- c("boundWasHit", "boundHit", "boundType", "stopTime")
stg2HE <- stg2highEff[ !names(stg2highEff) %in% dropComps ]
## Note: 'altDetected' result is included in 'finalStg1List'
out[[i]][[j]] <- c( finalStg1List,
#list( summNonEff = futRes$summObj ),
stg2HE )
} else {
## Stage 2 completed without hitting high-efficacy
finalStg1List$stage2complete <- TRUE
finalStg1List$stopTime <- max( datIall.j$exit)
out[[i]][[j]] <- c( finalStg1List,
list(
#summNonEff = futRes$summObj,
stg2highEff = stg2highEff$summObj ) )
}
}
}
##-----------------------------------------------------------------------
# bound stage1complete altDetected(stage1 tests - 2 of them)
# "NonEffInterim" FALSE FALSE
# "HighEff" FALSE TRUE/FALSE [could be either]
# "NonEffFinal" TRUE FALSE
# "Eff" TRUE TRUE/FALSE [could be either]
# "HighEff" TRUE (Same as for "Eff")
##-----------------------------------------------------------------------
## I left this in here, but not very useful - will excise later probably
if (verbose){
for (i in 1:nTrtArms) {
cat("Probabilities of reaching each possible conclusion:\n")
print( round( table(sapply( out, function(x) x[[i]]$boundHit ),
useNA="ifany")/nTrials, 4))
if (!is.null(altVE)){
altDetectedMatrix <-
do.call("rbind",
lapply(out, function(x) x[[i]]$altDetected ))
# this occurs if each trial stops for potential harm, and thus the value of
# 'altDetected' is NULL for each trial
if ( is.null( altDetectedMatrix ) ) {
designPower <- 0
} else {
designPower <- colSums(altDetectedMatrix, na.rm=TRUE)/nTrials
}
cat("\nUnconditional power to reject the specified null hypotheses =",
format(designPower, digits=3, nsmall=3), "\n\n")
}
}
}
## save monitoring output
if ( !is.null(saveDir) ) {
if ( is.null(saveFile) ) {
if ( is.list(dataFile) )
warning(
"The output of 'monitorTrial' will not be saved to a file.\n",
"You have not specified the argument 'saveFile' and a default\n",
"filename cannot be constructed when argument 'dataFile' is ",
"a list.\n\n", immediate.=TRUE)
## used in specifying saveFile name
neEst <- ifelse( !is.null(nonEffCohorts[[1]]$estimand),
paste0("_", nonEffCohorts[[1]]$estimand), "")
saveFile <- paste0("monitorTrial",
substr(dataFile, 9, nchar(dataFile)-6),
neEst, ".RData")
}
#save(out, file = file.path(saveDir, saveFile), compress="xz" )
save(out, file = file.path(saveDir, saveFile) )
if (verbose) {
cat("Output saved in:\n", file.path(saveDir, saveFile), "\n\n")
}
}
## it should not be an "either/or" decision whether you save or output the results
return( invisible( out ) )
}
######################################## end of 'monitorTrial'#########################################
mjuraska/seqDesign documentation built on Dec. 14, 2022, 4:26 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
################# Begin code for monitorTrial ################
## Some arguments described just below the argument list
#' Group Sequential Monitoring of Simulated Efficacy Trials for the Event of Potential Harm, Non-Efficacy, and High Efficacy
#'
#' \code{monitorTrial} applies a group sequential monitoring procedure to data-sets generated by \code{simTrial}, which may result in modification or termination of each simulated trial.
#'
#' @param dataFile if \eqn{\code{saveDir} = \code{NULL}}, a list returned by \code{simTrial}; otherwise a name (character string) of an \code{.RData} file created by \code{simTrial}
#' @param stage1 the final week of stage 1 in a two-stage fixed-follow-up trial
#' @param stage2 the final week of stage 2 in a two-stage fixed-follow-up trial, i.e., the maximum total follow-up time
#' @param harmMonitorRange a 2-component numeric vector specifying the range of the pooled number of events (pooled over the placebo and vaccine arm accruing events the fastest) over which the type I error rate, specified in \code{harmMonitorAlpha}, shall be spent (per vaccine arm). Note that \code{harmMonitorRange} does not specify a range over which potential-harm stopping boundaries will be computed; instead, it specifies when potential-harm monitoring will start, and the range over which the type I error rate \code{harmMonitorAlpha} will be spent. If \eqn{\code{nonEffStartMethod} = \code{"FKG"}} or \code{"fixed"}, then the second value is ignored and can be replaced with \code{NA}.
#' @param harmMonitorAlpha a numeric value (0.05 by default) specifying the overall type I error rate for potential-harm monitoring (per vaccine arm). To turn off potential-harm monitoring, set \code{harmMonitorAlpha} equal to 0.00001.
#' @param alphaPerTest a per-test nominal significance level for potential-harm monitoring. If \code{NULL} (default), a per-test significance level is calculated that yields the overall type I error rate of \code{harmMonitorAlpha} at the end of \code{harmMonitorRange}.
#' @param nonEffStartMethod a character string specifying the method used for determining when non-efficacy monitoring is to start. The default method of Freidlin, Korn, and Gray (2010) ("\code{FKG}") calculates the minimal pooled event count (pooled over the placebo and vaccine arm accruing events the fastest) such that a hazard-ratio-based VE point estimate of 0\% would result in declaring non-efficacy, i.e., the upper bound of the two-sided \eqn{(1-\code{alphaNoneff}) 100\%} confidence interval for VE based on the asymptotic variance of the log-rank statistic equals the non-efficacy threshold specified as component \code{upperVEnonEff} in the list \code{nonEffStartParams}. If this list component is left unspecified, the argument \code{upperVEnonEff} is used as the non-efficacy threshold. The alternative method ("\code{fixed}") starts non-efficacy monitoring at a fixed pooled event count (pooled over the placebo and vaccine arm accruing events the fastest) specified by component \code{N1} in the list \code{nonEffStartParams}.
#' @param nonEffStartParams a list with named components specifying parameters required by \code{nonEffStartMethod} (\code{NULL} by default)
#' @param nonEffIntervalUnit a character string specifying whether intervals between two adjacent non-efficacy interim analyses should be event-driven (default option "\code{counts}") or calendar time-driven (option "\code{time}")
#' @param nonEffInterval a numeric vector (a number of events or a number of weeks) specifying the timing of non-efficacy interim analyses. If a single numeric value is specified, then all interim looks are equidistant (in terms of the number of events or weeks), and the value specifies the constant increment of information or time between two adjacent interim looks. If a numeric vector with at least two components is specified, then, following the initial interim look, the timing of subsequent interim looks is determined by (potentially differential) increments of information or time specified by this vector.
#' @param nonEffCohorts a named list specifying all cohorts (for both timing and analysis) used in non-efficacy monitoring. The required list components characterize the 'timing cohort,' i.e., the cohort events in which determine the analysis timepoints in an event-driven design (components \code{times}, \code{timeUnit}, and \code{timingCohort}), and the analysis cohort(s), i.e., the cohort(s) in which inference is made about the null hypothesis of non-efficacy (component \code{cohort1} is required, and \code{cohort2}, etc. are optional for additional analysis cohorts if, e.g., non-efficacy monitoring is conducted in both the ITT and per-protocol cohorts). As for the timing cohort, \code{times} specifies the analysis timepoints in terms of event counts (\eqn{\code{timeUnit} = \code{"counts"}} is the only implemented option). \code{timingCohort} is a list characterizing the timing cohort by components \code{lagTime}, a lag time (in weeks) that controls event inclusion for timing (if no lag is desired, it can be set to \code{NULL} or 0), and \code{cohortInd}, a character string naming an indicator variable included in the data frames outputted by \code{\link{simTrial}}, which also controls event inclusion for the purpose of determining analysis timepoints. The other top-level list components \code{cohort1}, \code{cohort2}, etc. are each a list that must contain a component named \code{estimand}, which can take on values \code{"cox"}, \code{"cuminc"}, or \code{"combined"}. Optional list components in \code{cohort1}, \code{cohort2}, etc. are \code{lagTime}, \code{cohortInd}, \code{nullVE}, and \code{nominalAlphas}. \code{lagTime} specifies a lag time (in weeks) that controls event inclusion for the analysis cohort. If no lag time is desired, then this component can be ignored, set to \code{NULL}, or set to 0. \code{cohortInd} is a character string naming an indicator variable included in the data frames outputted by \code{\link{simTrial}} to be used to subset participants into the analysis cohort. This component allows inclusion of, e.g., a "per-protocol" variable (e.g., by setting \code{cohortInd} to \code{"pp1"}). \code{nullVE} specifies the one-sided null hypothesis as \eqn{H_0: VE \geq \code{nullVE}}. The rejection of \eqn{H_0} constitutes evidence for non-efficacy. \code{nominalAlphas} specifies nominal significance levels in a two-arm event-driven trial design.
#' @param effCohort a named list specifying the timing and analysis cohorts used in group-sequential efficacy monitoring (if part of the monitoring plan). The required list components characterize the 'timing cohort,' i.e., the cohort events in which determine the analysis timepoints in an event-driven design (components \code{times}, \code{timeUnit}, and \code{timingCohort}), and a single analysis cohort, i.e., the cohort in which inference is made about the null hypothesis of efficacy (components \code{estimand}, \code{lagTime}, \code{cohortInd}, \code{nullVE}, and \code{nominalAlphas}). As for the timing cohort, \code{times} specifies the analysis timepoints in terms of event counts (\eqn{\code{timeUnit} = \code{"counts"}} is the only implemented option). \code{timingCohort} is a list characterizing the timing cohort by components \code{lagTime}, a lag time (in weeks) that controls event inclusion for timing (if no lag is desired, it can be set to \code{NULL} or 0), and \code{cohortInd}, a character string naming an indicator variable included in the data frames outputted by \code{\link{simTrial}}, which also controls event inclusion for the purpose of determining analysis timepoints. As for the analysis cohort, \code{estimand} can take on values \code{"cox"}, \code{"cuminc"}, or \code{"combined"}. \code{lagTime} specifies a lag time (in weeks) that controls event inclusion for the analysis cohort. If no lag time is desired, then this component can be ignored, set to \code{NULL}, or set to 0. \code{cohortInd} is a character string naming an indicator variable included in the data frames outputted by \code{\link{simTrial}} to be used to subset participants into the analysis cohort. This component allows inclusion of, e.g., a "per-protocol" variable (e.g., by setting \code{cohortInd} to \code{"pp1"}). \code{nullVE} specifies the one-sided primary null hypothesis as \eqn{H_0: VE \leq \code{nullVE}}. The rejection of \eqn{H_0} constitutes evidence for clinically relevant efficacy. \code{nominalAlphas} specifies nominal significance levels in a two-arm event-driven trial design.
#' @param lowerVEnoneff specifies an additional criterion for declaring non-efficacy if the hypothesis test is based on Wald confidence interval(s). It requires that the lower bound of the two-sided Wald CI(s) for the VE estimand(s), at the confidence level determined by \code{nonEffCohorts$cohort1$nominalAlphas}, etc., lie(s) below \code{lowerVEnoneff} (typically set equal to 0). If \code{NULL} (default), this criterion is ignored.
#' @param highVE specifies a criterion for declaring high-efficacy: the lower bound of the two-sided \eqn{(1-\code{alphaHigh}) 100\%} confidence interval for the VE estimand lies above \code{highVE} (typically a number in the 0.5--1 range). To turn off high efficacy monitoring, set \code{highVE} equal to 1.
#' @param stage1VE specifies a criterion for advancement of a treatment's evaluation into Stage 2: the lower bound of the two-sided \eqn{(1-\code{alphaStage1}) 100\%} confidence interval for the VE estimand lies above \code{stage1VE} (typically set equal to 0)
#' @param lowerVEuncPower a numeric vector with each component specifying a one-sided null hypothesis \eqn{H_0: VE(0-\code{stage1}) \leq \code{lowerVEuncPower} \times 100\%}. Unconditional power (i.e., accounting for sequential monitoring) to reject each \eqn{H_0} is calculated, where the rejection region is defined by the lower bound of the two-sided \eqn{(1-\code{alphaUncPower}) 100\%} confidence interval for the VE estimand being above the respective component of \code{lowerVEuncPower} (typically values in the 0--0.5 range).
#' @param alphaHigh one minus the nominal confidence level of the two-sided confidence interval used for high efficacy monitoring
#' @param alphaStage1 one minus the nominal confidence level of the two-sided confidence interval used for determining whether a treatment's evaluation advances into Stage 2
#' @param alphaUncPower one minus the nominal confidence level of the two-sided confidence interval used to test one-sided null hypotheses \eqn{H_0: VE(0-\code{stage1}) \leq \code{lowerVEuncPower} \times 100\%} against alternative hypotheses \eqn{H_1: VE(0-\code{stage1}) > \code{lowerVEuncPower} \times 100\%}. The same nominal confidence level is applied for each component of \code{lowerVEuncPower}.
#' @param saveFile a character string specifying the name of the output \code{.RData} file. If \code{NULL} (default), a default file name will be used.
#' @param saveDir a character string specifying a path for \code{dataFile}. If supplied, the output is also saved as an \code{.RData} file in this directory; otherwise the output is returned as a list.
#' @param verbose a logical value indicating whether information on the output directory, file name, and monitoring outcomes should be printed out (default is \code{TRUE})
#'
#' @details All time variables use week as the unit of time. Month is defined as 52/12 weeks.
#'
#' Potential harm monitoring starts at the \code{harmMonitorRange[1]}-th infection pooled over the placebo group and the vaccine regimen that accrues infections the fastest. The potential harm analyses continue at each additional infection until the first interim analysis for non-efficacy. The monitoring is implemented with exact one-sided binomial tests of H0: \eqn{p \le p0} versus H1: \eqn{p > p0}, where \eqn{p} is the probability that an infected participant was assigned to the vaccine group, and \eqn{p0} is a fixed constant that represents the null hypothesis that an infection is equally likely to be assigned vaccine or placebo. Each test is performed at the same prespecified nominal/unadjusted alpha-level (\code{alphaPerTest}), chosen based on simulations such that, for each vaccine regimen, the overall type I error rate by the \code{harmMonitorRange[2]}-th arm-pooled infection (i.e., the probability that the potential harm boundary is reached when the vaccine is actually safe, \eqn{p = p0}) equals \code{harmMonitorAlpha}.
#'
#' Non-efficacy is defined as evidence that it is highly unlikely that the vaccine has a beneficial effect measured as VE(0--\code{stage1}) of \code{upperVEnoneff} x 100\% or more. The non-efficacy analyses for each vaccine regimen will start at the first infection (pooled over the vaccine and placebo arm) determined by \code{nonEffStartMethod}. Stopping for non-efficacy will lead to a reported two-sided (1-\code{alphaNoneff}) x 100\% CI for VE(0--\code{stage1}) with, optionally, the lower confidence bound below \code{lowerVEnoneff} and the upper confidence bound below \code{upperVEnoneff}, where \code{estimand} determines the choice of the VE(0--\code{stage1}) estimand. This approach is similar to the inefficacy monitoring approach of Freidlin, Korn, and Gray (2010). If \code{estimand = "combined"}, stopping for non-efficacy will lead to reported (1-\code{alphaNoneff}) x 100\% CIs for both VE parameters with, optionally, lower confidence bounds below \code{lowerVEnoneff} and upper confidence bounds below \code{upperVEnoneff}. If \code{laggedMonitoring = TRUE}, stopping for non-efficacy will lead to reported (1-\code{alphaNoneff}) x 100\% CIs for both VE(0--\code{stage1}) and VE(\code{lagTime}--\code{stage1}) with, optionally, lower confidence bounds below \code{lowerVEnoneff} and upper confidence bounds below \code{upperVEnoneff}.
#'
#' High efficacy monitoring allows early detection of a highly protective vaccine if there is evidence that VE(0--\code{stage2}) \eqn{>} \code{highVE} x 100\%. It is synchronized with non-efficacy monitoring during Stage 1, and a single high-efficacy interim analysis during Stage 2 is conducted halfway between the end of Stage 1 and the end of the trial. While monitoring for potential harm and non-efficacy restricts to \code{stage1} infections, monitoring for high efficacy counts all infections during \code{stage1} or \code{stage2}, given that early stopping for high efficacy would only be warranted under evidence for durability of the efficacy.
#'
#' The following principles and rules are applied in the monitoring procedure:
#' \itemize{
#' \item Exclude all follow-up data from the analysis post-unblinding (and include all data pre-unblinding).
#' \item The monitoring is based on modified ITT analysis, i.e., all subjects documented to be free of the study endpoint at baseline are included and analyzed according to the treatment assigned by randomization, ignoring how many vaccinations they received (only pre-unblinding follow-up included).
#' \item If a vaccine hits the harm boundary, immediately discontinue vaccinations and accrual into this vaccine arm, and unblind this vaccine arm (continue post-unblinded follow-up until the end of Stage 1 for this vaccine arm).
#' \item If a vaccine hits the non-efficacy boundary, immediately discontinue vaccinations and accrual into this vaccine arm, keep blinded and continue follow-up until the end of Stage 1 for this vaccine arm.
#' \item If and when the last vaccine arm hits the non-efficacy (or harm) boundary, discontinue vaccinations and accrual into this vaccine arm, and unblind (the trial is over, completed in Stage 1).
#' \item Stage 1 for the whole trial is over on the earliest date of the two events: (1) all vaccine arms have hit the harm or non-efficacy boundary; and (2) the last enrolled subject in the trial reaches the final \code{stage1} visit.
#' \item Continue blinded follow-up until the end of Stage 2 for each vaccine arm that reaches the end of \code{stage1} with a positive efficacy (as defined by \code{stage1VE}) or high efficacy (as defined by \code{highVE}) result.
#' \item If at least one vaccine arm reaches the end of \code{stage1} with a positive efficacy or high efficacy result, continue blinded follow-up in the placebo arm until the end of Stage 2.
#' \item Stage 2 for the whole trial is over on the earliest date of the two events: (1) all subjects in the placebo arm and each vaccine arm that registered efficacy or high efficacy in \code{stage1} have failed or been censored; and (2) all subjects in the placebo arm and each vaccine arm that registered efficacy or high efficacy in \code{stage1} have completed the final \code{stage2} visit.
#' }
#'
#' The above rules have the following implications:
#' \itemize{
#' \item If a vaccine hits the non-efficacy boundary but Stage 1 for the whole trial is not over, then one includes in the analysis all follow-up through the final \code{stage1} visit for that vaccine regimen, including all individuals accrued up through the date of hitting the non-efficacy boundary (which will be the total number accrued to this vaccine arm).
#' \item If a vaccine hits the harm boundary, all follow-up information through the date of hitting the harm boundary is included for this vaccine; no follow-up data are included after this date.
#' \item If and when the last vaccine arm hits the non-efficacy (or harm) boundary, all follow-up information through the date of hitting the non-efficacy (or harm) boundary is included for this vaccine; no follow-up data are included after this date.
#' }
#'
#' @return If \code{saveDir} (and, optionally \code{saveFile}) is specified, the output list (named \code{out}) is saved as an \code{.RData} file in \code{saveDir} (the path to \code{saveDir} is printed); otherwise it is returned. The output object is a list of length equal to the number of simulated trials, each of which is a list of length equal to the number of treatment arms, each of which is a list with (at least) the following components:
#' \itemize{
#' \item \code{boundHit}: a character string stating the monitoring outcome in this treatment arm, i.e., one of \code{"Harm"}, \code{"NonEffInterim"}, \code{"NonEffFinal"}, \code{"Eff"}, or \code{"HighEff"}. The first four outcomes can occur in Stage 1, whereas the last outcome can combine data over Stage 1 and Stage 2.
#' \item \code{stopTime}: the time of hitting a stopping boundary since the first subject enrolled in the trial
#' \item \code{stopInfectCnt}: the pooled number of infections at \code{stopTime}
#' \item \code{summObj}: a \code{data.frame} containing summary information from each non-/high efficacy interim analysis
#' \item \code{finalHRci}: the final CI for the hazard ratio, available if \code{estimand!="cuminc"} and there is at least 1 infection in each arm
#' \item \code{firstNonEffCnt}: the number of infections that triggered non-efficacy monitoring (if available)
#' \item \code{totInfecCnt}: the total number of \code{stage1} (\code{stage2} if \code{boundHit = "HighEff"}) infections
#' \item \code{totInfecSplit}: a table with the numbers of \code{stage1} (\code{stage2} if \code{boundHit = "HighEff"}) infections in the treatment and control arm
#' \item \code{lastExitTime}: the time between the first subject's enrollment and the last subject's exiting from the trial
#' }
#'
#' @references Freidlin B., Korn E. L., and Gray R. (2010), A general inefficacy interim monitoring rule for randomized clinical trials. \emph{Clinical Trials} 7(3):197-208.
#'
#' @examples
#' simData <- simTrial(N=c(1000, 1000), aveVE=c(0, 0.4),
#' VEmodel="half", vePeriods=c(1, 27, 79), enrollPeriod=78,
#' enrollPartial=13, enrollPartialRelRate=0.5, dropoutRate=0.05,
#' infecRate=0.06, fuTime=156, visitSchedule=seq(0, 156, by=4),
#' missVaccProb=0.05, VEcutoffWeek=26, nTrials=5,
#' stage1=78, randomSeed=300)
#'
#' ### trial design: fixed follow-up; cumulative incidence-based VE estimand;
#' no efficacy monitoring; non-efficacy monitoring using
#' the Freidlin et al. method; hypothesis tests based on Wald
#' confidence intervals
#' ### RIGHT NOW, THE BELOW CALL IS BROKEN
#' monitorData <- monitorTrial(dataFile=simData, stage1=78, stage2=156,
#' harmMonitorRange=c(10, NA), harmMonitorAlpha=0.05,
#' nonEffStartMethod="FKG", nonEffInterval=20,
#' nonEffCohorts=list(timeUnit="counts",
#' # right now 'timingCohort' is required
#' timingCohort=list(lagTime=0),
#' cohort1=list(estimand="cuminc",
#' nullVE=0.4,
#' nominalAlphas=0.025)),
#' # it appears that right now 'effCohort' must be specified
#' # even when there is no efficacy monitoring;
#' # this call of monitorTrial() still doesn't run
#' # because it requires event counts for timing
#' effCohort=list(timeUnit="counts",
#' timingCohort=list(lagTime=0)),
#' lowerVEnoneff=0, highVE=0.7, lowerVEuncPower=0,
#' alphaHigh=0.05, alphaUncPower=0.05)
#'
#' ### trial design: event-driven; hazard-based VE estimand;
#' harmonized efficacy and non-efficacy monitoring;
#' hypothesis tests using the score test in the Cox model
#' ### THE SCORE TEST OPTION NEEDS TO BE ADDED
#' monitorData <- monitorTrial(dataFile=simData, stage1=78, stage2=156,
#' harmMonitorRange=c(10, 50), harmMonitorAlpha=0.05,
#' nonEffCohorts=list(
#' times=c(50, 100, 150),
#' timeUnit="counts",
#' timingCohort=list(lagTime=0),
#' cohort1=list(estimand="cox",
#' nullVE=0.4,
#' nominalAlphas=c(0.0001, 0.0060, 0.0231))),
#' effCohort=list(times=c(50, 100, 150),
#' timeUnit="counts",
#' timingCohort=list(lagTime=0),
#' estimand="cox",
#' nullVE=0,
#' nominalAlphas=c(0.0001, 0.0060, 0.0231)),
#' highVE=1, lowerVEuncPower=0,
#' alphaHigh=0.05, alphaUncPower=0.05)
#'
#' ### alternatively, to save the .RData output file (no '<-' needed):
#' ###
#' ### simTrial(N=c(1000, 1000), aveVE=c(0, 0.4),
#' ### VEmodel="half", vePeriods=c(1, 27, 79), enrollPeriod=78,
#' ### enrollPartial=13, enrollPartialRelRate=0.5, dropoutRate=0.05,
#' ### infecRate=0.06, fuTime=156, visitSchedule=seq(0, 156, by=4),
#' ### missVaccProb=0.05, VEcutoffWeek=26, nTrials=5,
#' ### stage1=78, saveDir="./", randomSeed=300)
#' ###
#' ### THIS CALL NEEDS TO BE REVISED TO A SIMPLE BUT WORKING CODE
#' ### THE INTENT HERE IS ONLY TO ILLUSTRATE HOW OUTPUT FILES CAN BE READ IN
#' ### THE PURPOSE IS NOT TO SHOW ANY ADDITIONAL TRIAL DESIGNS
#' ### monitorTrial(dataFile=
#' ### "simTrial_nPlac=1400_nVacc=1000_1000_aveVE=0.2_0.4_infRate=0.04.RData",
#' ### stage1=78, stage2=156, harmMonitorRange=c(10,100), alphaPerTest=NULL,
#' ### nonEffStartMethod="FKG", nonEffInterval=20, lowerVEnoneff=0,
#' ### upperVEnoneff=0.4, highVE=0.7, stage1VE=0, lowerVEuncPower=0,
#' ### alphaNoneff=0.05, alphaHigh=0.05, alphaStage1=0.05, alphaUncPower=0.05,
#' ### estimand="cuminc", lagTime=26, saveDir="./")
#'
#' @seealso \code{\link{simTrial}}, \code{\link{censTrial}}, \code{\link{rankTrial}}, \code{\link{estHRbound}}, \code{\link{crossBoundCumProb}}, and \code{\link{decisionTimes}}
#'
#' @export
monitorTrial <- function (dataFile,
stage1,
stage2,
## range over which to "spend" the type-I error specified in argument
## 'harmMonitorRange'. It should be a vector of length 2 giving the
## (start, stop) infection-count range over which the type-I error
## will be spent. Please note that this argument does *NOT* dictate
## when harm monitoring will END. The range dictates when it will
## START, and over what range the type-I error will be spread. If
## you plan to use nonEffStartMethod "FKG" or "fixed" without lagged
## Monitoring (i.e. laggedMonitoring=FALSE) then the 'stop' value of
## the range, if provided, will not be used and need not be specified
## or may be replaced with an NA.
harmMonitorRange,
## Total Type I error for potential harm monitoring (per vacc. arm)
harmMonitorAlpha=0.05,
## 'harmMonitorContol' is a list containing other parameters controlling
## the harm monitoring. Currently the only one used in 'maxCnt'. This
## *MUST* be specified and represesents the set of infection counts for
## which harm bounds will be constructed. Harm monitoring ends when
## 'maxCnt' has been exceeded, whether or not you've reached the
## criteria for initiaton of non-efficacy monitoring, so you should:
## (a) choose a nonEff start method that guarantees starting
## at/before 'maxCnt', or
## (b) choose maxCnt to be a value larger than you're non-efficacy
## starting count will ever be (maybe use 2 or 3 times the
## upper range for harm monitoring (harmMonitorRange[2]).
## If maxCnt is left unspecified, the default value of
## 3 x harmMonitorRange[2] will be used for it.
#harmMonitorControl=list( maxCnt=NULL),
## if you have determined the constant alpha value to do the binomial
## 'potential-harm-monitoring' tests at, they you can pass it through
## this argument, in which case argument 'harmMonitorAlpha' need not
## be specified
alphaPerTest=NULL,
## character vector of methods to use to determine when to start non-efficacy
## monitoring. Each method requires different input information, and that
## information should be input via the 'nonEffStartParams' argument - which
## should be a list that the specifies the inputs for the method. Each
## method's input requirements should be specified in the documentation.
##
## Methods:
## FKG - the starting method suggested in Freidlin, Korn and Gray's 2010
## 'Clinical Trials' paper: "A general inefficacy interim
## monitoring rule for randomized clinical trials"
## It boils down to start monitoring at the earlier infection count
## such that an estimated effect <= 0 would cause the trial to
## stop. In our code that translates to a 95% CI (based on the
## asymptotic variance of the log-rank statistic) around an
## estimated VE of 0% would exclude the VE specified by parameter
## 'upperVEnonEff'. The expectation is that this will be harmonized
## with the non-efficacy monitoring, and so will use the arguments
## 'upperVEnonEff' 'alphaNoneff' from argument list. However, should
## you want to use different values, you may do say by passing them
## via the 'nonEffStartParams' argument list. The same naming should
## be used within the list.
## Parameters: upperVEnonEff alphaNoneff
##
## fixed - Starts at the infection count specified by paramter 'N1'.
## If 'N1' is set to 75, then non-efficacy always begins at the
## 75th infection
## Parameter: N1
##
nonEffStartMethod=c("FKG", "fixed"),
## A *list* (not a vector) of parameters needed by the method specified
## in argument 'nonEffStartMethod'. Some methods have defaults in place,
## for those you do not need to use 'nonEffStartParams' unless you want
## values other than the defaults.
nonEffStartParams=NULL,
nonEffIntervalUnit=c("counts","time"),
# a single numeric value (for constant increments) or a numeric vector
# (for variable increments)
nonEffInterval,
# if noneff monitoring is being used (which it always is), then
# nonEffCohorts must provide info on each cohort being monitored.
# the argument is a list of lists. The internal lists should be
# named "cohort1", "cohort2", ..., "cohort<N>".
# Each cohort's list must contain a component named 'estimand'
# specifying either 'cox','cuminc', or 'combined' (for both).
# Optional list components are:
# 'lagTime' - indicates the lag time that controls event inclusion
# for the cohort. If no lag time is desired, then this
# component can be excluded, set to NULL, or set to 0
# 'cohortInd'- Character string naming an indicator variable to be
# used to subset participants to the desired cohort.
# The variable must already exist in the simTrials
# output. This allows inclusion of e.g a "per-protocol"
# variable
#
# Example specifications:
# - 2 cohorts, nominal CI monitoring
# nonEffCohorts <- list(
# times=NULL, timeUnit="counts",
# timingCohort=list(lagTime=NULL, cohortInd=NULL),
# cohort1 = list(
# estimand="cox", lagTime=NULL, nullVE=0.5,
# nominalAlphas=0.05
# ),
# cohort2 = list(
# estimand="cox", lagTime=6, cohortInd="pp1",
# nullVE=0.5, nominalAlphas=0.05)
# )
# )
#
# - 1 cohort, sequential control of type-I error
# nonEffCohorts <- list(
# times=c(50,75,100), timeUnit="counts",
# timingCohort=list(lagTime=NULL, cohortInd=NULL),
# cohort1 = list(
# estimand="cox", lagTime=NULL, cohortInd="pp1",
# nullVE=0.6, nominalAlphas=c(0.0030, 0.0183, 0.0440)
# ),
# )
#
## The default set-up (below) monitors one cohort, with 'cox' estimand
## and lagTime defaulting to NULL
nonEffCohorts = list(
times=NULL, timeUnit="counts",
timingCohort=list(lagTime=NULL, cohortInd=NULL),
cohort1 = list(
estimand="cox", lagTime=NULL, cohortInd=NULL,
nullVE=NULL, nominalAlphas=NULL
)
),
# timeLag = lag to use in converting 'times' counts into time
# (not needed if timeUnit = "time" or if no time lag is used)
#
# -- These args are only needed if using sequential monitoring bounds
# + totalAlpha = total (two-sided) alpha for efficacy testing
# + totalEvents = number of events expected (or planned for) at the
# last analysis time. Needed to convert to 'information time'
# for use of spending function approach
# + spendingFunction = Tells which kind of bound to create
# Specification is either: (1) the text string "OBF", indicating
# use of O'Brien-Fleming bounds, or (2) an alpha spending function.
# The alpha spending function (ASF) 'f' must satisfy the
# following conditions:
# (a) takes a single argument with values in [0,1],
# (b) returns values in [0,1]
# (c) must be strictly increasing
# (d) f(0)=0 and f(1)=1
# We provide a function asf() which allows creation of such
# functions for some common families of bounds.
# asf("Pocock") creates a pocock ASF,
# asf("power", phi=x) creates an ASF from the power family
# with parameter phi set to 'x'.
# asf("HSD", phi=x) creates an ASF from the Hwang-Shih-DeCani
# Family with parameter phi set to 'x'.
#
# More optional components are:
# 'lagTime' - indicates the lag time that controls event inclusion
# for the cohort. If no lag time is desired, then this
# component can be excluded, set to NULL, or set to 0
# 'cohortInd'- Character string naming an indicator variable to be
# used to subset participants to the desired cohort.
# The variable must already exist in the simTrials
# output. This allows inclusion of e.g a "per-protocol"
# variable
# NOT ALL IMPLEMENTED YET
#
# Minimum required:
# times, timeUnit, estimand, lagTime, totalAlpha, nullVE,
# totalEvents
# effCohort = list( times, timeUnit, timeLag, timeCohortInd,
# estimand, lagTime, cohortInd,
# totalAlpha, nullVE, totalEvents,
# spendingFunction),
effCohort = list( times=NULL, timeUnit="counts",
timingCohort=list(lagTime=NULL, cohortInd=NULL),
nullVE = NULL,
estimand="cox", lagTime=NULL, cohortInd=NULL,
nominalAlphas=NULL),
## lowerVEnoneff is not required. Specify only if you want this
## condition as part of your monitoring.
lowerVEnoneff=NULL,
#upperVEnoneff,
highVE,
stage1VE,
lowerVEuncPower=NULL,
#alphaNoneff,
alphaHigh,
alphaStage1,
alphaUncPower=NULL,
saveFile= NULL,
saveDir = NULL,
verbose = TRUE ) {
## selected Arguments
##
## lowerVEnoneff:
## lower confidence bound to be below 'lowerVEnoneff' to meet criterion 1
## for declaring non-efficacy
##
## upperVEnoneff:
## upper confidence bound to be below 'upperVEnoneff' to meet criterion 2
## for declaring non-efficacy
##
## highVE:
## lower confidence bound to be above 'highVE' for declaring high efficacy
##
## stage1VE:
## lower confidence bound to be above 'stage1VE' for advancing into Stage 2
##
## lowerVEuncPower:
## vector of VE value for which the user wishes to determine unconditional
## power to reject H0: VE <= (lowerVEuncPower * 100)%
## at 2-sided alpha level 'alphaUncPower'
##
## alphaNoneff:
## One minus confidence level of 2-sided CI for non-efficacy monitoring
##
## alphaHigh:
## One minus confidence level of 2-sided CI for high efficacy monitoring
##
## alphaStage1:
## One minus confidence level of 2-sided CI for testing whether an arm
## advances into Stage 2
##
## alphaUncPower:
## One minus confidence level of 2-sided CI for unconditional power to
## reject H0: VE <= (lowerVEuncPower * 100)%
nonEffIntervalUnit <- match.arg(nonEffIntervalUnit)
#nonEffTimeUnit <- match.arg(nonEffTimeUnit)
nonEffStartMethod <- match.arg(nonEffStartMethod)
if (!is.null(nonEffCohorts$times) ) {
## disable start method, since not needed
nonEffStartMethod <- ""
}
if ( is.list(dataFile) ) {
trialObj <- dataFile
} else {
if ( !is.null(saveDir) ){
## load in RData object (a list named 'trialObj' )
load(file.path(saveDir, dataFile))
} else {
load(dataFile)
}
}
## check contents of 'trialData'
if ( !all( c("entry","exit","event","trt") %in% names(trialObj$trialData[[1]]) ) )
stop("Trial data must contain columns: entry, exit, event, trt\n")
nTrtArms <- as.integer( trialObj$nArms - 1 )
nTrials <- length(trialObj[["trialData"]])
## set altVE based on lowerVEuncPower, etc. if those were specified instead
alphaAltVE <- altVE <- NULL
if ( is.null(altVE) && !is.null(lowerVEuncPower)) {
altVE <- lowerVEuncPower
}
if (is.null(alphaAltVE) && !is.null(altVE)) {
alphaAltVE <- ifelse(!is.null(alphaUncPower), alphaUncPower, alphaStage1)
}
## ----------------------------------------------------------------------
## The following code is needed to determine infection count at which to
## begin the nonEff monitoring. Also needed for harm-monitoring
## derive the prob. of assignment to each vaccine arm (relative to placebo
## only, not relative to other vaccine arms)
null.p <- numeric(length=nTrtArms)
for (i in 1:nTrtArms) {
arms.i <- c(1, i+1)
assn.probs <- trialObj$trtAssgnProbs[ arms.i ]
null.p[i] <- assn.probs[2] / sum(assn.probs)
}
## if we have multiple vaccine groups and they are of different sizes,
## then we must stop. The code is not set up to handle that yet.
if (nTrtArms > 1 && ( max(null.p)-min(null.p) > sqrt(.Machine$double.eps) ) )
{
stop("The code is not currently set up to handle harm monitoring when\n",
"there are multiple treatment arms and their sizes differ. \n\n")
} else {
## all values of null.p are equal, just choose the first
null.p <- null.p[1]
}
## This will also determine when nonEfficacy monitoring starts
## (and therefore how long harm monitoring will last).
## Two methods are found here and two later on. These two do not
## require any specific information about the trial to determine,
## whereas the later methods do.
## Method implemented here:
## Start nonEff monitoring at the smallest infection count for which we
## would be able to reject the null hypothesis that VE = upperVEnoneff
## if our point estimate of VE = 0. This computation requires us to
## assume that our test statistic is asymptotically equivalent to the
## log-rank statistic, which (under the null of no diff - not the setting
## we're in, right?) has a distribution that is
## Normal( 0, sqrt(1/(p*(1-p)*D))), where
## D = number of endpoints between the placebo and trt arm being monitored,
## p = prob. of assignment to the vaccine arm vs. the placebo arm (all
## other arms ignored for computing 'p')
##
## if 'timeUnit' is not specified but 'timingCohort' is, assume timeUnit is 'counts'
if ( is.null(nonEffCohorts$timeUnit) ) {
if ( !is.null(nonEffCohorts$timingCohort) ) {
nonEffCohorts$timeUnit <- "counts"
} else {
stop("You must specify 'timeUnit' within 'nonEffCohorts' as 'time' or specify",
" 'timingCohort' to provide info needed to translate counts into times\n\n")
}
}
## Create indicators of whether either a lagTime or cohortInd has been
## specified for timing of non-efficacy analyses
if ( nonEffCohorts$timeUnit != "time" ) {
if ( !is.null(nonEffCohorts$timingCohort$lagTime) &&
nonEffCohorts$timingCohort$lagTime > 0 ) {
laggedMonitoring <- TRUE
} else {
laggedMonitoring <- FALSE
}
neTimeCohInd <- ifelse( !is.null(nonEffCohorts$timingCohort$cohortInd),
TRUE, FALSE )
## if either indicator is set (is "ne" timing diff. from Harm timing
neTimingDiff <- ifelse(laggedMonitoring || neTimeCohInd, TRUE, FALSE )
}
if ( !is.null(nonEffCohorts$times) && nonEffCohorts$timeUnit == "counts") {
nonEffCnts <- nonEffCohorts$times
if ( !neTimingDiff ) {
N1 <- nonEffCnts[1]
N1.ne <- N1
} else {
N1.ne <- nonEffCnts[1]
}
}
if ( nonEffStartMethod == "FKG") {
if (!is.null(nonEffStartParams)) {
## values specified directly by user (maybe - we don't
## know what's actually in 'nonEffStartMethod' yet
upperVE <- nonEffStartParams$upperVEnoneff
alphaNE <- nonEffStartParams$alphaNoneff
if ( is.null(upperVE) && is.null(alphaNE) )
warning("The argument 'nonEffStartParams' has been specified but",
"it does not contain\n", "parameters used by method FKG\n\n",
immediate.=TRUE)
} else {
upperVE <- NULL
alphaNE <- NULL
}
## Fill in NULL values for either parameter using the monitorTrial
## arguments upperVEnoneff and alphaNoneff
if ( is.null(upperVE) )
upperVE <- nonEffCohorts$cohort1$nullVE
if ( is.null(alphaNE) )
alphaNE <- nonEffCohorts$cohort1$nominalAlphas
if ( is.null(upperVE) || is.null(alphaNE) )
stop( "One or both of the arguments needed by nonEffStartMethod 'FKG'",
"are missing\n.", "They were not specified via argument ",
"'nonEffStartParams' and their values also could not be obtained\n",
"from the arguments 'nonEffCohorts$cohort1$nullVE' and ",
" 'nonEffCohorts$cohort1$nominalAlphas'\n",
"Exiting...\n" )
N1.ne <- ceiling(qnorm(1 - alphaNE/2)^2 /( null.p*(1-null.p)*log(1-upperVE)^2 ))
if ( !neTimingDiff ) {
N1 <- N1.ne
}
} else if ( nonEffStartMethod == "fixed" ) {
if (is.null( nonEffStartParams$N1 ) )
stop("Argument 'nonEffStartMethod' was specified as 'fixed'. ",
"This method requires\n", "you to provide an argument named",
"'N1' via the list argument 'nonEffStartParams'\n.",
"e.g. nonEffStartParams <- list(N1=60). Please fix.\n" )
N1.ne <- nonEffStartParams$N1
if (!neTimingDiff ) {
N1 <- N1.ne
}
}
## N1 should only exist if laggedMonitoring is FALSE, while N1.ne exists
## otherwise (N1.ne also exists if N1 does).
## Harm monitoring will end at N1 infections, so when "N1" exists, set the upper
## limit of 'harmMonitorRange' to N1 so boundaries will be defined based on it
if ( exists("N1") ) {
harmMonitorRange <- c(harmMonitorRange[1], N1)
}
## -------------------------- get harm bounds ---------------------------
## do a check on 'harmMonitorRange' to be sure it's correctly specified
if ( length(harmMonitorRange) != 2 || is.na(harmMonitorRange[2]) ||
is.infinite(harmMonitorRange[2]) ) {
if (neTimingDiff) {
stop("When the non-efficacy timing cohort differs from the harm monitoring\n",
"cohort, then the harmMonitoringRange must be fully specified - i.e. it\n",
"must be length 2 vector with both values being integers.\n\n")
} else {
stop("When nonEffStartMethod is NOT one of 'FKG' or 'fixed', then ",
"harmMonitoringRange must be fully specified.\n", "It should be a ",
"length 2 vector with both values being integers.\n\n")
}
}
# calculate stopping boundaries for harm
## NOTE: harmMonitorRange dictates that range over which the type-I error is
## spent, it does not specify the time at which harm monitoring stops. That is
## dictated by when non-eff monitoring begins
if ( is.null(alphaPerTest) ) {
## choose the value of alphaPerTest to spend the type I error over harmMonitorRange
alphaPerTest <- getAlphaPerTest(harmMonitorRange, null.p,
totalAlpha = harmMonitorAlpha)
}
## For generating the bounds, we need both 'N' and the upper limit of
## 'harmMonitorRange' to be an upper limit on the number of tests to do
## (i.e. an upper limit on when the first non-efficacy will be done)
if (exists("N1")) {
maxCnt = N1
} else if ( !exists("maxCnt") ) {
maxCnt = 2 * harmMonitorRange[2]
}
## get harm bounds
harmBounds <- getHarmBound(
N = maxCnt,
per.test = alphaPerTest,
harmBoundRange = c(harmMonitorRange[1], maxCnt),
null.p = null.p,
dataDir = saveDir,
verbose = verbose)
## -------------------------- end get harm bounds ---------------------------
## creates a list of length 'nTrials' each element of which is a list of
## length 'nTrtArms'
out <- rep( list(vector("list",nTrtArms)), nTrials )
## censor all trials to stage 1 (fastest to do all at once). Still need to
## keep access to uncensored data, for highEff monitoring
if ( !missing(stage1) ) {
stg1dat <- censorTrial(trialObj$trialData, times=stage1,
timeScale="follow-up")
}
## Make note of whether 'N1' exists (it will for certain types of non-eff. monitoring
## timing). If it does, then N1 will remain constant throughout the trial.
## Otherwise we'll need to reset it at the beginning of each new trial/arm
constant.N1 <- ifelse( exists("N1"), TRUE, FALSE )
## Start looping over trials
for (i in 1:nTrials ) {
## extract data for the i-th trial
if ( !missing(stage1) ) {
datI <- stg1dat[[i]]
datIall <- trialObj[["trialData"]][[ i ]]
} else {
datI <- trialObj[["trialData"]][[ i ]]
}
## create separate set of only infections ('events'), then
## order the events by trial time at which they are observed
eventDF <- datI[datI$event == 1, ]
eventDF <- eventDF[order(eventDF$exit), ]
#cat("Trial :", i, " \n")
## Now start comparisons - each active arm vs. placebo arm, one at a time
for (j in 1:nTrtArms) {
## subset *all data* down to the two arms being compared
datI.j <- datI[datI$trt %in% c(0,j), ]
## convert 'trt' to indicator variable before passing it to
## 'applyStopRules' (i.e. convert the non-zero values to 1)
datI.j$trt <- as.integer(datI.j$trt > 0 )
## same for the 'all' data.frame
if ( exists("datIall") ){
datIall.j <- datIall[datIall$trt %in% c(0,j), ]
datIall.j$trt <- as.integer(datIall.j$trt > 0 )
}
## 1. do harm monitoring for j-th treatment arm vs. placebo
## if we don't have a constant value of N1 (across trials), then
## remove the existing value of N1
if ( !constant.N1 && exists("N1") )
rm("N1")
## subset events in relevant arms: j-th active trt and placebo(trt=0)
E.j <- eventDF[eventDF$trt %in% c(0,j), ]
nInfec <- nrow(E.j) # counts infections through 'stage1'
E.j$nInf <- 1:nInfec
## set N1.ne for the cases for which N1 was just defined (immed. above)
if ( !exists("N1.ne") && exists("N1") )
N1.ne <- N1
## ---------------- Process nonEff timing args -------------------
## We need to process the timing information for non-efficacy monitoring now,
## because it's used to determine when harm monitoring will end
## The code is LONG :(
## determine end of stage 1 time, needed if any specified monitoring
## times/counts are not attained in stage 1.
## [ endStg1 == calendar time at which stage 1 follow-up ends ]
endStg1 <- max( datI.j$exit )
if ( neTimingDiff ) {
if (neTimeCohInd) {
timeCoh <- ( datI.j[[nonEffCohorts$timingCohort$cohortInd]]==1 )
nedatI.j <- datI.j[ timeCoh, ]
} else {
nedatI.j <- datI.j
}
if (laggedMonitoring) {
nedatI.j <- censorTrial( nedatI.j, times=nonEffCohorts$timingCohort$lagTime,
timeScale="follow-up", type="left")
}
nInfec.ne <- sum( nedatI.j$event == 1 )
}
## create new variable, 'nInfec_nonEff' that will contain the total number of
## infections under the chosen non-efficacy monitoring cohort
nInfec_nonEff <- ifelse(neTimingDiff, nInfec.ne, nInfec)
## Determine how many non-eff monitoring cohorts were defined
n.neCohorts <- sum( substr(names(nonEffCohorts),1,6) == "cohort" )
## initialize 'alphaNoneff.ij', a temp list containing the nominal alpha
## values for each non-efficacy cohort. This list we be modified, as req.
## based on the current simulated trial, which is why we reinitialize it
## for trial
alphaNoneff.ij <- lapply(1:n.neCohorts, function(idx, neC) {
neC[[paste0("cohort",idx)]]$nominalAlphas
}, neC=nonEffCohorts)
## determine times of non-efficacy testing (times, not counts), then monitor
if (is.null(nonEffCohorts$times)) {
if (nonEffIntervalUnit == "counts") {
## makes sequence of counts at which analyses will be done
if (length(nonEffInterval)==1){
# then a constant increment in the endpoint count is assumed
nonEffCnts <- seq(from = N1.ne, to = nInfec_nonEff, by = nonEffInterval)
} else {
# then 'nonEffInterval' is a vector of potentially variable increments
# specifying endpoint counts for subsequent interim looks following the
# initial look
nonEffCnts <- cumsum(c(N1.ne, nonEffInterval))
}
## Convert counts into times.
if ( neTimingDiff ) {
## subsetted and lagged (if required) dataset 'nedatI.j' was already created earlier
nonEffTimes.ij <- getInfectionTimes(nedatI.j, cnts=nonEffCnts, lagTime=NULL )
} else {
nonEffTimes.ij <- getInfectionTimes(datI.j, cnts=nonEffCnts)
}
} else {
## User specified a time-sequence for monitoring. Figure out when the first
## time and last times will be, then create the sequence
firstLastnonEffTimes <-
getInfectionTimes(datI.j, cnts=c(N1.ne, nInfec_nonEff),
lagTime=nonEffTimeLag )
if (length(nonEffInterval)==1){
nonEffTimes.ij <-
seq(from = firstLastnonEffTimes[1],
to = firstLastnonEffTimes[2],
by = nonEffInterval)
} else {
# then 'nonEffInterval' is a vector of potentially variable time increments
nonEffTimes.ij <- cumsum( c(firstLastnonEffTimes[1], nonEffInterval) )
## truncate so that nonEffTimes don't go beyond the trial duration
nonEffTimes.ij <- nonEffCohorts$times[ nonEffCohorts$times <= firstLastnonEffTimes[2] ]
}
}
} else {
## if nonEffTimes was specified
if (nonEffCohorts$timeUnit == "time") {
nonEffTimes.ij <- sort(nonEffCohorts$times)
## check if any specified times are beyond the extent of stage 1
if ( any(nonEffTimes.ij > endStg1) ) {
wTooLrg <- which( nonEffTimes.ij > endStg1 )
## if just one is too large, then set that last time to 'endStg1',
## no need to modify 'alpha's
if ( length(wTooLrg) == 1 ) {
nonEffTimes.ij[wTooLrg] <- endStg1
} else if (length(wTooLrg)>1) {
## If 'nominalAlphas' is length 1 for all cohorts, then we don't need
## modify nominal alphas. We set the time of the first missed test to
## the end of stage 1 fu time and drop the remaining tests
if ( all( sapply(alphaNoneff.ij, length) == 1 )) {
nonEffTimes.ij[ wTooLrg[1] ] <- endStg1
nonEffTimes.ij <- nonEffTimes.ij[ 1:wTooLrg[1] ]
} else {
## If not all are length 1, Issue warning but continue anyway
cat("Warning: More than one of the specified non-efficacy test times was ",
"not reached.\n", "The code cannot adjust the nominal alphas properly",
" for this.\n", "A final test will be inserted at the end of stage1 ",
"follow-up and will be carried\n", "out using the nominal alpha ",
"level associated with the final scheduled test\n")
## replace the first missed test with a test at the end of stage 1
## using the nominalAlpha level corresponding to final analysis
nonEffTimes.ij[ wTooLrg[1] ] <- endStg1
nonEffTimes.ij <- nonEffTimes.ij[ 1:wTooLrg[1] ]
## modify alpha for any cohorts that have more than one specified
for (idx in 1:n.neCohorts) {
if ( (L <- length(alphaNoneff.ij[[idx]]) ) > 1 )
alphaNoneff.ij[[idx]] <- alphaNoneff.ij[[idx]][ c(1:(wTooLrg[1]-1), L) ]
}
}
}
}
} else {
if ( neTimingDiff ) {
## subsetted and lagged (if required) dataset 'nedatI.j' was already created earlier
nonEffTimes.ij <- getInfectionTimes(nedatI.j, cnts=nonEffCnts)
} else {
nonEffTimes.ij <- getInfectionTimes(datI.j, cnts=nonEffCnts)
}
## check if we "lost" any tests, due to the counts not being reached
if ( length(nonEffTimes.ij) < length(nonEffCnts) ) {
nLost <- length(nonEffCnts) - length(nonEffTimes.ij)
## if we only "lost" 1 test, or if only one alpha specifed for each noneff
## cohort, then we don't need to modify the alphaNoneff.ij values
if (nLost==1 || all( sapply( alphaNoneff.ij, length) == 1) ) {
nonEffTimes.ij <- c(nonEffTimes.ij, endStg1)
} else {
## if we "lost" >1 test, or if have > 1 alpha specifed for any non-eff
## cohort, then we *do* need to modify some alphaNoneff.ij values
## Issue warning but continue anyway
cat("Warning: More than one of the specified non-efficacy test times was",
" not reached.\n", "The code cannot adjust the nominal alphas properly",
" for this.\n", "A final test will be inserted at the end of stage1 ",
"follow-up and will be carried\n", "out using the nominal alpha level",
" associated with the final scheduled test\n\n")
## add a test at the end stage 1 follow-up. Perform the test using
## using the 'alphaLevelNoneff' value corresponding to final analysis
nonEffTimes.ij <- c(nonEffTimes.ij, endStg1)
## modify alpha for any cohorts that have more than one specified
for (idx in 1:n.neCohorts) {
if ( (L <- length(alphaNoneff.ij[[idx]]) ) > 1 ){
alphaNoneff.ij[[idx]] <-
alphaNoneff.ij[[idx]][ c(1:(length(nonEffTimes.ij)-1), L) ]
}
}
}
}
}
}
## ---------------- End Processing nonEff timing args -------------------
## ----- Process efficacy analysis inputs here so we know what test --------
## ----- to perform if we stop early for harm or non-eff --------
## Now get 'N1' for the cases where we don't have it, which should only
## be when laggedMonitoring was specified
if ( !exists("N1") && neTimingDiff ) {
## determine the total infection count (based on full set of data) that corresponds
## to the first non-eff test time
N1 <- sum( datI.j$event==1 & datI.j$exit <= nonEffTimes.ij[1] )
}
## Ensure we have harm bounds up to 'N1', if not, regenerate them
if (N1 > max(harmBounds$N)) {
maxCnt <- round(1.5 * N1 )
## get harm bounds
harmBounds <- getHarmBound(
N = maxCnt,
per.test = alphaPerTest,
harmBoundRange = c(harmMonitorRange[1], maxCnt),
null.p = null.p,
dataDir = saveDir,
verbose = verbose)
}
## subset to restrict harm monitoring to the first 'N1' infections
harmBounds.j <- harmBounds[harmBounds$N <= N1, c("N","V","P")]
## doEvaluate harm monitoring bounds. Just based on infected ppts, so
## only uses 'E.j'
harmRes <- do_harm_monitoring(E.j, bounds = harmBounds.j )
## if there is indication of "harm" then need to create store output
## object then move to next trial
if ( harmRes$isHarm ) {
if ( !is.null( effCohort$cohortInd ) ) {
datI.j <- datI.j[ datI.j[[ effCohort$cohortInd ]] == 1, ]
}
fst <- finalStage1Test(
datI.j,
analysisType = "stopTime",
lowerVE = altVE,
alphaLevel = alphaAltVE,
estimand = effCohort$estimand,
lagTime = effCohort$lagTime,
time = harmRes$stopTime,
randFraction = null.p )
fest <- fst$stage1summ
VEcir <- !is.null(fest$VE_CIR)
VE_CI <- if (VEcir) {
c(fest$VE_CIR_loCI, fest$VE_CIR_upCI)
} else {
c(fest$VE_Cox_loCI, fest$VE_Cox_upCI)
}
out[[i]][[j]] <-
list(
boundHit = "Harm",
stopTime = harmRes$stopTime,
stopInfectCnt = harmRes$stopInfectCnt,
stopInfectSplit = harmRes$stopInfectSplit,
stopVE = ifelse(VEcir, fest$VE_CIR, fest$VE_Cox),
stopVE_CI= VE_CI,
VE_estimand = effCohort$estimand,
stage1complete = FALSE,
stage2complete = FALSE,
altDetected = fst$altDetected )
next
}
## If harm bounds not hit, move to nonEff monitoring
## 2. begin non-eff monitoring prep
## --------------------------------
## I don't think I need to skip this situation. The later code that deals with
## test times that are not reached should be used to handle it
#if (N1.ne > nInfec_nonEff ) {
## deal with annoying case where the trial didn't have enough infec.s
## to reach the start point of the non-eff. interim analyses
## For now we do nothing, except skip all the interim analysis code
## and go directly to the finalStage1 test
# nonEffTimes.ij <- NA
## create an empty object futRes needed later to make code work
# futRes <- list()
# effRes <- list()
## Set 'N1' to NA ?? Need to think about it
## N1 <- NA
#} else { # This 'else' goes on FOREVER!
## ------- run non-eff ---------
## We need to loop over the defined cohorts, then concatenate results across them
## The details on the non-eff cohorts exist within sublists named 'cohort1',
## 'cohort2'...
#n.neCohorts <- sum( substr(names(nonEffCohorts),1,6) == "cohort" )
cohortFutRes <- vector( "list", n.neCohorts )
for(coh.idx in 1:n.neCohorts) {
## name of component
coh.i <- paste0("cohort", coh.idx)
ind <- nonEffCohorts[[ coh.i ]]$cohortInd
## if cohortInd is specified, use it to subset the data
if (!is.null( ind ) ) {
cohI.j <- datI.j[ datI.j[[ind]] == 1, ]
} else {
cohI.j <- datI.j
}
cohortFutRes[[ coh.idx ]] <-
applyStopRules (
cohI.j,
testTimes = nonEffTimes.ij,
boundType = "nonEff",
boundLabel = "NonEffInterim",
lowerVE = lowerVEnoneff,
upperVE = nonEffCohorts[[ coh.i ]]$nullVE,
alphaLevel = alphaNoneff.ij[[coh.idx]],
lagTime = nonEffCohorts[[ coh.i ]]$lagTime,
estimand = nonEffCohorts[[ coh.i ]]$estimand,
returnAll = ifelse(n.neCohorts>1, TRUE, FALSE),
randFraction = null.p )
}
if (n.neCohorts > 1) {
## extract the 'summObj' components
summObj <- lapply(cohortFutRes, function(y) y$summObj)
names(summObj) <- paste0("cohort", 1:n.neCohorts)
## if all cohorts report a bound hit, then we'll need to dig in deeper
## to see if they all hit at the same timepoint ever.
## If they didn't all hit, then we're done here!
if ( all( sapply(cohortFutRes, function(y) y$boundWasHit) ) ) {
# extract the stopTimes from each cohort
stopIndxList <- lapply( cohortFutRes, function(x) x$stopIndx )
## find common values in all components
commonStopIndx <- unlist( Reduce(intersect, stopIndxList) )
if ( length(commonStopIndx) > 0 ) {
stopIndx <- min( commonStopIndx )
## time of test corresponding to stopIndx
stopTime <- nonEffTimes.ij[ stopIndx ]
## subset each set of results in summObj to rows 1:stopIndx
summObj <- lapply(summObj, function(y) y[1:stopIndx,] )
boundWasHit <- TRUE
boundHit <- "NonEffInterim"
} else {
boundWasHit <- FALSE
boundHit <- NA
stopTime <- NA
}
} else {
boundWasHit <- FALSE
boundHit <- NA
stopTime <- NA
}
futRes <- list( boundWasHit = boundWasHit,
boundHit = boundHit,
boundType = "NonEffInterim",
stopTime = stopTime,
summObj = summObj)
} else {
## just one cohort, so set futRes to it
futRes <- cohortFutRes[[1]][ c("boundWasHit", "boundHit",
"boundType", "stopTime", "summObj") ]
}
## if noneff hit, then store results and go to next arm
## (note: output will be a list because 'futRes' is a list
if ( futRes$boundWasHit ) {
if ( !is.null( effCohort$cohortInd ) ) {
datI.j <- datI.j[ datI.j[[ effCohort$cohortInd ]] == 1, ]
}
fst <- finalStage1Test(
datI.j,
analysisType = "stopTime",
lowerVE = altVE,
alphaLevel = alphaAltVE,
estimand = effCohort$estimand,
lagTime = effCohort$lagTime,
time = futRes$stopTime,
randFraction = null.p )
fest <- fst$stage1summ
VEcir <- !is.null(fest$VE_CIR)
VE_CI <- if (VEcir) {
c(fest$VE_CIR_loCI, fest$VE_CIR_upCI)
} else {
c(fest$VE_Cox_loCI, fest$VE_Cox_upCI)
}
out[[i]][[j]] <-
c( futRes,
list( stopVE = ifelse(VEcir, fest$VE_CIR, fest$VE_Cox),
stopVE_CI= VE_CI,
VE_estimand = effCohort$estimand,
firstNonEffCnt = N1,
stage1complete = FALSE,
stage2complete = FALSE,
altDetected = fst$altDetected ) )
next
}
## Insert efficacy analyses here. We'll need to rework later to
## allow it to be executed earlier, so parts of it can be included
## if some test are done prior to stopping for harm or non-eff
## There's no chance that we would have stopped for efficacy
## but instead ended up stopping for harm or non-eff, the two
## results are two disparate
#effCohort = list( times=c(52,77,103), timeUnit="counts", timeLag=2,
# nullVE = 0.2, totalAlpha=0.05, totalEvents=103,
# estimand="cox", lagTime=2,
# nominalAlphas=c(0.0030, 0.0183, 0.0440) ),
## clear out variable, if it exists
if ( exists("nominalAlphas.ij") )
rm( nominalAlphas.ij )
## if
## Process to derive 'effTimes.ij', etc.
if (effCohort$timeUnit == "time") {
effTimes.ij <- sort(effCohort$times)
## check if any specified times are beyond the extent of stage 1
if ( any(effTimes.ij > endStg1) ) {
wTooLrg <- which( effTimes.ij > endStg1 )
## if just one is too large, then set that last time to 'endStg1'
if ( length(wTooLrg) == 1 ) {
effTimes.ij[wTooLrg] <- endStg1
} else if ( length(effCohort$nominalAlphas) == 1 ) {
effTimes.ij[ wTooLrg[1] ] <- endStg1
effTimes.ij <- effTimes.ij[ 1:wTooLrg[1] ]
} else {
## Issue warning but continue anyway
cat("Warning: More than one of the specified efficacy test times was ",
"not reached.\n", "The code cannot adjust the nominalAlphas properly",
" for this.\n", "A final test will be inserted at the end of stage1 ",
"follow-up and will be carried\n", "out using the nominal alpha ",
"level associated with the final scheduled test\n")
## replace the first missed test with a test at the end of stage 1.
## test using 'alphaLevelNoneff' corresponding to final analysis
effTimes.ij[ wTooLrg[1] ] <- endStg1
effTimes.ij <- effTimes.ij[ 1:wTooLrg[1] ]
nominalAlphas.ij <-
effCohort$nominalAlphas[ c( 1:(wTooLrg[1]-1),
length(effCohort$nominalAlphas) ) ]
}
}
} else {
## extract info on timingCohort
if ( is.null( effCohort$timingCohort) )
stop("Component 'timingCohort' must be specified for the efficacy analysis\n.",
"It is needed to translate counts to calendar time for analyses\n")
if ( is.null(effCohort$timingCohort$lagTime) )
stop("Component 'lagTime' is a required component of 'timingCohort' in the\n",
"efficacy testing specification. Set it's value to zero if not needed\n")
if (!is.null( effCohort$timingCohort$cohortInd ) ) {
## subset data using 'cohortInd' and then determine the times
timingI.j <- datI.j[ datI.j[[ effCohort$timingCohort$cohortInd ]] == 1, ]
effTimes.ij <- getInfectionTimes(
timingI.j, cnts=sort(effCohort$times),
lagTime= effCohort$timingCohort$lagTime )
} else {
## cohortInd not specified, just use lagTime
effTimes.ij <- getInfectionTimes(
datI.j, cnts=sort(effCohort$times),
lagTime= effCohort$timingCohort$lagTime )
}
## check if we "lost" any tests, due to the counts not being reached
if ( length(effTimes.ij) < length(effCohort$times) ) {
nLost <- length(effCohort$times) - length(effTimes.ij)
## if we only lost 1, or if we're using a fixed nominal alpha level
## then add a test at the end of stage 1 follow-up
if (nLost==1 || length(effCohort$nominalAlphas)==1 ) {
effTimes.ij <- c(effTimes.ij, endStg1)
} else {
## Issue warning but continue anyway
cat("Warning: More than one of the specified efficacy test times was",
" not reached.\n", "The code cannot adjust the nominalAlphas properly",
" for this.\n", "A final test will be inserted at the end of stage1 ",
"follow-up and will be carried\n", "out using the nominal alpha level",
" associated with the final scheduled test\n\n")
## add a test at the end stage 1 follow-up. Perform the test using
## using the 'alphaLevelNoneff' value corresponding to final analysis
effTimes.ij <- c(effTimes.ij, endStg1)
nominalAlphas.ij <-
effCohort$nominalAlphas[ c( 1:(length(effTimes.ij)-1),
length(effCohort$nominalAlphas) ) ]
}
}
}
## If we've not created a customized set of values 'nominalAlphas.ij' in the
## code above, create it now
if ( !exists("nominalAlphas.ij") )
nominalAlphas.ij <- effCohort$nominalAlphas
if ( !is.null( effCohort$cohortInd ) ) {
effdatI.j <- datI.j[ datI.j[[ effCohort$cohortInd ]] == 1, ]
} else {
effdatI.j <- datI.j
}
effRes <-
applyStopRules(
effdatI.j,
testTimes = effTimes.ij,
boundType = "eff",
boundLabel = "Eff",
lowerVE = effCohort$nullVE,
alphaLevel = nominalAlphas.ij,
estimand= effCohort$estimand,
lagTime = effCohort$lagTime,
randFraction = null.p )
## 3. if nonEff NOT hit, then move on to high-eff monitoring
## We first do *only* at the stage1 times (i.e. same times as for
## nonEff monitoring ), but using all data
highEffRes <- applyStopRules(
datIall.j,
testTimes = nonEffTimes.ij,
boundType = "highEff",
boundLabel = "HighEff", ## assumed by later function
lowerVE = highVE,
alphaLevel = alphaHigh,
estimand = effCohort$estimand,
randFraction = null.p )
if ( highEffRes$boundWasHit ) {
fst <- finalStage1Test(
effdatI.j,
analysisType = "stopTime",
lowerVE = altVE,
lagTime = effCohort$lagTime,
time = highEffRes$stopTime,
randFraction = null.p )
fest <- fst$stage1summ
VEcir <- !is.null(fest$VE_CIR)
VE_CI <- if (VEcir) {
c(fest$VE_CIR_loCI, fest$VE_CIR_upCI)
} else {
c(fest$VE_Cox_loCI, fest$VE_Cox_upCI)
}
## store output, then determine value for 'altDetected'
## (note: output will be a list because 'highEffRes' is a list
out[[i]][[j]] <-
c( highEffRes,
list( summNonEff = futRes$summObj,
stopVE = ifelse(VEcir, fest$VE_CIR, fest$VE_Cox),
stopVE_CI= VE_CI,
VE_estimand = effCohort$estimand,
firstNonEffCnt = N1,
stage1complete = FALSE,
stage2complete = FALSE,
altDetected = fst$altDetected ) )
next
}
#} ## end of (N1.ne > nInfec_nonEff) if/else clause
##-----------------------------------------------------------------------
## --- skip this for now ---
if (FALSE) {
## Reached the end of stage 1 w/o hitting any boundaries (yee-ha!)
## celebrate by doing 'end-of-stage-1' dance
fst <- finalStage1Test(
datI.j,
analysisType = "final",
stage1VE = stage1VE,
lowerVE = altVE,
alphaLevel = alphaStage1,
estimand = effCohort$estimand,
lagTime = effCohort$lagTime,
boundLabels=c("Eff", "NonEffFinal"),
randFraction = null.p )
fest <- fst$stage1summ
VEcir <- !is.null(fest$VE_CIR)
VE_CI <- if (VEcir) {
c(fest$VE_CIR_loCI, fest$VE_CIR_upCI)
} else {
c(fest$VE_Cox_loCI, fest$VE_Cox_upCI)
}
## not storing into out[[i]][[j]] yet, we'll add on if not final-non-eff
finalStg1List <-
c( fst,
list( stopVE = ifelse(VEcir, fest$VE_CIR, fest$VE_Cox),
stopVE_CI= VE_CI,
VE_estimand = effCohort$estimand,
firstNonEffCnt = N1,
summNonEff = futRes$summObj,
stage1complete = TRUE,
stage2complete = FALSE ) )
## if we hit the bound at stage 1 analysis, then we're done with this cohort
if (fst$boundHit == "NonEffFinal") {
out[[i]][[j]] <- finalStg1List
next
}
}
## --- End of skipping for now ---
## -- Temp replacement for skipped code --
finalStg1List <-
c( effRes[ c("boundWasHit", "boundHit", "boundType",
"stopIndx", "stopTime") ] ,
altDetected = effRes$boundWasHit,
list(
# stopVE = ifelse(VEcir, fest$VE_CIR, fest$VE_Cox),
# stopVE_CI= VE_CI,
# VE_estimand = effEst,
firstNonEffCnt = N1,
summNonEff = futRes$summObj,
summEff = effRes$summObj,
stage1complete = TRUE,
stage2complete = FALSE ) )
## -- End temp replacement for skipped code --
## 5. If no stage1 highEff, and we "passed" the end-of-stg1 analysis, then do
## the stage2 highEff analysis
## (a) timing of it is half way between end of stage 1 and end of trial
## compute as end-of-stg1 (max(datI.j$exit)) plus half the remaining time
## (stage2- stage1)/2 until the last person from stage1 completes stg2
stg2highEffTime <- ( max(datI.j$exit) + (stage2 - stage1)/2 )
stg2highEff <-
applyStopRules(
datIall.j,
testTimes = stg2highEffTime,
boundType = "highEff",
boundLabel = "HighEff", ## assumed by later function
lowerVE = highVE,
alphaLevel = alphaHigh,
estimand=effCohort$estimand,
randFraction = null.p )
## Note: 'altDetected' result is included in 'finalStg1List', which is why
## you will not see it added to any of the following output objects
if ( stg2highEff$boundWasHit ) {
## alter some values in present in 'finalStg1List'
finalStg1List$boundHit <- stg2highEff$boundHit
finalStg1List$boundType <- stg2highEff$boundType
finalStg1List$stopTime <- stg2highEff$stopTime
## remove components from stg2Eff already present in 'finalStg1List'
dropComps <- c("boundWasHit", "boundHit", "boundType", "stopTime")
stg2HE <- stg2highEff[ !names(stg2highEff) %in% dropComps ]
## Note: 'altDetected' result is included in 'finalStg1List'
out[[i]][[j]] <- c( finalStg1List,
#list( summNonEff = futRes$summObj ),
stg2HE )
} else {
## Stage 2 completed without hitting high-efficacy
finalStg1List$stage2complete <- TRUE
finalStg1List$stopTime <- max( datIall.j$exit)
out[[i]][[j]] <- c( finalStg1List,
list(
#summNonEff = futRes$summObj,
stg2highEff = stg2highEff$summObj ) )
}
}
}
##-----------------------------------------------------------------------
# bound stage1complete altDetected(stage1 tests - 2 of them)
# "NonEffInterim" FALSE FALSE
# "HighEff" FALSE TRUE/FALSE [could be either]
# "NonEffFinal" TRUE FALSE
# "Eff" TRUE TRUE/FALSE [could be either]
# "HighEff" TRUE (Same as for "Eff")
##-----------------------------------------------------------------------
## I left this in here, but not very useful - will excise later probably
if (verbose){
for (i in 1:nTrtArms) {
cat("Probabilities of reaching each possible conclusion:\n")
print( round( table(sapply( out, function(x) x[[i]]$boundHit ),
useNA="ifany")/nTrials, 4))
if (!is.null(altVE)){
altDetectedMatrix <-
do.call("rbind",
lapply(out, function(x) x[[i]]$altDetected ))
# this occurs if each trial stops for potential harm, and thus the value of
# 'altDetected' is NULL for each trial
if ( is.null( altDetectedMatrix ) ) {
designPower <- 0
} else {
designPower <- colSums(altDetectedMatrix, na.rm=TRUE)/nTrials
}
cat("\nUnconditional power to reject the specified null hypotheses =",
format(designPower, digits=3, nsmall=3), "\n\n")
}
}
}
## save monitoring output
if ( !is.null(saveDir) ) {
if ( is.null(saveFile) ) {
if ( is.list(dataFile) )
warning(
"The output of 'monitorTrial' will not be saved to a file.\n",
"You have not specified the argument 'saveFile' and a default\n",
"filename cannot be constructed when argument 'dataFile' is ",
"a list.\n\n", immediate.=TRUE)
## used in specifying saveFile name
neEst <- ifelse( !is.null(nonEffCohorts[[1]]$estimand),
paste0("_", nonEffCohorts[[1]]$estimand), "")
saveFile <- paste0("monitorTrial",
substr(dataFile, 9, nchar(dataFile)-6),
neEst, ".RData")
}
#save(out, file = file.path(saveDir, saveFile), compress="xz" )
save(out, file = file.path(saveDir, saveFile) )
if (verbose) {
cat("Output saved in:\n", file.path(saveDir, saveFile), "\n\n")
}
}
## it should not be an "either/or" decision whether you save or output the results
return( invisible( out ) )
}
######################################## end of 'monitorTrial'#########################################
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.