R/CalcConPwrExpAndersen.R

In CP: Conditional Power Calculations

CalcConPwrExpAndersen <- function(theta.0,
                                  d1, o1, O1.star, lambda1.hat,
                                  d2, o2, O2.star,
                                  alpha) {

  ## Calculates the conditional power
  ## for the exponential model, i. e.
  ##  S(t) = exp(- lambda * t), lambda > 0, t >= 0,
  ## with the original formulae of the Andersen paper.
  ##
  ## Args:
  ##   Parameters from exponential power calculations.
  ##
  ## Results:
  ##   Returns conditional power values.

  # range of theta for conditional power function
  min <- exp(min(log(theta.0), - log(theta.0)) - 1)
  max <- exp(max(log(theta.0), - log(theta.0)) + 1)

  # choice of suitable stepwidth for theta's
  if ((max - min) / 0.01 <= 1000) {
    theta <- seq(from = min,
                 to   = max,
                 by   = 0.01)
  }
  else {
    theta <- seq(from       = min,
                 to         = max,
                 length.out = 1000)
  }

  # expected value under the null hypothesis
  mu.theta.null     <- (log((d2 + lambda1.hat * O2.star) / (o2 + O2.star))
                         - log((d1 + lambda1.hat * O1.star) / (o1 + O1.star)))
  # standard deviation under the null hypothesis
  sigma.theta.null  <- (sqrt((1 / (o1 + O1.star) + 1 / (o2 + O2.star))
                               / ((d1 + d2) / (o1 + o2))))
  # expected value under the alternative hypothesis
  mu.theta.alter    <- (log((d2 + theta * lambda1.hat * O2.star) / (o2 + O2.star))
                         - log((d1 + lambda1.hat * O1.star) / (o1 + O1.star)))
  # standard deviation under the alternative hypothesis
  sigma.theta.alter <- (sqrt(lambda1.hat * O1.star / (d1 + lambda1.hat * O1.star)^2
                              + theta * lambda1.hat * O2.star / (d2 + theta * lambda1.hat * O2.star)^2))
  # (asymptotically) conditional power function
  gamma.theta       <- (stats::pnorm((stats::qnorm(alpha / 2) * sigma.theta.null
                                - mu.theta.alter) / sigma.theta.alter)
                         + 1 - stats::pnorm((stats::qnorm(1 - alpha / 2) * sigma.theta.null
                                       - mu.theta.alter) / sigma.theta.alter))
  # values at theta.0
  mu.theta.0        <- (log((d2 + theta.0 * lambda1.hat * O2.star) / (o2 + O2.star))
                         - log((d1 + lambda1.hat * O1.star) / (o1 + O1.star)))
  sigma.theta.0     <- (sqrt(lambda1.hat * O1.star / (d1 + lambda1.hat * O1.star)^2
                              + theta.0 * lambda1.hat * O2.star / (d2 + theta.0 * lambda1.hat * O2.star)^2))
  gamma.theta.0     <- (stats::pnorm((stats::qnorm(alpha / 2) * sigma.theta.null
                                - mu.theta.0) / sigma.theta.0)
                         + 1 - stats::pnorm((stats::qnorm(1 - alpha / 2) * sigma.theta.null
                                       - mu.theta.0) / sigma.theta.0))

  return(list(theta, gamma.theta, gamma.theta.0))
}