Bioequivalence and Bioavailability Forum

Hi Achievwin,

❝ R PowerTOST has a neat code that gives tabular listing of N and power attained for Scaled ABE using the macro below.

❝

❝ Is there a similar code for ABE for 2x2 or parallel design studies. Any help is much appreciated.

Not out of the box. Consider the functions pa.ABE() and pa.scABE() – please in the R-console and not in RStudio.
If you are interested in a table, try the lengthy -script at the end, which works for all covered designs (ABE, and if applicable, RSABE or ABEL, or the FDA’s RSABE for NTIDs).

Examples:

power.table(CV = 0.25) # using all defaults

design        : 2x2x2
alpha         : 0.0500
method        : Average Bioequivalence (ABE)
power method  : exact
CVw           : 0.2500
theta0        : 0.9500
theta1        : 0.8000
theta2        : 1.2500
targetpower   : 0.8000
n             : 28 (achieved power 0.80744)

  n   power
 24 0.73912
 26 0.77606
 28 0.80744

power.table(CV = 0.4, design = "parallel")

design        : parallel
alpha         : 0.0500
method        : Average Bioequivalence (ABE)
power method  : exact
CVtotal       : 0.4000
theta0        : 0.9500
theta1        : 0.8000
theta2        : 1.2500
targetpower   : 0.8000
n             : 130 (achieved power 0.80351)

   n   power
 104 0.70576
 106 0.71491
 108 0.72375
 110 0.73229
 112 0.74053
 114 0.74850
 116 0.75620
 118 0.76365
 120 0.77085
 122 0.77781
 124 0.78456
 126 0.79108
 128 0.79740
 130 0.80351

power.table(CV = 0.46, theta0 = 0.93, target = 0.9, ABE = FALSE, design = "2x2x3",
            regulator = "FDA", lower = 0.6) # your example

design        : 2x2x3
alpha         : 0.0500
regulator     : FDA
method        : Reference-scaled Average Bioequivalence (RSABE)
simulations   : intra-subject contrasts
number of sims: 100,000
CVswitch      : 0.3000
Reg. constant : 0.8925742
PE constraints: 0.8000, 1.2500
CVw           : 0.4600
theta0        : 0.9300
targetpower   : 0.9000
n             : 40 (achieved power 0.91323)

  n   power
 24 0.72529
 26 0.76284
 28 0.79686
 30 0.82498
 32 0.84757
 34 0.86782
 36 0.88550
 38 0.89948
 40 0.91323

power.table(CV = 0.4, ABE = FALSE, design = "2x2x4") # default ABEL (EMA),
                                                     # additionally Type I Error

design        : 2x2x4
alpha         : 0.0500
regulator     : EMA
method        : Average Bioequivalence with Expanding Limits (ABEL)
simulations   : ANOVA
number of sims: 100,000
CVswitch      : 0.3000
Reg. constant : 0.76
PE constraints: 0.8000, 1.2500
Cap           : 0.5000
CVw           : 0.4000
theta0        : 0.9000
targetpower   : 0.8000
n             : 30 (achieved power 0.80656)

  n   power      TIE
 24 0.72851 0.059438
 26 0.75868 0.059873
 28 0.78286 0.059994
 30 0.80656 0.059119

power.table(CV = 0.4, ABE = FALSE, design = "2x2x4",
            regulator = "HC") # ABEL for Health Canada

design        : 2x2x4
alpha         : 0.0500
regulator     : HC
method        : Average Bioequivalence with Expanding Limits (ABEL)
simulations   : intra-subject contrasts
number of sims: 100,000
CVswitch      : 0.3000
Reg. constant : 0.76
PE constraints: 0.8000, 1.2500
Cap           : 0.5738
CVw           : 0.4000
theta0        : 0.9000
targetpower   : 0.8000
n             : 32 (achieved power 0.8134)

  n   power
 26 0.73791
 28 0.76636
 30 0.79108
 32 0.81340

power.table(CV = 0.4, ABE = FALSE, design = "2x2x4",
            regulator = "GCC") # the Gulf Cooperation Council’s variant

design        : 2x2x4
alpha         : 0.0500
regulator     : GCC
method        : Average Bioequivalence with widened limits
simulations   : ANOVA
number of sims: 100,000
CVswitch      : 0.3000
Reg. constant : 0.9799758
PE constraints: 0.8000, 1.2500
CVw           : 0.4000
theta0        : 0.9000
L             : 0.7500
U             : 1.3333
targetpower   : 0.8000
n             : 30 (achieved power 0.81111)

  n   power
 24 0.72505
 26 0.75653
 28 0.78532
 30 0.81111

power.table(CV = 0.15, theta1 = 0.9, design = "2x2x4",
            lower = 0.7) # narrower limits

design        : 2x2x4
alpha         : 0.0500
method        : Average Bioequivalence (ABE) (NTID: EMA and others)
power method  : exact
CVw           : 0.1500
theta0        : 0.9750
theta1        : 0.9000
theta2        : 1.1111
targetpower   : 0.8000
n             : 24 (achieved power 0.82627)

  n   power
 18 0.70436
 20 0.75296
 22 0.79299
 24 0.82627

power.table(CV = c(0.17, 0.15), NTID = TRUE, design = "2x2x4",
            lower = 0.7) # NTID for the FDA, CVwT > CVwR

design        : 2x2x4
alpha         : 0.0500
regulator     : FDA or CDE
method        : RSABE (NTID)
simulations   : intra-subject contrasts
number of sims: 100,000
Reg. constant : 1.053605
‘Cap’         : 0.2142
Implied limits: 0.8546, 1.1702
ABE limits    : 0.8000, 1.2500
CVw           : 0.1700, 0.1500 (T, R)
theta0        : 0.9750
targetpower   : 0.8000
n             : 18 (achieved power 0.83074)

  n   power pwr.RSABE pwr.ABE pwr.sratio
 14 0.68508   0.76537 0.99680    0.84246
 16 0.76718   0.83156 0.99908    0.88913
 18 0.83074   0.87845 0.99978    0.92548

---

Sorry, 248 lines of code. May I charge you for three hours of coding and testing?
power.table <- function(alpha = 0.05, CV = CV, theta0, theta1, theta2,
                        design = "2x2x2", ABE = TRUE, NTID = FALSE,
                        target = 0.80, method = "exact", regulator = "EMA",
                        nsims = 1e5, lower = 0.80, upper = 1, TIE = FALSE) {
  ###########################################################################
  # Arguments:                                                              #
  # alpha    : Nominal level of the test (default 0.05)                     #
  # CV       : within-subject (crossovers, paired design), total (parallel) #
  #            can be a two-element vector (for SABE only)                  #
  #            first element CVw of T, second element CVw of R              #
  # theta0   : assumed T/R-ratio (default 0.95 for ABE, 0.90 for SABE,      #
  #            0.975 for NTID and ABE if theta1 = 0.9 or theta1 = 1 / 0.9   #
  # theta1   : lower BE-limit (ABE/NTID), lower PE-constraint (SABE)        #
  #            default 0.80 for ABE/NTID, fixed PE-constraint 0.80 for SABE #
  # theta2   : upper BE-limit (ABE/NTID), upper PE-constraint (SABE)        #
  #            default 1.25 for ABE/NTID, fixed PE-constraint 1.25 for SABE #
  # design   : any supported; see known.designs()                           #
  # ABE      : TRUE (default) for ABE, FALSE for SABE                       #
  # NTID     : FALSE (default), TRUE for RSABE (FDA, CDE)                   #
  # target   : target (desired) power (default 0.8)                         #
  # method   : ABE only; power method (default "exact")                     #
  # nsims    : SABE only; number of simulations for SABE (default 1e5)      #
  # regulator: SABE only; any of "EMA", "FDA", "HC", "GCC" (default "EMA")  #
  # lower    : lower fraction of the estimated sample size (default 0.80)   #
  # upper    : upper fraction of the estimated sample size (default 1)      #
  # TIE      : SABE only; empiric Type I Error (default FALSE)              #
  # ----------------------------------------------------------------------- #
  # Note:                                                                   #
  # Only for the multiplicative model, i.e., in PowerTOST’s functions for   #
  # ABE logscale = TRUE is applied                                          #
  ###########################################################################
  require(PowerTOST)
  balance <- function(n, ns) {
    return(as.integer(ns * (n %/% ns + as.logical(n %% ns))))
  }
  # fair amount of error handling
  if (NTID) {
    ABE       <- FALSE
    regulator <- "FDA"
  }
  if (ABE & length(CV) > 1)
    stop("Give \"CV\" as a one element vector for ABE.")
  if (!ABE & !regulator %in% c("EMA", "FDA", "HC", "GCC"))
    stop("regulator \"", regulator, "\" not supported.")
  if (!ABE & !design %in% c("2x3x3", "2x2x4", "2x2x3"))
    stop("A replicate design is required for reference-scaling.")
  if (NTID & !design %in% c("2x2x4", "2x2x3"))
    stop("A full replicate design is required for reference-scaling.")
  if (missing(theta1) & missing(theta2)) theta1 <- 0.80
  if (missing(theta1)) theta1 <- 1/theta2
  if (missing(theta2)) theta2 <- 1/theta1
  # default theta0
  if (missing(theta0)) {
    if (NTID | theta1 == 0.9 | theta2 == 1 / 0.9) {
      theta0 <- 0.975
    } else {
      ifelse (ABE, theta0 <- 0.95, theta0 <- 0.90)
    }
  }
  if (!method %in% c("exact", "owenq", "mvt", "noncentral",
                     "nct", "central", "shifted"))
    stop("method \"", method, "\" not supported.")
  if (upper <= lower) stop("\"upper\" must be > \"lower\".")
  ifelse (length(CV) > 1, CVwR <- CV[2], CVwR <- CV)
  designs <- known.designs()[, c(2, 5)]
  n.step  <- designs[designs$design == design, "steps"]
  if (NTID) {
    if (!theta1 == 0.8) {
      msg <- paste("theta1 =", sprintf("%.4f", theta1),
                   "not compliant with the guidance.")
      message(msg)
    }
    if (!theta2 == 1.25) {
      msg <- paste("theta2 =", sprintf("%.4f", theta2),
                   "not compliant with the guidance.")
      message(msg)
    }
    n <- sampleN.NTID(alpha = alpha, CV = CV, theta0 = theta0,
                      theta1 = theta1, theta2 = theta2, targetpower = target,
                      design = design, nsims = nsims, details = FALSE,
                      print = FALSE)[["Sample size"]]
  } else {
    if (ABE) {
      if (length(CV) > 1)
        stop("Give \"CV\" as a one element vector for ABE.")
      n <- sampleN.TOST(alpha = alpha, CV = CV, theta0 = theta0,
                        theta1 = theta1, theta2 = theta2, targetpower = target,
                        design = design, method = method,
                        print = FALSE)[["Sample size"]]
    } else {
      if (!design %in% c("2x3x3", "2x2x4", "2x2x3"))
        stop("Replicate design required for reference-scaling.")
      if (regulator == "FDA") {
        n <- sampleN.RSABE(alpha = alpha, CV = CV, theta0 = theta0,
                           targetpower = target, design = design,
                           nsims = nsims, details = FALSE,
                           print = FALSE)[["Sample size"]]
      } else {
        if (!regulator == "HC" & design == "2x3x3" & CV[1] > CV[2]) {
          # subject simulations if partial replicate and heteroscedasticity
          n <- sampleN.scABEL.sdsims(alpha = alpha, CV = CV, theta0 = theta0,
                                     targetpower = target, design = design,
                                     regulator = regulator, nsims = nsims,
                                     details = FALSE,
                                     print = FALSE)[["Sample size"]]
        } else {
          n <- sampleN.scABEL(alpha = alpha, CV = CV, theta0 = theta0,
                              targetpower = target, design = design,
                              regulator = regulator, nsims = nsims,
                              details = FALSE, print = FALSE)[["Sample size"]]
        }
      }
    }
  }
  n.lower <- balance(n * lower, n.step)
  n.upper <- balance(n * upper, n.step)
  x       <- data.frame(n = seq(n.lower, n.upper, n.step), power = NA_real_)
  if (NTID) x <- cbind(x, pwr.RSABE = NA_real_, pwr.ABE = NA_real_,
                          pwr.sratio = NA_real_)
  if (!ABE & TIE) x <- cbind(x, TIE = NA_real_)
  for (j in 1:nrow(x)) {
    if (NTID) {
      tmp <- suppressMessages(
               as.numeric(
                 power.NTID(alpha = alpha, CV = CV, theta0 = theta0,
                            theta1 = theta1, theta2 = theta2, design = design,
                            n = x$n[j], nsims = nsims, details = TRUE)[1:4]))
      x[j, 2:5] <- tmp
    } else {
      if (ABE) {
        x$power[j] <- power.TOST(alpha = alpha, CV = CV, theta0 = theta0,
                                 theta1 = theta1, theta2 = theta2,
                                 design = design, method = method, n = x$n[j])
      } else {
       reg <- reg_const(regulator = regulator)
       U   <- scABEL(CV = CVwR, regulator = regulator)[["upper"]]
       if (regulator == "FDA") {
          x$power[j] <- power.RSABE(alpha = alpha, CV = CV, theta0 = theta0,
                                    design = design, n = x$n[j], nsims = nsims)
          if (TIE) {
            x$TIE[j] <- power.RSABE(alpha = alpha, CV = CV, theta0 = U,
                                    design = design, n = x$n[j], nsims = 1e6)
          }
        } else {
          if (!regulator == "HC" & design == "2x3x3" & CV[1] > CV[2]) {
            x$power[j] <- power.scABEL.sdsims(alpha = alpha, CV = CV,
                                              theta0 = theta0, design = design,
                                              regulator = regulator, n = x$n[j],
                                              nsims = nsims)
            if (TIE) {
              x$TIE[j] <- power.scABEL.sdsims(alpha = alpha, CV = CV,
                                              theta0 = U, design = design,
                                              regulator = regulator,
                                              n = x$n[j], nsims = 1e6)
            }
          } else {
            x$power[j] <- power.scABEL(alpha = alpha, CV = CV, theta0 = theta0,
                                       design = design, regulator = regulator,
                                       n = x$n[j], nsims = nsims)
            if (TIE) {
              x$TIE[j] <- power.scABEL(alpha = alpha, CV = CV, theta0 = U,
                                       design = design, regulator = regulator,
                                       n = x$n[j], nsims = 1e6)
            }
          }
        }
      }
    }
  }
  f1  <- "\n%-14s: %s"
  f2  <- "\n%-14s: %.4f"
  f3  <- "\n%-14s: %.7g"
  f4  <- "\n%-14s: %.4f, %.4f"
  txt <- paste("\ndesign        :", design,
               sprintf(f2, "alpha", alpha))
  if (ABE) {
    txt <- paste(txt, sprintf(f1, "method", "Average Bioequivalence (ABE)"))
    if (theta1 == 0.9) txt <- paste(txt, "(NTID: EMA and others)")
    txt <- paste(txt, sprintf(f1, "power method", method))
  } else {
    if (NTID) {
      r_const <- -log(0.9) / 0.1
      txt <- paste(txt, sprintf(f1, "regulator", "FDA or CDE"))
      txt <- paste(txt, sprintf(f1, "method", "RSABE (NTID)"))
      txt <- paste(txt, sprintf(f1, "simulations", "intra-subject contrasts"))
      txt <- paste(txt, sprintf("\n%-14s:", "number of sims"),
                        formatC(nsims, format = "d", big.mark = ","))
      txt <- paste(txt, sprintf(f3, "Reg. constant", r_const),
                        sprintf(f2, "‘Cap’", se2CV(log(theta2) / r_const)))
      txt <- paste(txt, sprintf(f4, "Implied limits",
                                     exp(-r_const * CV2se(CVwR)),
                                     exp(r_const * CV2se(CVwR))))
      txt <- paste(txt, sprintf(f4, "ABE limits", theta1, theta2))
    } else {
      txt <- paste(txt, sprintf(f1, "regulator", regulator))
      if (regulator == "FDA") {
        txt <- paste(txt, sprintf(f1, "method",
                          "Reference-scaled Average Bioequivalence (RSABE)"))
      } else {
        if (!regulator == "GCC") {
          txt <- paste(txt, sprintf(f1, "method",
                            "Average Bioequivalence with Expanding Limits (ABEL)"))
        } else {
          txt <- paste(txt, sprintf(f1, "method",
                            "Average Bioequivalence with widened limits"))
        }
      }
      if (reg$est_method == "ANOVA") {
        txt <- paste(txt, sprintf(f1, "simulations", "ANOVA"))
      } else {
        txt <- paste(txt, sprintf(f1, "simulations", "intra-subject contrasts"))
      }
      txt <- paste(txt, sprintf("\n%-14s:", "number of sims"),
                        formatC(nsims, format = "d", big.mark = ","))
      txt <- paste(txt, sprintf(f2, "CVswitch", reg$CVswitch),
                        sprintf(f3, "Reg. constant", reg$r_const))
      if (reg$pe_constr) {
        txt <- paste(txt, sprintf(f4, "PE constraints", theta1, theta2))
      }
      if (!regulator %in% c("FDA", "GCC")) {
        txt <- paste(txt, sprintf(f2, "Cap", reg$CVcap))
      }
    }
  }
  if (design == "parallel") {
    txt <- paste(txt, "\nCVtotal       :",
                 paste(sprintf("%.4f", CV), collapse = ", "))
  } else {
    txt <- paste(txt, "\nCVw           :",
                 paste(sprintf("%.4f", CV), collapse = ", "))
  }
  if (length(CV) > 1) {
    txt <- paste(txt, "(T, R)")
  }
  txt <- paste(txt, sprintf(f2, "theta0", theta0))
  if (ABE) {
    txt <- paste(txt, sprintf(f2, "theta1", theta1),
                      sprintf(f2, "theta2", theta2))
  }
  if (regulator == "GCC" & CVwR > 0.3) {
    txt <- paste(txt, sprintf(f2, "L", 0.75), sprintf(f2, "U ", 1 / 0.75))
  }
  txt <- paste(txt, sprintf(f2, "targetpower", target),
                    "\nn             :", n,
               sprintf("(achieved power %.5g)", x$power[x$n == n]), "\n\n")
  cat(txt)
  print(signif(x, 5), row.names = FALSE)
}