Bioequivalence and Bioavailability Forum

Dear Detlew,

❝ ❝ […] we have a data-driven decision – which might be false and result in an inflated type I error.

❝

❝ Correct.

Shit.

❝ ❝ I performed simulations acc. to my understanding of the GL: […] 10⁶ studies with θ₀ = 1.25: ~20.6% passed…

❝ Confirmed! At least in magnitude.

❝ My result: 20.46% passed.

❝ Quick and dirty captured code from power.scABEL.sds.

How, did you do that? That’s beyond me.

I misused your old subject sim code and added:

ow     <- options("digits")
nsims  <- 1e6
CVT    <- CVR <- 0.3
sWT    <- CV2se(CVT)
sWR    <- CV2se(CVR)
sBT    <- CV2se(CVT*2)
sBR    <- CV2se(CVR*2)
n      <- c(20, 20) # power ~0.81 for theta0 = 0.9
theta0 <- 1.25      # for TIE
options(digits = 12)
on.exit(options(ow))
set.seed(123456)
mvc    <- mean_vcov(c("TRTR", "RTRT"), muR = log(100), ldiff = log(theta0),
                      sWT = sWT, sWR = sWR, sBT = sBT, sBR = sBR, rho = 1)
res    <- data.frame(CVwR = rep(NA, nsims), PE = NA, lower = NA, upper = NA,
                     L = NA, U = NA, BE = FALSE)
pb     <- txtProgressBar(0, 1, 0, char="\u2588", width=NA, style=3)
st     <- proc.time()[[3]]
for (j in 1:nsims) {
  data        <- prep_data(seqs = c("TRTR", "RTRT"), n = n, metric = "PK",
                           dec = 5, mvc_list = mvc)
  data        <- data[, !(names(data) %in% c("seqno", "logval"))]
  cols        <- c("subject", "period", "sequence", "treatment")
  data[cols]  <- lapply(data[cols], factor)
  modBE       <- lm(log(PK) ~ sequence + subject + period + treatment,
                              data = data)
  res[j, 2:3] <- round(100*as.numeric(exp(confint(modBE, "treatmentT",
                                                  level = 1-2*0.05))), 2)
  res$PE[j]   <- round(100*as.numeric(exp(coef(modBE)[["treatmentT"]])), 2)
  modCVR      <- lm(log(PK) ~ sequence + subject + period,
                              data = data[data$treatment == "R", ])
  res$CVwR[j] <- mse2CV(anova(modCVR)["Residuals", "Mean Sq"])
  ifelse (res$CVwR[j] <= 0.3, res$L[j] <- 80, res$L[j] <- 75)
  res$U[j]    <- 100/(res$L[j]*0.01)
  if (res$lower[j] >= res$L[j] & res$upper[j] <= res$U[j] &
      res$PE[j] >= 80 & res$PE[j] <= 125) res$BE[j] <- TRUE
  setTxtProgressBar(pb, j/nsims)
}
et     <- proc.time()[[3]]
close(pb)
options(ow) # restore options
cat(paste0(j, " subj. simulations with theta0 = ", theta0, ": ",
           signif(sum(res$BE[!is.na(res$L)])/j, 5),
           " passed BE (empiric TIE)\nRuntime: ", signif((et-st)/60, 4),
           " minutes\n"))

And got 20.568% after eleven (!) hours.
PS: Only one of the sim’s failed on the PE restriction.
res[(res$lower >= res$L & res$upper <= res$U) & (res$PE <80 | res$PE >125), ]
            CVwR     PE  lower upper  L        U    BE
566300 0.3000039 125.03 117.36 133.2 75 133.3333 FALSE