Issue 89: Integer designs, vaccine efficacy designs

DESCRIPTION

@@ -1,5 +1,5 @@
 Package: gsDesign
-Version: 3.4.0.9000
+Version: 3.4.0.9001
 Title: Group Sequential Design
 Authors@R: person(given = "Keaven", family = "Anderson", email =
           "keaven_anderson@merck.com", role = c('aut','cre'))

NAMESPACE

@@ -85,6 +85,8 @@ export(spendingFunction)
 export(ssrCP)
 export(tEventsIA)
 export(testBinomial)
+export(toBinomialExact)
+export(toInteger)
 export(varBinomial)
 export(xprint)
 export(xtable)

NEWS.md

@@ -1,6 +1,12 @@
-# gsDesign 3.4.1 (July, 2023?)
+# gsDesign 3.5.0 (July, 2023)
 
-- Allowed wider parameter range (0,15] for sfPower 
+- Allowed wider parameter range (0, 15] for `sfPower()`.
+- Added function `toInteger()` to convert `gsDesign` or `gsSurv` classes
+  to integer sample size and event counts.
+- Added function `toBinomialExact()` to convert time-to-event bounds to
+  exact binomial for low event rate studies.
+- Added "A Gentle Introduction to Group Sequential Design" vignette for
+  an introduction to asymptotics for group sequential design.
 
 # gsDesign 3.4.0 (October, 2022)
 
@@ -16,7 +22,7 @@
   - Vaccine efficacy design using spending bounds and exact binomial boundary crossing probabilities
 - Minor fix to labeling in print.gsProbability
 - Fixed error in sfStep
-- Updates to reduce R CMD check and other minor issues 
+- Updates to reduce R CMD check and other minor issues
 
 # gsDesign 3.2.2 (January, 2022)
 
@@ -37,8 +43,8 @@
 - Substantially updated unit testing to increase code coverage above 80%
 - Updated error checking messages to print function where check fails
 - Removed dependencies on plyr packages
-- Updated github actions 
- 
+- Updated github actions
+
 # gsDesign 3.1.1 (May, 2020)
 
 - Vignettes updated
@@ -69,7 +75,7 @@
 
 - First Github-based release
 - Cleaned up documentation for `nBinomial1Sample()`
-- Updated documentation and code (including one default value for an argument) for `nBinomial1Sample()` to improve error handling and clarity 
+- Updated documentation and code (including one default value for an argument) for `nBinomial1Sample()` to improve error handling and clarity
 - Updated `sfLDOF()` to generalize with rho parameter; still backwards compatible for Lan-DeMets O'Brien-Fleming
 
 # gsDesign 3.0-3
@@ -90,7 +96,7 @@
 
 # gsDesign 3.0-0 (December, 2015)
 
-- Updated xtable extension to meet R standards for extensions. 
+- Updated xtable extension to meet R standards for extensions.
 - Fixed `xtable.gsSurv` and `print.gsSurv` to work with 1-sided designs
 - Update to calls to ggplot to replace show_guide (deprecated) with `show.legend` arguments where used in `ggplot2::geom_text` calls; no user impact
 - Minor typo fixed in `sfLogistic` help file

R/gsSurv.R

@@ -903,12 +903,11 @@ KT <- function(alpha = .025, sided = 1, beta = .1,
 #' follow-up duration, or "Power" if accrual rates and duration as well as
 #' follow-up duration was specified and \code{beta=NULL} was input.}
 #'
-#' \code{gsSurv()} returns much of the above plus an object of class
-#' \code{gsDesign} in a variable named \code{gs}; see \code{\link{gsDesign}}
+#' \code{gsSurv()} returns much of the above plus variables in the class
+#' \code{gsDesign}; see \code{\link{gsDesign}}
 #' for general documentation on what is returned in \code{gs}.  The value of
 #' \code{gs$n.I} represents the number of endpoints required at each analysis
 #' to adequately power the trial. Other items returned by \code{gsSurv()} are:
-#' \item{gs}{A group sequential design (\code{gsDesign}) object.}
 #' \item{lambdaC}{As input.} \item{etaC}{As input.} \item{etaE}{As input.}
 #' \item{gamma}{As input unless none of the following are \code{NULL}:
 #' \code{T}, \code{minfup}, \code{beta}; otherwise, this is a constant times
@@ -926,10 +925,7 @@ KT <- function(alpha = .025, sided = 1, beta = .1,
 #' extension for the \code{nSurv} class.} \item{eDC}{A vector of expected
 #' number of events by stratum in the control group under the alternate
 #' hypothesis.} \item{eDE}{A vector of expected number of events by stratum in
-#' the experimental group under the alternate hypothesis.} \item{eDC0}{A vector
-#' of expected number of events by stratum in the control group under the null
-#' hypothesis.} \item{eDE0}{A vector of expected number of events by stratum in
-#' the experimental group under the null hypothesis.} \item{eNC}{A vector of
+#' the experimental group under the alternate hypothesis.} \item{eNC}{A vector of
 #' the expected accrual in each stratum in the control group.} \item{eNE}{A
 #' vector of the expected accrual in each stratum in the experimental group.}
 #' \item{variable}{A text string equal to "Accrual rate" if a design was
@@ -949,7 +945,7 @@ KT <- function(alpha = .025, sided = 1, beta = .1,
 #' enrollment in the control group at the output time \code{T}.} \item{eNE}{The
 #' expected enrollment in the experimental group at the output time \code{T}.}
 #'
-#' \code{tEventsIA()} returns
+#' \code{tEventsIA()} returns the same structure as \code{nEventsIA(..., simple=TRUE)} when
 #' @author Keaven Anderson \email{keaven_anderson@@merck.com}
 #' @seealso \code{\link{gsBoundSummary}}, \code{\link{xprint}},
 #' \link{gsDesign package overview}, \link{plot.gsDesign},
@@ -1158,7 +1154,7 @@ nEventsIA <- function(tIA = 5, x = NULL, target = 0, simple = TRUE) {
   if (simple) {
     if (class(x)[1] == "gsSize") {
       d <- x$d
-    } # OLD      else d <- sum(x$eDC[length(x$eDC)]+x$eDE[length(x$eDE)])
+    }
     else if (!is.matrix(x$eDC)) {
       d <- sum(x$eDC[length(x$eDC)] + x$eDE[length(x$eDE)])
     } else {
@@ -1198,6 +1194,8 @@ gsSurv <- function(k = 3, test.type = 4, alpha = 0.025, sided = 1,
   z <- gsnSurv(x, y$n.I[k])
   eDC <- NULL
   eDE <- NULL
+  eDC0 <- NULL
+  eDE0 <- NULL
   eNC <- NULL
   eNE <- NULL
   T <- NULL
@@ -1206,12 +1204,18 @@ gsSurv <- function(k = 3, test.type = 4, alpha = 0.025, sided = 1,
     T <- c(T, xx$T)
     eDC <- rbind(eDC, xx$eDC)
     eDE <- rbind(eDE, xx$eDE)
+    # Following is a placeholder
+    # eDC0 <- rbind(eDC0, xx$eDC0) # Added 6/15/2023
+    # eDE0 <- rbind(eDE0, xx$eDE0) # Added 6/15/2023
     eNC <- rbind(eNC, xx$eNC)
     eNE <- rbind(eNE, xx$eNE)
   }
   y$T <- c(T, z$T)
   y$eDC <- rbind(eDC, z$eDC)
   y$eDE <- rbind(eDE, z$eDE)
+  # Following is a placeholder
+  # y$eDC0 <- rbind(eDC0, z$eDC0) # Added 6/15/2023
+  # y$eDE0 <- rbind(eDE0, z$eDE0) # Added 6/15/2023
   y$eNC <- rbind(eNC, z$eNC)
   y$eNE <- rbind(eNE, z$eNE)
   y$hr <- hr

R/toBinomialExact.R

@@ -0,0 +1,128 @@
+#' 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.
+#'
+#' @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)
+toBinomialExact <- function(x) {
+  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
+  k <- xx$k
+  counts <- xx$n.I
+
+  # 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
+  atem <- a
+  alpha_spend <- x$upper$sf(alpha = x$alpha, t = xx$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, lower.tail = TRUE)
+    btem <- b
+    # Compute target beta-spending
+    beta_spend <- xx$lower$sf(alpha = xx$beta, t = xx$timing, param = xx$lower$param)$spend
+  } else {
+    b <- counts + 1 # test.type = 1 means no futility bound
+    nbupperprob <- 0
+  }
+  for (j in 1:x$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 = atem[1:max(j, 2)], b = counts[1:max(j, 2)] + 1
+    )$lower$prob[1:j])
+    atem <- a # Work space for updating efficacy bound, if needed
+    # 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]
+        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 && j < xx$k) {
+      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] == a[j] + 1) break # Cannot make a and b bounds the same
+          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] <- 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 = b[1:max(j, 2)]
+          )$upper$prob[1:j])
+        }
+      }
+    }
+  }
+  b[xx$k] <- a[xx$k] + 1 # Final upper bound = lower bound + 1
+  xxxx <- gsBinomialExact(k = k, theta = c(p0, p1), n.I = counts, a = a, b = b)
+  return(xxxx)
+}