Achievwin
★    

US,
2022-09-16 06:31
(581 d 02:15 ago)

Posting: # 23300
Views: 2,888
 

 R code for Average BE or parallel design N com­pu­ta­tions [🇷 for BE/BA]

Post reply
Dear All:

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.
(may be a simple thing but I don't know)

Thanks and regards,

Achievwin

sampleN.RSABE(alpha = 0.05, targetpower = 0.80, theta0 = 0.9, theta1 = 0.8, theta2 = 1.25, CV = c(0.46, 0.46),
design = c("2x2x4"), regulator = c("FDA"),
nsims = 1e+05, nstart = 24, imax=100,
print = TRUE, details = TRUE, setseed=TRUE).
---------------------------------------------
Study design: 2x2x3 (3 period full replicate)
log-transformed data (multiplicative model)
1e+05 studies for each step simulated.
alpha  = 0.05, target power = 0.9
CVw(T) = 0.46; CVw(R) = 0.46
True ratio = 0.93
ABE limits / PE constraints = 0.8 ... 1.25
FDA regulatory settings
- CVswitch            = 0.3
- regulatory constant = 0.8925742
- pe constraint applied
Sample size search
 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
 
Helmut
★★★
avatar
Homepage
Vienna, Austria,
2022-09-16 15:37
(580 d 17:09 ago)

@ Achievwin
Posting: # 23301
Views: 2,366
 

 Script for ABE and SABE (any applicable design)

Post reply
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 [image]-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)
}

—
Dif-tor heh smusma 🖖🏼 Довге життя Україна! [image]
Helmut Schütz
[image]

The quality of responses received is directly proportional to the quality of the question asked. 🚮
Science Quotes
 
Helmut
★★★
avatar
Homepage
Vienna, Austria,
2022-10-24 18:51
(542 d 13:54 ago)

@ Achievwin
Posting: # 23346
Views: 2,104
 

 New article / wacky script

Post reply
Hi Achievwin,

you inspired me to write an article.
Cave: It contains a script with 604 (‼) lines of code. ;-)
Every method implemented in PowerTOST for estimating sample sizes for equivalence (ABE, ABEL, RSABE, RSABE for NTIDs) is covered.
—
Dif-tor heh smusma 🖖🏼 Довге життя Україна! [image]
Helmut Schütz
[image]

The quality of responses received is directly proportional to the quality of the question asked. 🚮
Science Quotes
 
Achievwin
★    

US,
2022-11-11 13:51
(524 d 17:54 ago)

@ Helmut
Posting: # 23360
Views: 1,903
 

 New article / wacky script

Post reply
Helmut:

As always you are a source a for folks across the globe....

Regards,

❝ you inspired me to write an article.

❝ Cave: It contains a script with 604 (‼) lines of code. ;-)

 
UA Flag
Activity
 Admin contact
22,987 posts in 4,824 threads, 1,664 registered users;
82 visitors (0 registered, 82 guests [including 3 identified bots]).
Forum time: 08:46 CEST (Europe/Vienna)

The only way to comprehend what mathematicians mean by Infinity
is to contemplate the extent of human stupidity.    Voltaire

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