Inflated type I error with fixed widened limits? [RSABE / ABEL]

posted by Helmut Homepage – Vienna, Austria, 2020-11-26 01:09 (1692 d 07:51 ago) – Posting: # 22086
Views: 11,897

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: […] 106 studies with θ0 = 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

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
23,436 posts in 4,932 threads, 1,677 registered users;
59 visitors (0 registered, 59 guests [including 12 identified bots]).
Forum time: 10:00 CEST (Europe/Vienna)

You can’t fix by analysis
what you bungled by design.    Richard J. Light, Judith D. Singer, John B. Willett

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