Bioequivalence and Bioavailability Forum

Dear all,

if you want to use iteratively adjusted α, new R-code and examples at the end of the post. Substantially faster than my pre­vi­ous clumsy attempts. Execution takes generally less than ten seconds on my machine. Results for GMR 0.9, target power 80%, and the most critical CV_WR (30%):

     design  reg      n   TIE1    pwr1  adj. α   %CI    TIE2    pwr2
────────────────────────────────────────────────────────────────────
  RTRT|TRTR  EMA     34  0.0816  80.28  0.0285  94.30  0.0499  72.48
             ANVISA  40  0.0821  85.25  0.0282  94.36  0.0499  78.37
             FDA     32  0.1471  81.67  0.0114  97.72  0.0497  60.50
    RTR|TRT  EMA     50  0.0877  80.16  0.0257  94.86  0.0500  70.84
             ANVISA  58  0.0882  84.56  0.0255  94.90  0.0499  76.41
             FDA     46  0.1460  80.64  0.0112  97.76  0.0500  59.09
RRT|RTR|TRR  EMA     54  0.0719  81.59  0.0338  93.24  0.0499  76.41
             ANVISA  60  0.0719  84.87  0.0339  93.22  0.0500  80.22
             FDA     45  0.1434  80.34  0.0120  97.60  0.0498  59.91

n   : sample size for ≥ target power
TIE1: empiric Type I Error for α 0.05
pwr1: power for α 0.05
TIE2: TIE for the adjusted α
pwr2: power for the adjusted α

The nasty thing is that even if we don’t scale the TIE will be inflated. With

ad.alpha(reg="FDA", des="2x2x4", CV=0.28, n=32, GMR=0.9)

I got

Regulator: FDA (method: RSABE)
CVwr 28%, n 32, design 2x2x4 (RTRT|TRTR)
empiric TIE for α 0.05: 0.1013 (rel. change of risk: +103%)
power for GMR 0.9: 81.39%
iteratively adjusted α: 0.0209 (95.82% CI)
empiric TIE: 0.0499
power for GMR 0.9: 68.81% (rel. loss: -15.5%)

For FDA’s method and CVs close to 30% not even Bon­fer­roni would be enough (i.e., we have to go below 0.025). No ad­just­ment is required for FDA’s method if the CV is lower than ~22% or >30%.
The adjustment range for EMA and ANVISA is ~25% to ~44%.

Is it important for the patient’s risk? I think so. Sponsors seemingly don’t care. A common reply to my con­cerns was “Since it is neither stated in the GL nor the Q&A we keep it as it is.”
Does it have an impact on power? Depends. For low to moderate CVs (EMA, ANVISA) or low CVs (FDA) either one looses power if the study was planned for α 0.05 or – if adjusted alphas are intended to use – one has to increase the sample size ac­cord­ing­ly. For “true” HVDs/HVDPs it doesn’t matter at all.

The question arises whether regulators would accept a method which adapts α in face of the data. Why not? The adjusted α is always ≤0.05. In the context of TSDs the EMA states “[…] using an adjusted co­verage probability which will be higher than 90%”. The devil is in the details. Further down we read “pre­specified in the protocol along with the adjusted significance levels”. This statement effectively pre­vents adaptive methods (where α for the final analysis is also adjusted) from entering the scene.
If you care about the patient’s risk (as I hope) but want to pre-specify an α (see my suggestions for the EMA above), explore the argument alpha of the function. Will always be more conservative than nec­es­sary. In order not to adjust (where no adjustment is required; thus maintaining power), state in the pro­to­col some­thing like “In order to prevent an inflation of the Type I Error (patient’s risk)¹ the following pro­ce­dure will be employed: If CV_WR is in the range 25–44%, a con­ser­va­tive α of 0.025 (95% CI) will be used. Other­wise, the conventional 0.05 (90% CI) will be used.”

Example: GMR 0.9 (planned for 80% power at CV 30%). Comparison of power (iteratively adjusted and pre-specified α)

          RTRT|TRTR (n=34)              RRT|RTR|TRR (n=54)     
    ────────────────────────────   ────────────────────────────
CV%   α     power    α     power     α     power    α     power
───────────────────────────────────────────────────────────────
25  0.0487  86.96  0.0250  79.29   0.0494  88.54  0.0302  83.52
30  0.0285  72.48  0.0250  70.54   0.0338  76.41  0.0302  74.82
35  0.0362  77.24  0.0250  72.23   0.0435  82.29  0.0302  77.70
40  0.0404  82.43  0.0250  76.98   0.0449  86.94  0.0302  82.94
43  0.0472  85.85  0.0250  79.22   0.0484  89.48  0.0302  85.20

I would prefer the former and happily defend it by the oracle’s pre­vious own wis­dom²

“Only statistical pro­cedures, which do not exceed the nominal risk of 5% can be accepted,
and among them the one with the smallest risk of erro­ne­ously rejecting bio­equi­va­lence
should be selected.”

This beautiful sentence tells me two things:

1. Don’t use the method as stated in the GLs. Type I Error might be >5%. Adjust the α.
   
2. Don’t use a pre-specified α. Iteratively adjusted α is more powerful (lower Type II Error).

4-period full replicate, GMR 0.9, n 18–48, CV 20–50%. Left panels unadjusted, right panels adjusted.
Contours enclose combinations of n/CV which are expected to lead to a significant inflation of the TIE.

EMA



FDA

Maximum inflation of the TIE (at CV 30%) for the EMA’s method 0.0825 and for the FDA’s 0.1708. Lower inflation for EMA’s, but IMHO unacceptable (relative increase of the patient’s risk up to 65%). Nasty for the FDA’s if no scaling is applied (risk ≤+241%). At higher CVs the con­ser­va­tism of the test cuts in (TIE drops below nominal α). Not surprising; we see a similar behavior in TOST. For ABEL (dependent on n) this “border” is at CV 41–44%. Due to the discontinuity of RSABE the border is at 30%. Another observation: The sample size has limited in­flu­ence on ABEL’s TIE whereas RSABE’s increases substantially with n.

 n   ABEL    RSABE
18 0.07951 0.12279
48 0.08244 0.17076

---

1. Wonnemann M, Frömke C, Koch A. Inflation of the Type I Error: Investigations on Regulatory Recommendations for Bioequivalence of Highly Variable Drugs. Pharm Res. 2015;32(1):135–43.
   doi 10.1007/s11095-014-1450-z.
   
2. Commission of the European Communities (1991)
   CPMP guideline: Investigation of Bioavailability and Bioequivalence.
   online

---

R-code (with a fair degree of error-handling):

ad.alpha <-function(reg = "EMA", des = "2x3x3", CV = 0.3, n = 54,
                    print = TRUE, details = FALSE,
                    GMR = NA, alpha = 0.05) {
  ## Arguments:
  ##   reg     Regulator: "EMA", "FDA", "ANVISA"
  ##   des     design ("2x2x4", "2x2x3", "2x3x3")
  ##   CV      intra-subject of the reference obtained in a
  ##           replicate design (ratio, /not/ percent)
  ##   n       total sample size or a vector of subjects / sequences
  ##   print   Boolean (FALSE returns a list of results)
  ##   details Boolean (runtime, number of simulations)
  ##   GMR     if given, power is estimated
  ##   alpha   nominal level. Might be useful in assessing pre-
  ##           specified values. Experimental!
  ## Returns:
  ##   alpha.adj  Iteratively adjusted alpha which does not inflate the
  ##              TIE. The computed value is down-rounded at the 4th
  ##              decimal for additional conservatism.
  ##   CI.adj     The adjusted confidence interval in percent, where
  ##              CI.adj = 100*(1-  2*alpha.adj)
  ##   TIE.unadj  The empiric Type I Error based on nominal alpha 0.05.
  ##   TIE.adj    TIE based on adjusted alpha.
  ##   rel.change Relative change in risk (%) compared to nominal 0.05.
  ##   If details = TRUE:
  ##              Runtime in seconds & number of simulations performed.
  ##   If GMR is given:
  ##   pwr.unadj  Power for alpha 0.05.
  ##   pwr.adj    Power for adjusted alpha.
  ##   rel.loss   Relative loss in power if the sample size was planned
  ##              for alpha 0.05 and will be evaluated with alpha.adj,
  ##              where rel.loss = 100*(pwr.adj - pwr.unadj)/pwr.unadj
  ##   If alpha is given:
  ##   Assessment of TIE; alpha is justified if n.s. > 0.05
  ######################################################################
  if (missing(reg)) warning("Default reg='EMA' used.")
  regs <- c("EMA", "FDA", "ANVISA")
  if (reg %in% regs != TRUE)
    stop("reg must be any of 'EMA', 'FDA', or 'ANVISA'.")
  if (missing(des)) warning("Default des='2x3x3' used.")
  design <- c("2x2x4", "2x2x3", "2x3x3")
  if (des %in% design != TRUE)
    stop("des must be any of '2x2x4', '2x2x3', or '2x3x3'.")
  if (missing(CV)) warning("Default CV=0.3 used.")
  if (length(CV) > 1) {
    warning("Vectorized CV not implemented. CV[2] used.")
    CV <- tail(CV, 1)
  }
  if (CV <= 0) stop("Positive CV must be given.")
  if (missing(n)) warning("Default n=54 used.")
  if (sum(n) < 4) stop("Sample size too low.")
  if (((des == "2x2x4" || des == "2x2x3") &&
    (n %% 2 != 0) && length(n) == 1) ||
    (des == "2x3x3" && n %% 3 != 0 && length(n) == 1))
    warning("Unbalanced sequences. Give n as a vector.")
  if (!is.na(GMR) && (GMR < 0.8 || GMR > 1.25))
    stop("GMR must be within 0.80-1.25.")
  n   <- as.integer(n)
  pwr <- rep(NA, 2)
  TIE <- rep(NA, 2)
  adj <- NA
  require(PowerTOST) # what else?
  fun1 <- function(x) power.scABEL(alpha=x, CV=CV, theta0=theta0,
                                   n=n, design=des, nsims=1e6)-0.05
  fun2 <- function(x) power.RSABE(alpha=x, CV=CV, theta0=theta0,
                                  n=n, design=des, nsims=1e6)-0.05
  type <- c("RTRT|TRTR", "RTR|TRT", "RRT|RTR|TRR") # clear words
  sig  <- binom.test(x=0.05*1e6, n=1e6, alternative="less",
                     conf.level=1-0.05)$conf.int[2] # binom test limit
  if (details) ptm <- proc.time()
  theta0 <- scABEL(CV=CV, regulator=reg)[["upper"]]
  if (reg != "FDA") { # EMA or ANVISA
    TIE[1] <- power.scABEL(alpha=alpha, CV=CV, theta0=theta0,
                           n=n, design=des, nsims=1e6)
    if (!is.na(GMR)) { # only if GMR given
      pwr[1] <- power.scABEL(alpha=alpha, CV=CV, theta0=GMR,
                             n=n, design=des)
    }
    if (TIE[1] > sig) { # adjust only if needed
      x      <- uniroot(fun1, interval=c(0.01, 0.05), tol=1e-8)
      adj    <- floor(x$root*1e4)/1e4 # round down
      TIE[2] <- power.scABEL(alpha=adj, CV=CV, theta0=theta0,
                             n=n, design=des, nsims=1e6)
      if (!is.na(GMR)) { # power for adj. alpha
        pwr[2] <- power.scABEL(alpha=adj, CV=CV, theta0=GMR,
                               n=n, design=des)
      }
    }
  } else { # FDA
    TIE[1] <- power.RSABE(alpha=alpha, CV=CV, theta0=theta0,
                          n=n, design=des, nsims=1e6)
    if (!is.na(GMR)) { # only if GMR given
      pwr[1] <- power.RSABE(alpha=alpha, CV=CV, theta0=GMR,
                            n=n, design=des)
    }
    if (TIE[1] > sig) { # adjust only if needed
      x      <- uniroot(fun2, interval=c(0.005, 0.05), tol=1e-8)
      adj    <- floor(x$root*1e4)/1e4 # round down
      TIE[2] <- power.RSABE(alpha=adj, CV=CV, theta0=theta0,
                            n=n, design=des, nsims=1e6)
      if (!is.na(GMR)) { # power for adj. alpha
        pwr[2] <- power.RSABE(alpha=adj, CV=CV, theta0=GMR,
                              n=n, design=des)
      }
    }
  }
  if (details) run.time <- proc.time()-ptm
  no  <- 1e6*(1 + x$iter)
  if (!is.na(GMR)) no <- no + 2e5
  if (print) { # fetch and print results
    txt <- paste0("\nRegulator: ", reg, " (method: ")
    ifelse(reg != "FDA",
      txt <- paste0(txt, "ABEL)\n"),
      txt <- paste0(txt, "RSABE)\n"))
    txt <- paste0(txt, "CVwr ", 100*CV, "%, n ", paste0(n, collapse=", "),
                  ", design ", des, " (",
                   type[match(des, design)], ")\n", "empiric TIE ",
                  "for \u03B1 ", alpha, ": ", sprintf("%.4f", TIE[1]))
    if (TIE[1] > sig || alpha != 0.05) {
      rel.change <- 100*(TIE[1]-0.05)/0.05
      txt <- paste0(txt, " (rel. change of risk: ",
                    sprintf("%+1.3g%%", rel.change), ")")
    }
    if (!is.na(pwr[1]) || alpha != 0.05) {
      pwr.unadj <- pwr[1]*100
      txt <- paste0(txt, "\npower for GMR ", GMR, ": ",
                    sprintf("%.2f%%", pwr.unadj))
    }
    if (TIE[1] > sig) {
      CI.adj <- sprintf("%.2f%%", 100*(1-2*adj))
      txt <- paste0(txt, "\niteratively adjusted \u03B1: ",
                    sprintf("%.4f", adj), " (", CI.adj,
                    " CI)\nempiric TIE: ", sprintf("%.4f", TIE[2]))
    if (!is.na(pwr[2])) {
      pwr.adj <- pwr[2]*100
      txt <- paste0(txt, "\npower for GMR ", GMR, ": ",
                    sprintf("%.2f%%", pwr.adj), " (rel. loss: ",
                    sprintf("%+1.3g%%", 100*(pwr[2]-pwr[1])/pwr[1]),
                    ")\n\n")
    } else {
      txt <- paste0(txt, "\n\n")
    }
      if (details) {
        txt <- paste0(txt, "Runtime              : ",
                      signif(run.time[3], 3), " seconds",
                      "\nNumber of simulations: ",
                      formatC(no, format="d", big.mark=",",
                        decimal.mark="."), "\n\n")
      }
    } else {
      txt <- paste0(txt, "\nNo significant inflation of the TIE ",
                         "expected;")
      ifelse(alpha == 0.05,
        txt <- paste0(txt, "\nno adjustment of \u03B1 is required.\n\n"),
        txt <- paste0(txt, "\nthe chosen pre-specified \u03B1 ",
                           "is justified.\n\n"))
    }
  cat(txt)
  } else { # prepare and return list of results
    res <- list(regulator=reg, method=ifelse(reg!="FDA", "ABEL", "RSABE"),
                design=des, type=type[match(des, design)], alpha=alpha,
                CV=CV, n=n, GMR=GMR, TIE.unadj=TIE[1],
                rel.change=ifelse(!is.na(adj), 100*(TIE[1]-0.05)/0.05, NA),
                pwr.unadj=pwr[1]*100, alpha.adj=adj,
                CI.adj=ifelse(!is.na(adj), 100*(1-2*adj), NA),
                TIE.adj=TIE[2], pwr.adj=pwr[2]*100,
                rel.loss=ifelse(!is.na(pwr[2]),
                  100*(pwr[2]-pwr[1])/pwr[1], NA), sims=no)
    return(res)
  }
}

Examples:
ad.alpha(reg="EMA", des="2x2x4", CV=0.3, n=34)
ad.alpha(reg="FDA", des="2x3x3", CV=0.3, n=45, details=TRUE)
ad.alpha(reg="FDA", des="2x3x3", CV=c(0.25, 0.3), n=45)
ad.alpha(CV=0.3, n=54)            # some of the defaults
ad.alpha(CV=0.3, n=52)            # unbalanced; tries to guess…
ad.alpha(CV=0.3, n=c(18, 18, 16)) # better!
ad.alpha(CV=0.3, n=51)            # unbalanced; code doesn’t “know” that (assumes 17|17|17)
ad.alpha(CV=0.3, n=c(18, 17, 16)) # only you know!

If you are interested in only some of the results, use print=FALSE and assign to a variable. The function returns a list containing regulator, method, design, type, alpha, CV, n, TIE.unadj, rel.change, pwr.unadj, alpha.adj, CI.adj, TIE.adj, pwr.adj, and rel.loss.
Examples:

x <- ad.alpha(reg="EMA", des="2x2x4", CV=0.3, n=34, print=FALSE)
x$TIE.unadj
[1] 0.081626
x$rel.change
[1] 63.252
x$alpha.adj
[1] 0.0285
x$TIE.adj
[1] 0.04989
unlist(x[12:14])
alpha.adj    CI.adj   TIE.adj
 "0.0285"  "94.30%" "0.04989"

If you give the optional argument GMR, you can explore the impact on power by the adjustment:

ad.alpha(reg="EMA", des="2x2x4", CV=0.3, n=34, GMR=0.9)

Regulator: EMA (method: ABEL)
CVwr 30%, n 34, design 2x2x4 (RTRT|TRTR)
empiric TIE for α 0.05: 0.0816 (rel. change of risk: +63.3%)
power for GMR 0.9: 80.28%
iteratively adjusted α: 0.0285 (94.30% CI)
empiric TIE: 0.0499
power for GMR 0.9: 72.48% (rel. loss: -9.72%)

Pre-specified lower α:

ad.alpha(reg="EMA", des="2x2x4", CV=0.3, n=34, alpha=0.025)

Regulator: EMA (method: ABEL)
CVwr 30%, n 34, design 2x2x4 (RTRT|TRTR)
empiric TIE for α 0.025: 0.0444 (rel. change of risk: -11.1%)
No significant inflation of the TIE expected;
the chosen pre-specified α is justified.

ad.alpha(reg="EMA", des="2x3x3", CV=0.3, n=54, alpha=0.0302)

Regulator: EMA (method: ABEL)
CVwr 30%, n 54, design 2x3x3 (RRT|RTR|TRR)
empiric TIE for α 0.0302: 0.0448 (rel. change of risk: -10.4%)
No significant inflation of the TIE expected;
the chosen pre-specified α is justified.

The impact on power might be severe for the FDA’s method:

ad.alpha(reg="FDA", des="2x2x4", CV=0.3, n=32, print=FALSE, GMR=0.9)$pwr.adj
[1] "60.505"

Fancy:

x <- ad.alpha(reg="EMA", des="2x2x4", CV=0.3, n=34, print=FALSE, GMR=0.9)
y <- matrix(data="", nrow=2, ncol=12, byrow=TRUE, dimnames=list(NULL,
              c("regulator", "method", "design", "type", "CV", "n", "GMR",
                "alpha", "TIE", "rel.change", "power", "rel.loss")))
y[1, 1:4]   <- as.character(x[1:4])
y[1, 5]     <- as.numeric(x[6])
y[1, 6]     <- as.numeric(sum(unlist(x[7])))
y[1, 7]     <- as.numeric(x[8])
y[1, 8]     <- sprintf("%.4f", x[5])
y[1, 9]     <- sprintf("%.5f", x[9])
y[1, 10]    <- sprintf("%+1.3g%%", x[10])
y[1, 11]    <- sprintf("%.2f%%", x[11])
y[2, 8]     <- sprintf("%.4f", x[12])
y[2, 9]     <- sprintf("%.5f", x[14])
y[2, 11:12] <- sprintf("%.2f%%", x[15:16])
print(as.data.frame(y), row.names=FALSE)

regulator method design      type  CV  n GMR  alpha     TIE rel.change  power rel.loss
      EMA   ABEL  2x2x4 RTRT|TRTR 0.3 34 0.9 0.0500 0.08163     +63.3% 80.28%         
                                             0.0285 0.04989            72.48%   -9.72%