Helmut
★★★
avatar
Homepage
Vienna, Austria,
2021-01-27 17:50
(36 d 04:19 ago)

Posting: # 22194
Views: 465
 

 Sensitivity analysis (ABE, ABEL, RSABE) [R for BE/BA]

Dear all,

following this post another [image]-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 (CVwT ≡ CVwR) 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 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 (θ0 > 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 θ0, it would need a substantially lower CV to compensate for a worse θ0. 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’ θ0s allow for higher CVs. That’s shown in the upper right quadrants. However, if the CV gets too large, even a θ0 of 1 might not always give the target power.
In the upper left quadrants are the worst case combinations (θ0 < 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 θ0 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")

    [image]

  • 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")


    [image]

  • ABEL (Gulf Cooperation Council)
    sensitivity(CV = 0.4, theta0 = 0.9, do.rate = 0.15, design = "2x2x4",
                SABE = TRUE, regulator = "GCC")


    [image]

  • RSABE (U.S. FDA, China CDE)
    sensitivity(CV = 0.4, theta0 = 0.9, do.rate = 0.15, design = "2x2x4",
                SABE = TRUE, regulator = "FDA")


    [image]
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            #
#####################################################


Dif-tor heh smusma 🖖
Helmut Schütz
[image]

The quality of responses received is directly proportional to the quality of the question asked. 🚮
Science Quotes
Activity
 Admin contact
21,366 posts in 4,462 threads, 1,493 registered users;
online 7 (0 registered, 7 guests [including 4 identified bots]).
Forum time: Thursday 22:10 UTC (Europe/Vienna)

Being really good at C++ is like being really good
at using rocks to sharpen sticks.    Thant Tessman

The Bioequivalence and Bioavailability Forum is hosted by
BEBAC Ing. Helmut Schütz
HTML5