R/dove2.R

In DOVE: Durability of Vaccine Efficacy

#' Durability of Vaccine Efficacy 
#' 
#' Estimates the potentially waning long-term efficacy of vaccines in
#'  randomized, placebo-controlled clinical trials with staggered
#'  enrollment of participants and sequential crossover of placebo recipients.
#'  The log hazard ratio for the vaccine effect is assumed to be a 
#'  piecewise linear function of time since vaccination, whereas the
#'  baseline hazard function is nonparametric.
#' 
#'  In dove2(), the log hazard ratio for the vaccine effect is a piecewise 
#'  linear function of time since vaccination.
#'  Its sister function, dove(), allows the hazard ratio to be a 
#'  nonparametric function of time for which confidence intervals are not 
#'  provided. The approach implemented in dove2() provides more precise and 
#'  more stable estimates of vaccine efficacy on the hazard rate and includes 
#'  proper confidence intervals.
#'
#' The information required for an analysis is 
#'   \describe{
#'     \item{Entry Time:}{Calendar time when the participant enters the 
#'       trial.}
#'     \item{Event Time:}{Calendar time when the participant experiences
#'       the clinical event of interest (e.g., symptomatic COVID-19) or the
#'       follow-up ends, whichever occurs first.}
#'     \item{Event Status:}{Binary indicator taking value 1 if the clinical
#'       event of interest occurs before the end of follow-up and 0 otherwise.}
#'     \item{Vaccination Status:}{Binary indicator taking value 1 if
#'       vaccination occurs before the end of follow-up and 0 otherwise.}
#'     \item{Vaccination Time:}{Calendar time when vaccination takes place,
#'       with NA, Inf, or an arbitrary non-negative value
#'       if the participant is not vaccinated.}
#'     \item{Covariates:}{Baseline covariates (e.g., priority group, age, 
#'       ethnicity).}
#'    }
#'
#' Note that all the time variables are measured from the start of the 
#' clinical trial and are specified in units of days. 
#' For each individual, the entry_time, event_time, and vaccination_time must satisfy
#' entry_time \eqn{\le} event_time and entry_time \eqn{\le} vaccination_time. 
#' If entry_time > event_time or entry_time > vaccination_time, the case will be 
#' removed from the analysis and a message will be generated.  
#' 
#' The general structure of the formula input is
#'   \preformatted{
#'   Surv(event_time, event_status) ~ covariates + 
#'     vaccine(entry_time, vaccination_status, vaccination_time)
#'   }
#' 
#' The response variable must be a survival analysis object as returned by the
#' 'Surv()' function of package \pkg{survival}, where event_time is the
#' observation time (formal argument 'time') and event_status is the status
#' indicator input (formal argument 'event'). Specifically, 
#' \preformatted{Surv(time = event_time, event = event_status)}
#' 
#' The covariates can be either numerical or categorical.
#' If categorical covariates are provided, all other categories are compared to the first category.
#' A model without covariates is also allowed.
#' 
#' @rdname dove2
#' @name dove2 
#' 
#' @param formula A formula object, with the response on the left-hand side of a
#'   '~' operator, and the covariates and vaccine() function on the right-hand side.  
#'   The response must be a survival analysis object as returned by the 'Surv'
#'   function of the \pkg{survival} package. See Details for further information.
#'   The vaccine() function must be used to specify the entry time and
#'   vaccination information. See ?vaccine and Details for further
#'   information. 
#'   
#' @param data A data.frame object. The data.frame in which to interpret the
#'   variable names in formula. Must contain the entry time, the event time,
#'   the event status, the vaccination status,
#'   the vaccination time, and the covariates. See Details.
#'   
#' @param plots A logical object. If TRUE (default), plots of the estimated 
#'   VE in reducing the attack rate, the estimated VE in reducing the hazard rate, 
#'   and their 95\% confidence intervals will be automatically generated. 
#'   If FALSE, plots will not be generated.
#'   
#' @param changePts A numerical vector object or NULL. The potential change
#'   points (in days) of the piece-wise log-linear hazard ratio for the
#'   vaccine effect. See Details 
#'   for further information. If NULL, one change point will automatically be 
#'   selected among \{28, 35, 42, 49, 56\} (Weeks 4, 5, 6, 7, 8) by the 
#'   Akaike information criterion (AIC).
#'   
#' @param constantVE A logical object. If FALSE (default), VE is allowed to 
#'   vary after the last change point. If TRUE, VE is assumed to be 
#'   constant after the last change point.
#' 
#' @param timePts A numerical vector object or NULL. The endpoints (in days) 
#'   of the time periods for which the VE in reducing the attack rate are to
#'   be estimated. If NULL, a default sequence 
#'   \eqn{t_1, 2t_1, 3t_1, \dots} will be used, where \eqn{t_1} is the first 
#'   change point. The sequence ends at the maximum of the event times from all 
#'   participants. This input is ignored when constantVE = TRUE.  
#'   
#' @returns An S3 object of class DOVE containing a list with elements
#'
#'   \item{call}{The unevaluated call.}
#'
#'   \item{changePts}{The changePts of the analysis.}
#'
#'   \item{covariates}{A matrix containing the estimated (log) hazard ratio of
#'     each covariate, together with the estimated standard error, the 95\%
#'     confidence interval, and the two-sided p-value for testing no covariate
#'     effect. NA if no covariate is included in the model.}
#'     
#'   \item{vaccine}{A list containing one or three elements, depending on the 
#'     value of constantVE. 
#'     If constantVE = TRUE, the only element is named 'VE' and is a vector 
#'     containing the estimate of constant VE, its standard error estimate, 
#'     and the 95\% confidence interval.
#'     If constantVE = FALSE, three matrices are returned. The first matrix 
#'     named 'VE_a' contains the daily VE estimates in reducing the attack 
#'     rate, together with the 95\% confidence intervals. The second matrix 
#'     named 'VE_h' contains the daily VE estimates in reducing the hazard rate,
#'     together with the 95\% confidence intervals.
#'     The third matrix named 'VE_period' contains the estimates of VE in
#'     reducing the attack rate over successive time periods according to
#'     timePts, together with the 95\% confidence intervals.}
#'
#' @references Lin, D-Y, Gu, Y., Zeng, D., Janes, H. E., and Gilbert, P. B. (2021). 
#'   Evaluating vaccine efficacy against SARS-CoV-2 infection. 
#'   Clinical Infectious Diseases, ciab630, https://doi.org/10.1093/cid/ciab630.
#'
#' @export
#' @import methods
#'
#' @useDynLib DOVE
#'
#' @include verifyInputs.R CoxReg.R 
#' 
#' @examples
#' data(doveData)
#'
#' set.seed(1234)
#' smp <- sample(1L:nrow(x = doveData), size = 2500L)
#' 
#' # NOTE: This sample size is chosen for example only -- larger data sets
#' # should be used.
#' # See the vignette for a full analysis of the doveData dataset
#'
#' # Fit the model with default settings
#' dove2(formula = Surv(event.time, event.status) ~ priority + sex + 
#'                 vaccine(entry.time, vaccine.status, vaccine.time), 
#'       data = doveData[smp,])
#' 
#' # Specify Week 4 as the change point
#' # Assume a potentially waning VE after 4 weeks
#' # Estimate VE_a over 0-4, 4-16, 16-28, 28-40 weeks
#' dove2(formula = Surv(event.time, event.status) ~ priority + sex + 
#'                 vaccine(entry.time, vaccine.status, vaccine.time), 
#'       data = doveData[smp,],
#'       changePts = 4*7,
#'       timePts = c(4, 16, 28, 40)*7)
#'       
#' # Specify multiple change points at Weeks 4 and 8
#' # Assume a constant VE after 8 weeks
#' dove2(formula = Surv(event.time, event.status) ~ priority + sex + 
#'                 vaccine(entry.time, vaccine.status, vaccine.time), 
#'       data = doveData[smp,],
#'       changePts = c(4, 8)*7,
#'       constantVE = TRUE)
# 6/30/2021 introduced in DOVE v1.7
dove2 <- function(formula,  
                  data, 
                  plots = TRUE, 
                  changePts = NULL,
                  constantVE = FALSE,
                  timePts = NULL) {

  # matched call to include with returned object
  cl <- match.call()
  
  # formula and data must be provided

  if (missing(x = formula)) {
    stop("a formula argument must be provided", call. = FALSE)
  }
  
  if (missing(x = data)) {
    stop("a data argument must be provided", call. = FALSE)
  }
  
  # process inputs to obtain time and covariate information and to
  # verify timePts
  inputs <- .verifyInputs(formula = formula, data = data, timePts = timePts)

  # notify user if default timePts will be used
  if (is.null(x = inputs$timePts) && !constantVE) {
    message("timePts not given; default values will be used")
  }

  # verify changePts
  if (!is.null(x = changePts)) {

    if (!is.numeric(x = changePts)) {
      stop("changePts must be a numeric vector or NULL", call. = FALSE)
    }
    
    if (length(x = changePts) == 0L) {
      message("changePts has zero length; default values will be used")
      changePts <- NULL
    } else if (any(changePts <= 0.0)) {
      stop("all changePts must be positive", call. = FALSE)
    } else if (any(changePts >= inputs$tau)) {
      message("changePts > tau have been removed")
      changePts <- changePts[changePts < inputs$tau]
      if (length(x = changePts) == 0L) {
        message("changePts has zero length; default values will be used")
        changePts <- NULL
      }
    }
    
  }

  # knots (in days) of linear splines 
  if (is.null(x = changePts)) {
    def <- (4:8)*7
    use <- def < inputs$tau
    knots <- as.list(x = def[use])
    message("changePts not given; using AIC to select from {",
            paste(def[use],collapse=","), "}")
  } else {
    knots <- list(changePts)
  }
    
  ### main analysis
  
  res <- CoxReg(dt = inputs$data,  
                X = inputs$X,  
                knots = knots,  
                constantVE = constantVE,  
                timePts = inputs$timePts,  
                plots = plots,
                tau = inputs$tau)

  res[["call"]] <- cl

  class(x = res) <- "DOVE"
  attr(x = res, which = "type") <- 2L + constantVE

  return( res )
}