R/scABEL.ad.R

In PowerTOST: Power and Sample Size for (Bio)Equivalence Studies

#-----------------------------------------------------------------------
# Iteratively adjust alpha to maintain the TIE <= nominal alpha for
# partial and full replicate design and scaled ABE via simulated
# (empirical) power
#
# Author: Helmut Schuetz
#-----------------------------------------------------------------------
scABEL.ad <-function(alpha = 0.05, theta0, theta1, theta2, CV,
                     design = c("2x3x3", "2x2x4", "2x2x3"), regulator, n,
                     alpha.pre = 0.05, imax = 100, tol, print = TRUE,
                     details = FALSE, setseed = TRUE, nsims = 1e6,
                     sdsims = FALSE, progress)
{
  env <- as.character(Sys.info()[1]) # get info about the OS
  ifelse ((env == "Windows") || (env == "Darwin"), flushable <- TRUE,
          flushable <- FALSE) # suppress flushing on other OS's
  # acceptance range defaults
  if (missing(theta1) && missing(theta2)) theta1 <- 0.8
  if (missing(theta1)) theta1 = 1/theta2
  if (missing(theta2)) theta2 = 1/theta1
  # check theta0
  if (missing(theta0)) theta0 <- 0.9
  if (theta0 < theta1 || theta0 > theta2)
    stop("theta0 must be within [theta1, theta2]")
  # check regulator arg
  if(missing(regulator)) regulator <- "EMA"
  reg <- reg_check(regulator, choices=c("EMA", "HC", "GCC", "FDA"))
  if (reg$name %in% c("HC", "FDA") && sdsims) {
    stop("Subject data simulations are not supported for regulator =\n",
         "       \"HC\" or \"FDA\".")
  }
  # set iteration tolerance for uniroot().
  if (missing(tol)) tol <- 1e-6
  design <- match.arg(design)
  if (missing(CV)) stop("CV must be given!")
  CVwT <- CV[1]
  if (length(CV) == 2) CVwR <- CV[2] else CVwR <- CVwT
  if(missing(progress)) progress <- FALSE
  no <- 0 # simulation counter
  if (details) ptm <- proc.time()
  if (missing(n) || any(is.na(n))) {
    if (is.na(alpha.pre) || (alpha.pre != alpha)) {
      al <- alpha.pre # If pre-specified, use alpha.pre
    } else {
      al <- alpha    # If not, use alpha (commonly 0.05)
    }
    # sample size for targetpower 0.8
    if (!reg$name == "FDA") { # EMA or HC
      n <- sampleN.scABEL(alpha = al, CV = CV, theta0 = theta0,
                          theta1 = theta1, theta2 = theta2,
                          design = design, regulator = reg,
                          imax = imax, print = FALSE, details = FALSE,
                          setseed = setseed)[["Sample size"]]
    } else {                 # FDA
      n <- sampleN.RSABE(alpha = al, CV = CV, theta0 = theta0,
                         theta1 = theta1, theta2 = theta2,
                         design = design, regulator = reg,
                         imax = imax, print = FALSE, details = FALSE,
                         setseed = setseed)[["Sample size"]]
    }
    if (is.na(n))
      stop("Sample size search in sampleN.scABEL() or sampleN.RSABE() failed.\n",
           "       Restart with an explicitly high n (>1000).")
    no <- 1e5
  }
  if (sum(n) < 6) stop("Sample size too low.")
  if (alpha.pre > alpha) {
    warning("alpha.pre > alpha does not make sense.\n",
            "alpha.pre was set to alpha.")
    alpha.pre <- alpha
  }
  seqs <- as.numeric(substr(design, 3, 3)) # subjects / sequence
  if (length(n) == 1) n <- nvec(n, seqs)   # vectorize n
  
  # here we go!
  TIE <- pwr <- rep(NA, 2) # initialize vectors: TIE and pwr
  alpha.adj <- NA          # adjusted alpha
  # Finds adjusted alpha which gives TIE as close as possible to alpha.
  # Simulate underlying statistics (if sdsims=FALSE)
  opt1 <- function(x) {
    if (reg$name == "FDA") { # FDA
      power.RSABE(alpha = x, CV = CV, theta0 = U, n = n,
                  regulator = reg, design = design,
                  nsims = nsims, setseed = setseed) - alpha
    } else {                # EMA, HC, or GCC
      power.scABEL(alpha = x, CV = CV, theta0 = U, n = n,
                   regulator = reg, design = design,
                   nsims = nsims, setseed = setseed) - alpha
    }
  }
  # Simulate subject data (if sdsims=TRUE)
  opt2 <- function(x) {
    power.scABEL.sdsims(alpha = x, CV = CV, theta0 = U, n = n,
                        regulator = reg, design = design,
                        nsims = nsims, setseed = setseed,
                        progress = progress) - alpha
  }
  sig <- binom.test(x = round(alpha*nsims, 0), n = nsims,
                    alternative = "less",
                    conf.level = 1 - alpha)$conf.int[2]
  method <- "ABE"
  if (CVwR > reg$CVswitch) {
    if (reg$name == "FDA") {
      method <- "RSABE"
    } else {
      method <- "ABEL" 
    }
  }
  limits <- as.numeric(scABEL(CV = CVwR, regulator = reg$name))
  U <- limits[2] # Simulate at the upper (expanded) limit. For CVwR
  # 30% that's 1.25. Due to the symmetry simulations
  # at the lower limit (0.80) should work as well.
  if (alpha.pre != alpha) {
    al <- alpha.pre # If pre-specified, use alpha.pre.
  } else {
    al <- alpha     # If not, use alpha (commonly 0.05).
  }
  designs <- c("2x2x4", "2x2x3", "2x3x3")
  type    <- c("4 period full replicate",
               "3 period full replicate",
               "partial replicate")
  if (print) { # Show input to keep the spirits of the user high.
    if (sdsims) cat("Be patient. Simulating subject data will take a good while!\n\n")  
    cat("+++++++++++ scaled (widened) ABEL ++++++++++++\n")
    cat("         iteratively adjusted alpha\n")
    if (reg$name %in% c("EMA", "GCC")) cat("   (simulations based on ANOVA evaluation)\n")
    if (reg$name %in% c("HC", "FDA")) cat("(simulations based on intra-subject contrasts)\n")
    cat("----------------------------------------------\n")
    cat("Study design: ")
    cat(paste0(design, " (", type[match(design, designs)], ")\n"))
    cat("log-transformed data (multiplicative model)\n")
    cat(formatC(nsims, format = "d", big.mark = ",", decimal.mark = "."),
        "studies in each iteration simulated.\n\n")
    txt <- paste0("CVwR ", sprintf("%.4g", CVwR))
    txt <- paste0(txt, ", CVwT ", sprintf("%.4g", CVwT), ", ")
    cat(paste0(txt, "n(i) ", paste0(n, collapse = "|"), " (N ", sum(n),
               ")\n"))
    txt <- paste0("Nominal alpha                 : ", signif(alpha, 5))
    if (!is.na(alpha.pre) && (alpha.pre != alpha)) {
      txt <- paste0(txt, ", pre-specified alpha ", alpha.pre, "\n")
    } else {
      txt <- paste(txt, "\n")
    }
    cat(txt)
    cat("True ratio                    :", sprintf("%.4f", theta0), "\n")
    cat(paste0("Regulatory settings           : ", reg$name, " (",
               method, ")\n"))
    
    # better theta1, theta2 as BE limits, PE constraint?
    if (CVwR <= reg$CVswitch) {
      cat("Switching CVwR                : ", reg$CVswitch, "\n",
          "BE limits                     : 0.8000 ... 1.2500\n", sep="")
    } else {
      cat(paste("Switching CVwR                :", reg$CVswitch,
                "\nRegulatory constant           :", signif(reg$r_const, 5), "\n"))
      if (reg$name == "FDA") { # FDA
        cat(sprintf("%s             : %.4f%s%.4f%s", "Implied BE limits",
                    limits[1], " ... ", limits[2], "\n"))
      } else {                 # EMA or HC
        if (reg$name %in% c("EMA", "HC")) {
          cat(sprintf("%s               : %.4f%s%.4f%s", "Expanded limits",
                      limits[1], " ... ", limits[2], "\n"))
        } else {               # GCC
          cat(sprintf("%s               : %.4f%s%.4f%s", "Widened limits ",
                      limits[1], " ... ", limits[2], "\n"))
        }
      }
    }
    if (!regulator %in% c("FDA", "GCC"))
      cat("Upper scaling cap             : CVwR >", reg$CVcap, "\n")
    cat("PE constraints                : 0.8000 ... 1.2500\n")
    if (progress) cat("Progress of each iteration:\n")
    if (flushable) flush.console() # advance console output.
  }
  if (!sdsims) { # simulate underlying statistics
    if (reg$name == "FDA") { # FDA
      TIE[1] <- power.RSABE(alpha = al, CV = CV, theta0 = U, n = n,
                            design = design, regulator = reg,
                            nsims = nsims, setseed = setseed)
    } else {                 # EMA, HC, GCC
      TIE[1] <- power.scABEL(alpha = al, CV = CV, theta0 = U, n = n,
                             design = design, regulator = reg,
                             nsims = nsims, setseed = setseed)
    }
    no <- no + nsims
    if (reg$name == "FDA") { # FDA
      pwr[1] <- power.RSABE(alpha = al, CV = CV, theta0 = theta0,
                            n = n, design = design, regulator = reg,
                            setseed = setseed)
    } else {                  # EMA, HC, GCC
      pwr[1] <- power.scABEL(alpha = al, CV = CV, theta0 = theta0,
                             n = n, design = design, regulator = reg,
                             setseed = setseed)
    }
    no <- no + 1e5
  } else { # simulate subject data
    TIE[1] <- power.scABEL.sdsims(alpha = al, CV = CV, theta0 = U, n = n,
                                  design = design, regulator = reg,
                                  nsims = nsims, setseed = setseed,
                                  progress = progress)
    no <- no + nsims
    pwr[1] <- power.scABEL.sdsims(alpha = al, CV = CV, theta0 = theta0,
                                  n = n, design = design, regulator = reg,
                                  setseed = setseed, progress = progress)
    no <- no + 1e5
  }
  if (TIE[1] > alpha) { # adjust only if needed (> nominal alpha)
    if (!sdsims) { # simulate underlying statistics
      x <- uniroot(opt1, interval = c(0, alpha), tol = tol)
    } else { # simulate subject data
      x <- uniroot(opt2, interval = c(0, alpha), tol = tol)
    }
    alpha.adj <- x$root
    if (!sdsims) { # simulate underlying statistics
      if (reg$name == "FDA") { # FDA
        TIE[2] <- power.RSABE(alpha = alpha.adj, CV = CV, theta0 = U, n = n,
                              design = design, regulator = reg,
                              nsims = nsims, setseed = setseed)
        pwr[2] <- power.RSABE(alpha = alpha.adj, CV = CV, theta0 = theta0,
                              n = n, design = design, regulator = reg,
                              setseed = setseed)
      } else {                # EMA, HC, GCC
        TIE[2] <- power.scABEL(alpha = alpha.adj, CV = CV, theta0 = U, n = n,
                               design = design, regulator = reg,
                               nsims = nsims, setseed = setseed)
        pwr[2] <- power.scABEL(alpha = alpha.adj, CV = CV, theta0 = theta0,
                               n = n, design = design, regulator = reg,
                               setseed = setseed)
      }
    } else { # simulate subject data
      TIE[2] <- power.scABEL.sdsims(alpha = alpha.adj, CV = CV, theta0 = U,
                                    n = n, design = design, regulator = reg,
                                    nsims = nsims, setseed = setseed,
                                    progress = progress)
      pwr[2] <- power.scABEL.sdsims(alpha = alpha.adj, CV = CV,
                                    theta0 = theta0, n = n, design = design,
                                    regulator = reg, setseed = setseed,
                                    progress = progress)
    }
  }
  if (!is.na(alpha.adj)) no <- no + nsims*x$iter
  if (details) run.time <- proc.time() - ptm
  if (print) { # fetch and print results
    txt <- paste0("Empiric TIE for alpha ", sprintf("%.4f", al), "  : ",
                  sprintf("%.5f", TIE[1]))
    if (TIE[1] > alpha || alpha.pre != alpha) {
      rel.change <- 100*(TIE[1] - alpha)/alpha
      if (details) {
        txt <- paste0(txt, " (rel. change of risk: ",
                      sprintf("%+1.3g%%", rel.change), ")")
      }
    }
    if (!is.na(pwr[1])) {
      pwr.unadj <- pwr[1]
      txt <- paste0(txt, "\nPower for theta0 ", sprintf("%.4f", theta0),
                    "       : ", sprintf("%.3f", pwr.unadj))
    }
    if (TIE[1] > alpha) {
      txt <- paste0(txt, "\nIteratively adjusted alpha    : ",
                    sprintf("%.5f", alpha.adj),
                    "\nEmpiric TIE for adjusted alpha: ",
                    sprintf("%.5f", TIE[2]))
      if (!is.na(pwr[2])) {
        pwr.adj <- pwr[2]
        txt <- paste0(txt, "\nPower for theta0 ",sprintf("%.4f", theta0),
                      "       : ", sprintf("%.3f", pwr.adj))
        if (details) {
          txt <- paste0(txt, " (rel. impact: ", sprintf("%+1.3g%%",
                                                        100*(pwr[2] - pwr[1])/pwr[1]), ")")
        }
        txt <- paste(txt, "\n\n")
      } else {
        txt <- paste0(txt, "\n\n")
      }
      if (details) {
        txt <- paste0(txt, "Runtime    : ", signif(run.time[3], 3),
                      " seconds\nSimulations: ",
                      formatC(no, format = "d", big.mark = ",",
                              decimal.mark = "."), " (",
                      (no - no %% nsims - nsims) / nsims,
                      " iterations)\n\n")
      }
    } else {
      txt <- paste0(txt, "\nTIE \u2264 nominal alpha; ")
      ifelse(alpha.pre == alpha,
             txt <- paste0(txt, "no adjustment of alpha is required.\n\n"),
             txt <- paste0(txt, "the chosen pre-specified alpha is ",
                           "justified.\n\n"))
    }
    cat(txt)
  } else { # print=FALSE
    # Prepare and return list of results.
    res <- list(regulator = reg$name, method = method, design = design,
                type = type[match(design, designs)], eval = reg$est_method,
                alpha = alpha,
                alpha.pre = ifelse(alpha.pre == alpha, NA, alpha.pre),
                CVwT = CV[1], CVwR = ifelse(length(CV)==1, CV[1], CV[2]),
                N = sum(n), theta0 = theta0, TIE.unadj = signif(TIE[1], 5),
                rel.change = ifelse(!is.na(alpha.adj),
                                    signif(100*(TIE[1] - alpha)/alpha, 5), NA),
                pwr.unadj = signif(pwr[1], 5), alpha.adj = signif(alpha.adj, 5),
                TIE.adj = signif(TIE[2], 5), pwr.adj = signif(pwr[2], 5),
                rel.loss = ifelse(!is.na(pwr[2]),
                                  signif(100*(pwr[2] - pwr[1])/pwr[1], 5), NA),
                sims = no)
    return(res)
  }
}