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

posted by Helmut Homepage – Vienna, Austria, 2020-11-26 00:09 (146 d 12:17 ago) – Posting: # 22086
Views: 2,104

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 🖖
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,419 posts in 4,475 threads, 1,510 registered users;
online 9 (0 registered, 9 guests [including 2 identified bots]).
Forum time: Wednesday 13:26 CEST (Europe/Vienna)

In the Middles Ages the lingua franca of science was Latin.
Nowadays the language of science is bad English.    Anonymous

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