R/sequential_pval.R

In gsDesign2: Group Sequential Design with Non-Constant Effect

#' @title Sequential p-value computation for gsDesign2
#' @description \code{sequential_pval} computes a sequential p-value for a group sequential design
#' created with \code{gs_design_ahr()} using a spending function approach.
#' It is the minimum of repeated p-values computed at each analysis (Jennison and Turnbull, 2000).
#' This is particularly useful for multiplicity methods such as the graphical method for group sequential designs
#' where sequential p-values for multiple hypotheses can be used as nominal p-values to plug into a multiplicity graph.
#' A sequential p-value is described as the minimum alpha level at which a one-sided group sequential bound would
#' be rejected given interim and final observed results.
#' Mild restrictions are required on spending functions used, but these are satisfied for commonly used spending functions
#' such as the Lan-DeMets spending function approximating an O'Brien-Fleming bound or a Hwang-Shih-DeCani spending function.
#'
#' @param gs_design Group sequential design generated by \code{gs_design_ahr()}.
#' @param event Event counts; numeric vector with increasing, positive values with at most one
#' value greater than or equal to largest value in \code{gs_design$analysis$event};
#' NOTE: if NULL, planned \code{event} will be used (\code{gs_design$analysis$event})
#' @param z z-value tests corresponding to analyses in \code{info}; positive values indicate a positive finding; must have the same length as \code{info}.
#' @param ustime Spending time for upper bound at specified analyses; specify default: \code{NULL} if this is to be based on information fraction;
#' if not \code{NULL}, must have the same length as \code{info}; increasing positive values with at most 1 greater than or equal to 1.
#' @param interval Interval for search to derive p-value; Default: \code{c(1e-05, 0.9999)}. Lower end of interval must be >0 and upper end must be < 1.
#' The primary reason to not use the defaults would likely be if a test were vs a Type I error <0.0001.
#' @return Sequential p-value (single numeric one-sided p-value between 0 and 1).
#' Note that if the sequential p-value is less than the lower end of the input interval,
#' the lower of interval will be returned.
#' Similarly, if the sequential p-value is greater than the upper end of the input interval,
#' then the upper end of interval is returned.
#' @references
#' Jennison C and Turnbull BW (2000), \emph{Group Sequential
#' Methods with Applications to Clinical Trials}. Boca Raton: Chapman and Hall.
#'
#' Liu, Qing, and Keaven M. Anderson. "On adaptive extensions of group sequential trials for clinical investigations." \emph{Journal of the American Statistical Association} 103.484 (2008): 1621-1630.
#'
#' Maurer, Willi, and Frank Bretz. "Multiple testing in group sequential trials using graphical approaches." \emph{Statistics in Biopharmaceutical Research} 5.4 (2013): 311-320.
#' @details
#' Solution is found with a search using \code{uniroot}.
#' This finds the maximum alpha-level for which an efficacy bound is crossed,
#' completely ignoring any futility bound.
#' @export
#'
#' @examples
#' # Derive Group Sequential Design using gsDesign2
#' x <- gs_design_ahr(
#'   enroll_rate = define_enroll_rate(duration = c(2, 2, 2, 6),
#'                                    rate = c(2.5, 5, 7.5, 10)),
#'   fail_rate = define_fail_rate(duration = Inf, fail_rate = log(2) / 6,
#'                                hr = 0.6, dropout_rate = .001),
#'   info_frac = c(.5, .65, .8, 1),
#'   analysis_time = 30,
#'   upper = gs_spending_bound,
#'   upar = list(sf = "sfLDOF", total_spend = 0.025),
#'   lower = gs_spending_bound,
#'   lpar = list(sf = "sfHSD", total_spend = 0.1, param = 2),
#'   binding = FALSE
#' )
#'
#' # Analysis at interim analysis 2
#' sequential_pval(gs_design = x, event = c(100, 160), z = c(1.5, 2))
#'
#' # Use planned spending for interim and final analyses
#' sequential_pval(gs_design = x,
#'                 event = c(100, 160, 190),
#'                 z = c(1.5, 2, 2.5))
sequential_pval <- function(gs_design,
                            event = NULL,
                            z = NULL,
                            ustime = NULL,
                            interval = c(.00001, .9999)) {
  q_e <- gs_design$input$ratio / (1 + gs_design$input$ratio)

  # gs_design should be object returned from gs_design_ahr or gs_power_ahr with >=2 analyses
  if (gs_design$design != "ahr") stop("gs_design must be an object created with gs_design_ahr() or gs_power_ahr()")
  if (nrow(gs_design$analysis) == 1) stop("gs_design must be an GSD object with more than 1 analysis")

  # check if gs_design is 1-sided, or has non-binding futility bounds
  one_sided <- all(gs_design$bound$bound == "upper")
  nonbinding <- !gs_design$input$binding
  if (!(one_sided || nonbinding)) stop("gs_design must be one-sided or have non-binding futility bounds")

  # if event is specified, check if it is an increasing numerical vector
  if (!is.null(event)) check_increasing(event)

  # if ustime is specified, check if it is an increasing numerical vector ranging within (0, 1]
  if (!is.null(ustime)) check_info_frac(ustime)

  # if event not specified, assume planned values used
  if (is.null(event)) event <- gs_design$analysis$event

  # if ustime not specified, use information fraction from design
  if (is.null(ustime)) ustime <- event/max(gs_design$analysis$event)

  # if z not specified, stop
  if (is.null(z)) stop("z must be provided")

  # check that z is numeric vector
  if (!is.vector(z, mode = "numeric")) stop("z must be numeric vector")

  # check lengths match
  if (length(event) != length(z)) stop("event and z must have the same length")
  if (length(event) != length(ustime)) stop("event and ustime must have the same length")

  # check that ustime has at most 1 value >= 1
  if (sum(ustime >= 1) > 1) stop("No more than one value of ustime >= 1")

  # check the interval is 2 values in (0,1), exclusive of endpoints
  if (!is.vector(interval, mode = "numeric") || length(interval) != 2 || min(interval) <= 0 || max(interval) >= 1) {
    stop("interval must be 2 distinct values strictly between 0 and 1")
  }

  # Get spending function and parameters from gs_design
  # sf is a character string, so we need to get the actual function from gsDesign namespace
  sf_upper <- get_sf(gs_design$input$upar$sf)
  sf_param <- gs_design$input$upar$param

  # check upper end of p-value interval input
  probhi <- sf_upper(alpha = max(interval), t = ustime, param = sf_param)$spend
  if (length(z) > 1) probhi <- probhi - c(0, probhi[1:(length(z) - 1)])
  if (min(gsDesign::gsBound1(I = event * (1 - q_e) * q_e, theta = 0, a = rep(-20, length(z)),
                             probhi = probhi)$b - z) > 0) return(max(interval))

  # check lower end of p-value interval input
  probhi <- sf_upper(alpha = min(interval), t = ustime, param = sf_param)$spend
  if (length(z) > 1) probhi <- probhi - c(0, probhi[1:(length(z) - 1)])
  if (min(gsDesign::gsBound1(I = event * (1 - q_e) * q_e, theta = 0, a = rep(-20, length(z)),
                             probhi = probhi)$b - z) < 0) return(min(interval))

  # if answer is between interval bounds, find it with root-finding
  x <- try(uniroot(sequential_zdiff, interval = -qnorm(interval), gs_design = gs_design,
                   event = event, z = z, ustime = ustime))

  if (inherits(x, "try-error")) stop("Failed to find root for sequential p-value")

  pnorm(-x$root)
}

#' This function calculates the difference between:
#' Computed efficacy bound at a given alpha level (using the spending function)
#' Observed Z-statistic
#' @noRd
sequential_zdiff <- function(x,
                             gs_design,
                             event = (1:3) * 50,
                             z = 1:3,
                             ustime = NULL) {
  q_e <- gs_design$input$ratio / (1 + gs_design$input$ratio)
  alpha <- pnorm(-x)

  sf_upper <- get_sf(gs_design$input$upar$sf)
  sf_param <- gs_design$input$upar$param

  probhi <- sf_upper(alpha = alpha, t = ustime, param = sf_param)$spend
  if (length(z) > 1) probhi <- probhi - c(0, probhi[1:(length(z) - 1)])

  return(min(gsDesign::gsBound1(I = event * (1 - q_e) * q_e, theta = 0, a = rep(-20, length(z)),
                                probhi = probhi)$b - z))
}