R/powlgrnk6.R

In raredd/desmon: Functions for design and monitoring of clinical trials

#' Calculate the Power of a Group Sequential Logrank Test
#'
#' Performs a simulation to calculate the power of a group sequential logrank
#' test using an error spending rate function upper boundary with repeated
#' confidence intervals for early stopping in favor of the null hypothesis.
#' Assumes exponential distributions for the failure times.
#'
#' @details
#' Performs a simulation to estimate the rejection probability under the
#' specified failure rates of a group sequential logrank test. \code{2*Npg}
#' subjects will be accrued uniformly over \code{acc.per} time units, and
#' followed until \code{total.inf} failures are available. The first interim
#' analysis is performed when at least \code{total.inf*min.inf} failures are
#' available.  Following that, interim analyses are performed every
#' \code{anal.int} units of time. The final analysis is performed at
#' \code{total.inf} failures, even if that does not occur at an otherwise
#' scheduled analysis. It is extremely important that the time units used in
#' the values of \code{control.rate}, \code{test.rate}, \code{acc.per}, and
#' \code{anal.int} be the same (eg all in years and number events per year).
#' The upper boundary for rejecting the null hypothesis is obtained from the
#' \code{\link{sequse}} function, using the boundary of type \code{use}.  The
#' study will be stopped early in favor of the null hypothesis if the upper RCI
#' on the hazard ratio (control/experimental) is \code{ < ratio} (that is, if
#' the target difference falls outside the RCI).  The confidence level and the
#' shape of the boundary of the RCI can be different in the RCI than in the
#' upper boundary.
#'
#' The truncated boundary with \code{use=6} is different (and probably
#' preferred) to the default settings with \code{use=1} and
#' \code{trunc=alpha/50}.  The latter takes the standard O'Brien-Fleming
#' boundary and truncates it at early analyses with no other adjustment, giving
#' a test of slightly inflated size.  The amount of the inflation is generally
#' negligible, though (see Freidlin, Korn and George, 1999). Using \code{use=6}
#' truncates the boundary at \code{alpha/50} (do not specify a larger value of
#' trunc), but adjusts the rest of the boundary to give the correct size.
#'
#' Failure distributions other than the exponential can be used by providing a
#' function to generate samples from the desired distribution via the
#' \code{gft} argument.  This function will be called as
#' \code{gft(rx,control.rate,test.rate,...)}, where \code{rx} is a vector of
#' length \code{2*Npg}, coded 0 for the control group and 1 for the
#' experimental treatment group, \code{control.rate} and \code{test.rate} are
#' the values given in the call to \code{powlgrnk6}, and \code{...} can be
#' used to pass additional arguments to \code{gft()}.  While any distributions
#' can be used, it should be noted that the logrank test and RCI methodology
#' are probably only appropriate when the study is targeting a proportional
#' hazards alternative.
#'
#' @param Npg Planned number of subjects in the each treatment group (control
#' and experimental)
#' @param control.rate Exponential failure rate for the control group
#' @param test.rate Exponential failure rate for the experimental treatment
#' group
#' @param acc.per Length of time over which the \code{2*Npg} subjects will be
#' accrued
#' @param total.inf Total number of failures at full information
#' @param anal.int The length of time between planned interim analyses
#' @param min.inf The minimum information that must be available at the first
#' interim analysis
#' @param nsamp number of samples to generate in the simulation
#' @param alpha One-sided significance level for rejecting the null hypothesis,
#' expressed as a proportion (rather than a percent)
#' @param conf The two-sided confidence level (proportion) for the repeated
#' confidence interval used for early stopping in favor of the null. If
#' \code{conf <=0} or \code{>=1} then the RCI is not used.
#' @param ratio The study will be stopped in favor of the null when the upper
#' RCI on the hazard ratio (control/experimental) \code{< ratio}
#' @param use The type of boundary used in \code{\link{sequse}} for the upper
#' boundary
#' @param use.rci The type of boundary used in \code{\link{sequse}} for the RCI
#' @param trunc If 0<\code{trunc}<1, then the upper boundary is truncated at
#' the significance level given by \code{trunc}; that is, if the nominal
#' significance level of the upper boundary is < \code{trunc}, the critical
#' value corresponding to \code{trunc} is used instead.  The default is .0005
#' when \code{alpha}=.025.
#' @param gft Function to generate a sample of failure times.  Default uses
#' exponential distributions with rates \code{control.rate} for the control
#' group and \code{test.rate} for the experimental treatment group
#' @param ... Additional arguments to \code{gft()}
#'
#' @return
#' Returns a \code{5 x nsamp} matrix giving the results for each sample in
#' the corresponding column, of class `powlgrnk6'.  The first row gives a
#' code indicating the type of stopping that occurred (1=crossed upper
#' boundary, 2=stopped in favor of the null, 3=reached full information without
#' crossing either boundary), the second gives the logrank test statistic at
#' the analysis where the study stopped (on the standard normal scale), the
#' third gives the information time when the study stopped, the 4th the number
#' of time units from the start of the study until the study was stopped, and
#' the 5th gives the total number of interim analyses performed (including at
#' full information, if reached).
#'
#' @seealso
#' \code{\link{sequse}}; \code{\link{rci}}; \code{\link{print.powlgrnk6}};
#' \code{\link{seqopr}}; \code{\link{seqss}}
#'
#' @references
#' Freidlin, Korn and George, 1999. \emph{Controlled Clinical  Trials}
#' \strong{20}:395-407.
#'
#' Jennison and Turnbull, 1990, Statistical Science 5:299-317.
#'
#' @keywords design survival
#'
#' @examples
#' ## power under the alternative; 7 years accrual, analyses every 1/2 year
#' ## one-sided 0.025 test, with an 80% RCI used for early stopping in
#' ## favor of H0
#' ## use larger nsamp (examples for illustration only)
#' out <- powlgrnk6(210, 0.11889, 0.11889 / 1.5, 7, 233, 0.5, nsamp = 10, conf = 0.80)
#' table(out[1, ]) / ncol(out)
#' apply(out[3:5, ], 1, mean)
#' out
#'
#' ## size under H0; note specification of ratio to keep the stopping rule
#' ## the same as under the alternative
#' out <- powlgrnk6(210, 0.11889, 0.11889, 7, 233, 0.5, nsamp = 10,
#'                  conf = 0.80, ratio = 1.5)
#' table(out[1, ]) / ncol(out)
#' apply(out[3:5, ], 1, mean)
#'
#' ## power without RCI
#' out <- powlgrnk6(210, 0.11889, 0.11889 / 1.5, 7, 233, 0.5, nsamp = 10, conf = 1)
#' table(out[1, ]) / ncol(out)
#' apply(out[3:5, ], 1, mean)
#'
#' ## size without truncation
#' out <- powlgrnk6(210, 0.11889, 0.11889, 7, 233, 0.5, nsamp = 10,
#'                  conf = 0.80, ratio = 1.5, trunc=0)
#' table(out[1, ]) / ncol(out)
#' apply(out[3:5, ], 1, mean)
#'
#' @export

powlgrnk6 <- function(Npg, control.rate, test.rate, acc.per, total.inf,
                      anal.int, min.inf = 0.25, nsamp = 5000, alpha = 0.025,
                      conf = 1 - 2 * alpha, ratio = control.rate / test.rate,
                      use = 1, use.rci = use, trunc = alpha / 50,
                      gft = function(rx, cr, tr, ...)
                        rexp(length(rx)) / (c(cr, tr)[rx + 1]), ...) {
  trcv <- if (trunc > 0 & trunc < 1)
    -qnorm(trunc) else 1000
  nn <- 2 * Npg
  rx <- c(rep(0, Npg), rep(1, Npg))
  out <- matrix(0, 5L, nsamp)
  for (i in 1:nsamp) {
    # generate data
    enter <- runif(nn) * acc.per
    # for non-exponential data, replace the following line
    ftime <- gft(rx, control.rate, test.rate, ...)
    atm <- quantile(enter + ftime, probs = total.inf / nn + 1e-8)
    atmin <- quantile(enter + ftime, probs = min.inf * total.inf / nn + 1e-8)
    at <- atmin - anal.int
    i.times <- NULL
    stop <- FALSE
    styp <- 0
    while (!stop) {
      at <- min(at + anal.int, atm)
      sub <- enter <= at
      enters <- enter[sub]
      rxs <- rx[sub]
      time <- enters + ftime[sub]
      status <- ifelse(time <= at, 1, 0)
      time <- pmin(time, at) - enters
      it <- sum(status) / total.inf
      if (it >= min.inf) {
        i.times <- c(i.times, it)
        bndry <- min(rev(sequse(pmin(i.times, 1), alpha, use))[1L], trcv)
        z <- do.call('survdiff', list(formula = Surv(time, status) ~ rxs,
                                      data = data.frame(time, status, rxs)))
        z <- (z$obs[1L] - z$exp[1L]) / sqrt(z$var[1L, 1L])
        if ( z>= bndry) {
          styp <- 1
          stop <- TRUE
        } else {
          if (1 > conf & conf > 0) {
            u <- rci(time, status, rxs, pmin(i.times, 1), conf = conf, use = use.rci)
            if (u[3L] < ratio) {
              styp <- 2
              stop <- TRUE
            }
          }
        }
        if (max(i.times) >= 1 & styp < 1) {
          styp <- 3
          stop <- TRUE
        }
      }
    }
    out[, i] <- c(styp, z, it, at, length(i.times))
  }
  class(out) <- 'powlgrnk6'
  out
}

#' Print a summary of the output from powlgrnk6
#'
#' Prints a summary of the output from \code{\link{powlgrnk6}}.
#'
#' @details
#' Prints the proportion of samples rejecting the null, the proportion stopping
#' in favor of the null, the average information and calendar times at
#' termination, and the distribution of the total number of analyses.
#'
#' @param x A matrix of class \code{\link{powlgrnk6}}
#' @param ... Not used
#'
#' @return
#' No value is returned - used for side-effect of printing to console.
#'
#' @seealso
#' \code{\link{powlgrnk6}}
#'
#' @keywords design survival
#'
#' @export

print.powlgrnk6 <- function(x, ...) {
  rejp <- mean(x[1L, ] == 1)
  cat('Total rejection probability:', format(rejp), '\n')
  cat('Standard error:', format(sqrt(rejp * (1 - rejp) / ncol(x))), '\n')
  h0 <- mean(x[1L, ] == 2)
  cat('Stopped in favor of H0:', format(h0), '\n')
  cat('Standard error:', format(sqrt(h0 * (1 - h0) / ncol(x))), '\n')
  cat('Ave inf time at termination:', format(mean(x[3L, ])), '\n')
  cat('Ave calendar time at termination:', format(mean(x[4L, ])), '\n')
  cat('Distribution of number of analyses:\n')
  print(table(x[5L, ]))
  invisible(NULL)
}