Bioequivalence and Bioavailability Forum

Dear all,

following this post another -script (a.k.a. jack of all trades device). Now reference-scaling is supported (including the approach of the Guld Cooperation Council). Note that you need at least version 1.5-3 of package PowerTOST.

• ABE
  • If less than 12 eligible subjects might remain (due to low CV and/or theta0 close to 1), the sample size is increased accordingly.
• SABE (ABEL, GCC, RSABE)
  • Like in the function pa.scABE() only homoscedasticity (CV_wT ≡ CV_wR) is implemented.
    
  • If for the EMA’s ABEL design = "2x2x3" (3-period full replicate designs TRT|RTR or TRR|RTT) less than 12 eligible subjects might remain in the sequence RTR or TRR, the sample size is increased accordingly.
    
  • If for the FDA’s RSABE the estimated sample size n is less than 24, it is increased to 24 (al­though the number of eligible subjects might still be lower).

Examples (all assumed θ₀ 0.9, CV 40%, target power 80%, 4-period 2-sequence full replicate design, anticipated dropout-rate 15%).
The lower right quadrants of each panel show ‘nice‘ combinations (θ₀ > assumed and CV < assumed). Higher power than desired, great. Sometimes you will see a statistically significant treatment effect. No worries about that (not clinically relevant).
The other combinations are tricky. Since power is most sensitive to θ₀, it would need a substantially lower CV to compensate for a worse θ₀. Have a look at the 0.80 contour lines in the lower left quadrant of the first panels (no dropouts).
On the other hand, ‘better’ θ₀s allow for higher CVs. That’s shown in the upper right quadrants. However, if the CV gets too large, even a θ₀ of 1 might not always give the target power.
In the upper left quadrants are the worst case combinations (θ₀ < assumed and CV > assumed). It might still be possible to show BE though with a lower chance (power < 0.80).
Dropouts don’t hurt that much.

The basic idea of reference-scaling was to maintain power independent from the CV. Hence, ideally for any θ₀ and sample size power would be straight vertical lines. Not that bad for ABEL and RSABE (at high CVs the PE constraint – and for ABEL the upper cap of scaling – bend the curves to the right).

• ABE (all agencies)
  sensitivity(CV = 0.4, theta0 = 0.9, do.rate = 0.15, design = "2x2x4")
  
  
  
  
• ABEL (EMA, the WHO, ASEAN States, Australia, Brazil, Chile, the East African Community, Egypt, the Eurasian Economic Union, New Zealand, the Russian Federation)
  sensitivity(CV = 0.4, theta0 = 0.9, do.rate = 0.15, design = "2x2x4",
            SABE = TRUE, regulator = "EMA")
  
  
  
  
• ABEL (Gulf Cooperation Council)
  sensitivity(CV = 0.4, theta0 = 0.9, do.rate = 0.15, design = "2x2x4",
            SABE = TRUE, regulator = "GCC")
  
  
  
  
• RSABE (U.S. FDA, China CDE)
  sensitivity(CV = 0.4, theta0 = 0.9, do.rate = 0.15, design = "2x2x4",
            SABE = TRUE, regulator = "FDA")
  
  

If you discover bugs or have suggestions for improvement, let me know.

---

check.packages <- function() {
  # Check whether the packages are installed in the minimum required version.
  # If yes, attach them. If not, stop and ask the user for download/installation.
  packs <- c("PowerTOST", "lattice", "latticeExtra") # Required packages
  inst <- installed.packages()[packs, "Version"]     # Are they installed?
  if (length(inst) == 0) {                           # No
    stop ("Please download/install the packages \'PowerTOST\',",
          "\n       \'lattice\', and \'latticeExtra\' from CRAN.")
  } else {                                           # Yes
    suppressMessages(require(PowerTOST))             # Attach them
    if (inst[1] < "1.5-3")                           # Check version
      stop ("At least version 1.5-3 of package \'PowerTOST\' required.",
            "\n       Please download/install the current version from CRAN.")
    suppressMessages(require(lattice))
    suppressMessages(require(latticeExtra))
  }
}
balance <- function(n, sequences) {
  # Round up to get balanced sequences for potentially unbalanced case.
  return(as.integer(sequences * (n %/% sequences + as.logical(n %% sequences))))
}
adjust.dropouts <- function(n, do.rate, sequences = 2) {
  # To be dosed subjects which should result in n eligible subjects based
  # on the anticipated droput-rate.
  n <- as.integer(balance(n / (1 - do.rate), sequences = sequences))
  return(n)
}
power.TOST.vectorized <- function(CVs, theta0s, ...) {
  # Supportive function, since only theta0 in power.TOST()
  # can be vectorized (not the CV).
  # Returns a matrix of power with rows CV and columns theta0.
  power <- matrix(ncol = length(CVs), nrow = length(theta0s),
                  dimnames = list(c(paste0("theta0.", 1:length(theta0s))),
                                  c(paste0("CV.", 1:length(CVs)))))
  for (j in seq_along(CVs)) {
    power[, j] <- suppressMessages( # for unbalanced cases
                    power.TOST(CV = CVs[j], theta0 = theta0s, ...))
  }
  return(power)
}
power.SABE.vectorized <- function(CVs, theta0s, regulator, ...) {
  # Supportive function, since power.scABEL() and power.RSABE()
  # are not vectorized.
  # Returns a matrix of power with rows CV and columns theta0.
  power <- matrix(ncol = length(CVs), nrow = length(theta0s),
                  dimnames = list(c(paste0("theta0.", 1:length(theta0s))),
                                  c(paste0("CV.", 1:length(CVs)))))
  pb    <- txtProgressBar(min = 0, max = 1, char = "\u2588", style = 3)
  i     <- 0
  for (j in seq_along(theta0s)) {
    for (k in seq_along(CVs)) {
      i <- i + 1
      if (!regulator == "FDA") {
        power[j, k] <- suppressMessages( # for unbalanced cases
                         power.scABEL(CV = CVs[k], theta0 = theta0s[j],
                                      regulator = regulator, ...))
      } else {
        power[j, k] <- suppressMessages( # for unbalanced cases
                         power.RSABE(CV = CVs[k], theta0 = theta0s[j],
                                     regulator = regulator, ...))
      }
      setTxtProgressBar(pb, i/(length(theta0s)*length(CVs)))
    }
  }
  close(pb)
  return(power)
}
sensitivity <- function(alpha = 0.05, CV, CV.lo, CV.hi, theta0,
                        theta0.lo, do.rate, target = 0.80,
                        design = "2x2", theta1, theta2, SABE = FALSE,
                        regulator = "EMA", nsims = 1e5, mesh = 25) {
  check.package()
  # use alpha = 0.5 for assessing only the PE (Health Canada: Cmax)
  if (alpha <= 0 | alpha > 0.5)
    stop("alpha ", alpha, " does not make sense.")
  if (missing(CV))
    stop("CV must be given.")
  if (length(CV) > 1)
    stop("CV must be scalar (only homoscedasticity supported).")
  if (missing(CV.lo))
    CV.lo <- CV * 0.8
  if (missing(CV.hi))
    CV.hi <- CV * 1.25
  if (missing(theta0)) {
    if (!SABE) {
      theta0 <- 0.95
    } else {
      theta0 <- 0.90
    }
  }
  if (theta0 > 1)
    stop("theta0 >1 not implemented.")
  if (missing(theta1) & missing(theta2))
    theta1 <- 0.8
  if (!missing(theta1) & missing(theta2))
    theta2 <- 1/theta1
  if (missing(theta1) & !missing(theta2))
    theta1 <- 1/theta2
  if (missing(theta0.lo))
    theta0.lo <- theta0 * 0.95
  if (theta0.lo < theta1) {
    message("theta0.lo ", theta0.lo, " < theta1 does not make sense. ",
            "Changed to theta1.")
    theta0.lo <- theta1
  }
  theta0.hi <- 1
  if (missing(do.rate))
    stop("do.rate must be given.")
  if (do.rate < 0)
    stop("do.rate", do.rate, " does not make sense.")
  if (target <= 0.5)
    stop("Target ", target, " does not make sense. Toss a coin instead.")
  if (target >= 1)
    stop("Target ", target, " does not make sense.")
  d.no <- PowerTOST:::.design.no(design)
  if (is.na(d.no))
    stop("design '", design, "' unknown.")
  steps <- known.designs()[d.no, "steps"]   if (SABE) {
    if (!design %in% c("2x3x3", "2x2x4", "2x2x3"))
      stop("design '", design, "' not implemented for SABE.")
    if (!regulator %in% c("EMA", "HC", "GCC", "FDA"))
      stop("regulator must be any of \"EMA\", \"HC\", \"GCC\", \"FDA\".")
    if (nsims < 1e5)
      message("Use <100,000 simulations only for a preliminary look.")
  }
  if (mesh <= 10) {
    message("Too wide mesh is imprecise. Increased to 25.")
    mesh <- 25
  }
  CVs     <- seq(CV.lo, CV.hi, length.out = mesh)
  theta0s <- seq(theta0.lo, theta0.hi, length.out = mesh)
  if (!SABE) {                  # ABE
    n <- sampleN.TOST(alpha = alpha, CV = CV, theta0 = theta0,
                      theta1 = theta1, theta2 = theta2,
                      targetpower = target, design = design,
                      print = FALSE, details = FALSE)[["Sample size"]]
    if (n < 12) n <- 12         # acc. to guidelines
  } else {                      # SABE
    if (!regulator == "FDA") {  # ABEL (EMA, ..., Health Canada, GCC)
      n <- sampleN.scABEL(alpha = alpha, CV = CV, theta0 = theta0,
                          theta1 = theta1, theta2 = theta2,
                          targetpower = target, design = design,
                          regulator = regulator, nsims = nsims,
                          print = FALSE, details = FALSE)[["Sample size"]]
    } else {                    # RSABE (FDA, China)
      n <- sampleN.RSABE(alpha = alpha, CV = CV, theta0 = theta0,
                         theta1 = theta1, theta2 = theta2,
                         targetpower = target, design = design,
                         regulator = regulator, nsims = nsims,
                         print = FALSE, details = FALSE)[["Sample size"]]
      if (n < 24) n <- 24 # acc. to guidance
    }
  }
  ns      <- as.integer(n:floor(n * (1 - do.rate)))
  if (!SABE & min(ns) < 12) {
    n  <- adjust.dropouts(n, do.rate, steps) # acc. to guidelines
    ns <- as.integer(n:floor(n * (1 - do.rate)))
    message("Increased sample size preventing <12 eligible subjects.")
  }
  if (SABE & regulator == "EMA" & design == "2x2x3" & min(ns) < 24) {
    n  <- adjust.dropouts(n, do.rate, steps) # acc. to Q & A document
    ns <- as.integer(n:floor(n * (1 - do.rate)))
    message("Increased sample size preventing <12 eligible subjects ",
            "in replicated R-sequence!")
  }
  if (SABE) {
     cat(paste0("Preparing ", length(ns), " panels for SABE: regulator = ",
         regulator, ", design = ", design, "\n"))
  }
  grid   <- expand.grid(x = theta0s, y = CVs, z = NA, n = as.factor(ns))
  grid$n <- factor(grid$n, levels = rev(levels(grid$n)))
  for (j in seq_along(ns)) {
    if (!SABE) {
      grid$z[grid$n == ns[j]] <- as.vector(
                                   power.TOST.vectorized(
                                     alpha = alpha, CV = CVs, theta0 = theta0s,
                                     n = ns[j], design = design))
    } else {
      cat(sprintf("  Panel %i (n = %i)%s", j, ns[j], "\n"))
      grid$z[grid$n == ns[j]] <- as.vector(
                                   power.SABE.vectorized(
                                     alpha = alpha, CV = CVs, theta0 = theta0s,
                                     n = ns[j], design = design,
                                     regulator = regulator, nsims = nsims))
    }
  }
  res <- data.frame(n = ns, power = grid$z[grid$x == theta0 & grid$y == CV],
                    do.rate = abs(1 - ns / n))
  custom                               <- trellis.par.get()
  custom$layout.widths$left.padding    <-  0.3
  custom$layout.widths$right.padding   <- -1.85
  custom$layout.heights$top.padding    <- -3
  custom$layout.heights$bottom.padding <- -0.8
  trellis.par.set(custom)
  trellis.device(device = windows, theme = custom,
                 width = 6.5, height = 6.5, record = TRUE)
  st <- strip.custom(strip.names = TRUE, strip.levels = TRUE, sep = " = ")
  ct <- length(pretty(grid$z, 10))
  fg <- contourplot(z ~ x * y | n, data = grid, as.table = TRUE, strip = st,
                    labels = list(cex = 0.75), label.style = "align",
                    xlab = expression(italic(theta)[0]), ylab = "CV", cuts = ct)
  plot(fg)
  cat("Results for theta0 =", theta0,
      "and CV =", CV, "\n"); print(signif(res, 5), row.names = FALSE)
}
#####################################################
# Specification of the study (mandatory values)     #
#####################################################
CV      <- 0.25 # assumed CV                        #
do.rate <- 0.10 # anticipated dropout-rate (10%)    #
#####################################################
# defaults (if not provided in named arguments)     #
#   alpha      = 0.05         common                #
#   CV.lo      = CV*0.80      best case             #
#   CV.hi      = CV*1.25      worst case            #
#   theta0                    assumed T/R-ratio     #
#                             0.95 for ABE and      #
#                             0.90 for SABE         #
#   theta0.lo  = theta0*0.95  worst case            #
#   target     = 0.80         target power          #
#   theta1     = 0.80         lower BE limit        #
#   theta2     = 1.25         upper BE limit        #
#   design     = "2x2"        in known.designs()    #
#   SABE       = FALSE        conventional ABE      #
#                             TRUE for ABEL/RSABE   #
#   regulator  = "EMA"        only observed if      #
#                             SABE == TRUE          #
#                             possible "EMA", "HC", #
#                             "GCC", "FDA"          #
#   nsims      = 1e5          only observed if      #
#                             SABE == TRUE          #
#   mesh       = 25           resolution            #
#####################################################