CI inclusion ≠ TOST [Regulatives / Guidelines]

posted by Helmut Homepage – Vienna, Austria, 2020-12-06 22:09 (91 d 00:03 ago) – Posting: # 22119
Views: 950

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 🖖
Helmut Schütz
[image]

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

Complete thread:

Activity
 Admin contact
21,371 posts in 4,463 threads, 1,496 registered users;
online 6 (0 registered, 6 guests [including 4 identified bots]).
Forum time: Sunday 22:12 CET (Europe/Vienna)

When people learn no tools of judgment
and merely follow their hopes,
the seeds of political manipulation are sown.    Stephen Jay Gould

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