CI inclusion ≠ TOST [BE/BA News]

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, 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 🖖
Helmut Schütz The quality of responses received is directly proportional to the quality of the question asked. 🚮
Science Quotes Ing. Helmut Schütz 