#' 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)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.