Bioequivalence and Bioavailability Forum

Dear Detlew,

❝ if the CI inclusion rule is really something different than TOST why do we calculate power / samplesize based on TOST .

Cause you decided to baptize the package PowerTOST and not “Eierlegende Wollmilchsau” back in 2009.

Given, Donald used the phrase “operationally identical” on p.661 (right column, 2^nd paragraph).

However, for me (!) those are two different “operations”. Results of an example:

theta1 = 0.8000
theta2 = 1.2500
alpha  = 0.05
90% CI: lower CL = 0.8448
        upper CL = 1.1003
        CI within 0.8000 and 1.2500: passed BE
TOST  : p(<0.8000) = 0.01239
        p(>1.2500) = 0.001565
        p(<0.8000) <0.05 and p(>1.2500) <alpha: passed BE

Checks with library PowerTOST
CI.BE(alpha = 0.05, ...):
  lower   = 0.8448
  upper   = 1.1003
pvalues.TOST():
  p.left  = 0.01239
  p.right = 0.001565

---

set.seed(123456789)
alpha  <- 0.05
theta0 <- 0.95
theta1 <- 0.80
theta2 <- 1.25
CV     <- 0.20
subj   <- 1:20
per    <- 1:2
seq    <- c("TR", "RT")
trt    <- c("T", "R")
data   <- data.frame(subj = as.factor(rep(subj, each = 2)),
                     per  = as.factor(rep(per, 2)),
                     seq  = as.factor(rep(seq, each = length(subj))),
                     trt  = as.factor(c(rep(trt, length(subj)/2),
                                        rep(rev(trt), length(subj)/2))))
data$PK[data$trt == "R"] <- rlnorm(n = length(subj),
                                   meanlog = log(1)-0.5*log(CV^2+1),
                                   sdlog = sqrt(log(CV^2+1)))
data$PK[data$trt == "T"] <- rlnorm(n = length(subj),
                                   meanlog = log(theta0)-0.5*log(CV^2+1),
                                   sdlog = sqrt(log(CV^2+1)))
mod    <- lm(log(PK) ~ seq + subj + per + trt, data = data)
#########################
# CI inclusion (base R) #
#########################
pe     <- exp(coef(mod)[["trtT"]])
ci     <- exp(as.numeric(confint(mod, "trtT", level = 1-2*alpha)))
names(ci) <- c("lower", "upper")
if (ci[["lower"]] >= 0.80 & ci[["upper"]] <=1.25) {
  BE.ci <- sprintf("CI within %.4f and %.4f: %s",
                   theta1, theta2, "passed BE")
}
if (ci[["lower"]] < 0.80 & ci[["upper"]] <= 1.25) {
  BE.ci <- sprintf("lower CL <%.4f, upper CL \u2264%.4f: %s",
                   theta1, theta2, "failed BE (inconclusive)")
}
if (ci[["lower"]] >= 0.80 & ci[["upper"]] > 1.25) {
  BE.ci <- sprintf("lower CL \u2265%.4f, upper CL >%.4f: %s",
                   theta1, theta2, "failed BE (inconclusive)")
}
if (ci[["lower"]] < 0.80 & ci[["upper"]] > 1.25) {
  BE.ci <- sprintf("lower CL <%.4f and upper CL >%.4f: %s",
                   theta1, theta2, "inequivalent")
}
#################
# TOST (base R) #
#################
s      <- sqrt(anova(mod)["Residuals", "Mean Sq"])
nu     <- anova(mod)["Residuals", "Df"]
n1     <- length(as.numeric(data$subj[data$seq == "TR"]))/2
n2     <- length(as.numeric(data$subj[data$seq == "RT"]))/2
muR    <- mean(log(data$PK[data$trt == "R"]))
muT    <- mean(log(data$PK[data$trt == "T"]))
delta  <- muT-muR
se     <- s*sqrt(0.5*(1/n2+1/n2))
t1     <- (delta-log(theta1))/se
t2     <- (log(theta2)-delta)/se
p.TOST <- c(p1 = pt(t1, df = nu, lower.tail = FALSE),
            p2 = pt(t2, df = nu, lower.tail = FALSE))
if (max(p.TOST) < alpha) {
  BE.TOST <- sprintf("p(<%.4f) <%.2g and p(>%.4f) <%.2g: %s",
                     theta1, alpha, theta2, alpha, "passed BE")
}
if (p.TOST[["p1"]] < alpha & p.TOST[["p2"]] >= alpha) {
  BE.TOST <- sprintf("p(<%.4f) <%.2g and p(>%.4f) \u2265%.2g: %s",
                     theta1, alpha, theta2, alpha, "failed BE (inconclusive)")
}
if (p.TOST[["p1"]] >= alpha & p.TOST[["p2"]] < alpha) {
  BE.TOST <- sprintf("p(<%.4f) \u2265%.2g and p(>%.4f) <%.2g: %s",
                     theta1, alpha, theta2, alpha, "failed BE (inconclusive)")
}
if (p.TOST[["p1"]] >= alpha & p.TOST[["p2"]] >= alpha) {
  BE.TOST <- sprintf("p(<%.4f) \u2265%.2g and p(>%.4f) \u2265%.2g: %s",
                     theta1, alpha, theta2, alpha, "inequivalent")
}
#####################
# library PowerTOST #
#####################
library(PowerTOST)
pt.ci <- CI.BE(alpha = alpha, pe = pe, CV = se2CV(s), n = c(n1, n2))
pt.p  <- pvalues.TOST(pe = pe, CV = se2CV(s), n = c(n1, n2),
                      theta1 = theta1, theta2 = theta2)
cat(paste0("\ntheta1 = ", sprintf("%.4f", theta1),
           "\ntheta2 = ", sprintf("%.4f", theta2),
           "\nalpha  = ", alpha, sprintf("\n%g%% CI: ", 100*(1-2*alpha)),
           "lower CL = ", sprintf("%.4f", ci["lower"]),
           "\n        upper CL = ", sprintf("%.4f", ci["upper"]),
           "\n        ", BE.ci,
           "\nTOST  : p(<", sprintf("%.4f) = %.4g", theta1, p.TOST[["p1"]]),
           "\n        p(>", sprintf("%.4f) = %.4g", theta2, p.TOST[["p2"]]),
           "\n        ", BE.TOST,
           "\n\nChecks with library PowerTOST",
           "\nCI.BE(alpha = ", alpha, ", ...):",
           "\n  lower   = ", sprintf("%.4f", pt.ci["lower"]),
           "\n  upper   = ", sprintf("%.4f", pt.ci["upper"]),
           "\npvalues.TOST():",
           "\n  p.left  = ", sprintf("%.4g", pt.p[["p.left"]]),
           "\n  p.right = ", sprintf("%.4g", pt.p[["p.right"]])), "\n")