CI inclusion ≠ TOST [BE/BA News]

posted by Helmut Homepage – Vienna, Austria, 2020-12-06 23:09 (720 d 10:45 ago) – Posting: # 22119
Views: 2,870

Dear Detlew,

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


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, 2nd 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")


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

Complete thread:

UA Flag
Activity
 Admin contact
22,417 posts in 4,692 threads, 1,596 registered users;
21 visitors (0 registered, 21 guests [including 11 identified bots]).
Forum time: 09:54 CET (Europe/Vienna)

To know that we know what we know,
and to know that we do not know what we do not know,
that is true knowledge.    Nicolaus Copernicus

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