R/toBinomialExact.R

In gsDesign: Group Sequential Design

#' Translate survival design bounds to exact binomial bounds
#'
#' @param x An object of class \code{gsSurv}; i.e., an object generated by
#'   the \code{gsSurv()} function.
#' @param observedEvents If NULL (default), targeted timing of analyses will come from \code{x$n.I}. 
#' Otherwise, this should be vector of increasing positive integers with at most 1 value \code{>= x$n.IPlan} and of length at least 2.
#' Only one value can be greater than or equal to \code{x$maxn.IPlan}. 
#' This determines the case count at each analysis performed. 
#' Primarily, this is used for updating a design at the time of analysis.
#' 
#' @details
#' The exact binomial routine \code{gsBinomialExact} has requirements that may not be satisfied
#' by the initial asymptotic approximation. 
#' Thus, the approximations are updated to satisfy the following requirements of \code{gsBinomialExact}:
#' \code{a} (the efficacy bound) must be positive, non-decreasing, and strictly less than n.I
#' \code{b} (the futility bound) must be positive, non-decreasing, strictly greater than a
#' \code{n.I - b} must be non-decreasing and >= 0
#' 
#' @return An object of class \code{gsBinomialExact}.
#'
#' @seealso \code{\link{gsBinomialExact}}
#'
#' @export
#'
#' @examples
#' # The following code derives the group sequential design using the method
#' # of Lachin and Foulkes
#'
#' x <- gsSurv(
#'   k = 3,                 # 3 analyses
#'   test.type = 4,         # Non-binding futility bound 1 (no futility bound) and 4 are allowable
#'   alpha = .025,          # 1-sided Type I error
#'   beta = .1,             # Type II error (1 - power)
#'   timing = c(0.45, 0.7), # Proportion of final planned events at interims
#'   sfu = sfHSD,           # Efficacy spending function
#'   sfupar = -4,           # Parameter for efficacy spending function
#'   sfl = sfLDOF,          # Futility spending function; not needed for test.type = 1
#'   sflpar = 0,            # Parameter for futility spending function
#'   lambdaC = .001,        # Exponential failure rate
#'   hr = 0.3,              # Assumed proportional hazard ratio (1 - vaccine efficacy = 1 - VE)
#'   hr0 = 0.7,             # Null hypothesis VE
#'   eta = 5e-04,           # Exponential dropout rate
#'   gamma = 10,            # Piecewise exponential enrollment rates
#'   R = 16,                # Time period durations for enrollment rates in gamma
#'   T = 24,                # Planned trial duration
#'   minfup = 8,            # Planned minimum follow-up
#'   ratio = 3              # Randomization ratio (experimental:control)
#' )
#' # Convert bounds to exact binomial bounds
#' toBinomialExact(x)
#' # Update bounds at time of analysis
#' toBinomialExact(x, observedEvents = c(20,55,80))
toBinomialExact <- function(x, observedEvents = NULL) {
  if (!inherits(x, "gsSurv")) stop("toBinomialExact must have class gsSurv as input")
  if (x$test.type != 1 && x$test.type != 4) stop("toBinomialExact input test.type must be 1 or 4")
  # Round interim sample size (or events for gsSurv object)
  xx <- if (max(round(x$n.I) != x$n.I)) toInteger(x) else x
  if(is.null(observedEvents)){
    counts <- xx$n.I
    k <- xx$k
  }else{
    if (!isInteger(observedEvents)) stop("toBinomialExact: observedEvents must be a vector of increasing positive integers")
    if (min(observedEvents - dplyr::lag(observedEvents, default = 0)) < 1) 
      stop("toBinomialExact: observedEvents must be a vector of increasing positive integers")
    counts <- observedEvents
    if (sum(observedEvents >= xx$maxn.IPlan) > 1) stop("toBinomialExact: at most 1 value in observedEvents can be >= maximum planned (x$maxn.IPlan)")
    k <- length(observedEvents)
    if (k < 2) stop("toBinomialExact: must have at least 2 values in observedEvents")
    xx <- gsDesign(
      k = k,
      test.type = x$test.type,
      alpha = x$alpha,
      beta = x$beta,
      astar = x$astar,
      sfu = x$upper$sf,
      sfupar = x$upper$param,
      sfl = x$lower$sf,
      sflpar = x$lower$param,
      n.I = observedEvents,
      maxn.IPlan = xx$n.I[x$k],
      delta = x$delta,
      delta1 = x$delta1,
      delta0 = x$delta0
    )
    xx$hr0 <- x$hr0
  }
  # Translate vaccine efficacy to exact binomial probabilities

  p0 <- x$hr0 * x$ratio / (1 + x$hr0 * x$ratio)
  p1 <- x$hr * x$ratio / (1 + x$hr * x$ratio)
  # Lower bound probabilities are for efficacy and Type I error should be controlled under p0
  a <- qbinom(p = pnorm(-xx$upper$bound), size = counts, prob = p0) - 1
  init_approx <- list(a = a) # save initial efficacy bound approximation

  # check that a is non-decreasing, >= -1, and < n.I
  a <- pmax(a, -1)
  a <- pmax(a, dplyr::lag(a, def = -1))
  a <- pmin(a, counts - 1)

  atem <- a
  timing <- counts / xx$maxn.IPlan
  alpha_spend <- x$upper$sf(alpha = x$alpha, t = timing, param = x$upper$param)$spend
  if (x$test.type != 1) {
    # Upper bound probabilities are for futility
    # Compute nominal p-values under H0 for futility and corresponding inverse binomial under H1
    b <- qbinom(p = pnorm(-xx$lower$bound), size = counts, prob = p0)
    init_approx$b <- b # save initial futility bound approximation

    # check that b is non-decreasing, > a, and n.I - b is non-decreasing
    b <- pmin(b, counts + 1)
    b <- pmax(a + 1, b)
    b <- pmin(b, counts - dplyr::lag(counts, def = 0) + dplyr::lag(b, def = 1))
    # Compute target beta-spending
    beta_spend <- xx$lower$sf(alpha = xx$beta, t = timing, param = xx$lower$param)$spend
  } else {
    b <- counts + 1 # test.type = 1 means no futility bound
    nbupperprob <- 0
  }
  for (j in 1:k) {
    # Non-binding bound assumed.
    # Compute spending through analysis j.
    # Upper bound set to > counts so that it cannot be crossed;
    # this is to compute lower bound spending with non-binding futility bound.
    # NOTE: cannot call gsBinomialExact with k == 1, so make it at least 2
    # cumulative spending through analysis j
    nblowerprob <- sum(gsBinomialExact(
      k = max(j, 2), theta = p0, n.I = counts[1:max(j, 2)],
      a = a[1:max(j, 2)], b = counts[1:max(j, 2)] + 1
    )$lower$prob[1:j])
    atem <- a # Work space for updating efficacy bound
    btem <- b # Work space for updating futility bound
    # Set range for possible changes to a[j]
    amin <- ifelse(j > 1, a[j - 1], -1)
    amax <- ifelse(j < k, counts[j] - 1, counts[j])
    a[j] <- ifelse(a[j] > amax, amax, a[j])
    a[j] <- ifelse(a[j] < amin, amin, a[j])
    if (amin > amax) stop(paste("toBinomialExact: amin > amax: amin =", amin, "amax =", amax, "j =", j))
    atem[j] <- a[j]
    # If less than allowed spending, check if bound can be increased
    if (nblowerprob < alpha_spend[j]) {
      while (nblowerprob < alpha_spend[j]) {
        a[j] <- atem[j]
        if (a[j] >= amax - 1) break # keep in allowable range 
        atem[j] <- atem[j] + 1
        nblowerprob <- sum(gsBinomialExact(
          k = max(j, 2), theta = p0, n.I = counts[1:max(j, 2)],
          a = atem[1:max(j, 2)], b = counts[1:max(j, 2)] + 1
        )$lower$prob[1:j])
      }
      # If > allowed spending, reduce bound appropriately
    } else if (nblowerprob > alpha_spend[j]) {
      while (nblowerprob > alpha_spend[j]) {
        a[j] <- a[j] - 1
        nblowerprob <- sum(gsBinomialExact(
          k = max(j, 2), theta = p0, n.I = counts[1:max(j, 2)],
          a = a[1:max(j, 2)], b = counts[1:max(j, 2)] + 1
        )$lower$prob[1:j])
      }
    }
    # beta-spending, if needed
    if (x$test.type == 4)
      # Set range for possible values of b[j]
      bmin <- a[j] + 1 # must be strictly > a[j]
      bmin <- ifelse(j == 1, bmin, max(bmin, b[j - 1])) # must be non-decreasing 
      bmax <- counts[j] + 1
      bmax <- ifelse(j == 1, bmax, min(bmax, counts[j] - counts[j - 1] + b[j - 1]))
      if (bmin > bmax) stop(paste("bmin > bmax: bmin =", bmin, "bmax =", bmax, "j =", j))
      b[j] <- ifelse(b[j] > bmax, bmax, b[j])
      b[j] <- ifelse(b[j] < bmin, bmin, b[j])
      btem[j] <- b[j]
      upperprob <- sum(gsBinomialExact(
        k = max(j, 2), theta = p1, n.I = counts[1:max(j, 2)],
        a = a[1:max(j, 2)], b = b[1:max(j, 2)]
      )$upper$prob[1:j])
      if (upperprob < beta_spend[j]) {
        while (upperprob < beta_spend[j]) {
          b[j] <- btem[j]
          if (btem[j] == bmin) break # only lower if range allows
          btem[j] <- btem[j] - 1
          upperprob <- sum(gsBinomialExact(
            k = max(j, 2), theta = p1, n.I = counts[1:max(j, 2)],
            a = a[1:max(j, 2)], b = btem[1:max(j, 2)]
          )$upper$prob[1:j])
        }
        
      } else if (upperprob > beta_spend[j]) {
        while (upperprob > beta_spend[j] && 
               b[j] < bmax) {
          b[j] <- b[j] + 1
          upperprob <- sum(gsBinomialExact(
            k = max(j, 2), theta = p1, n.I = counts[1:max(j, 2)],
            a = a[1:max(j, 2)], b = b[1:max(j, 2)]
          )$upper$prob[1:j])
        }
      }
  }
  xxxx <- gsBinomialExact(k = k, theta = c(p0, p1), n.I = counts, a = a, b = b)
  xxxx$init_approx <- init_approx
  return(xxxx)
}