## 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 