Bioequivalence and Bioavailability Forum

Dear Ohlbe,

❝ ❝ Didn’t know that the documents are in the public domain in the meantime.

❝ I found them by chance, linked to this paper and that one from the same author.

THX, interesting.

❝ ❝ However, as I wrote above it’s a pretty trivial exercise for anybody to critically inspect the data.

❝ "Trivial" if you have the coding skills or know someone who does, yes, definitely.

Your wish is my command. script at the end. To play with the script, download my data. There is another one, where I multiplied all original values of T up to subject 60 with 0.85 and then by 1.30.

Is this manipulation really obvious? By chance the T/R-ratios and PEs start to increase with subject 40. Would the original data already trigger an alarm?*

❝ And yes, it definitely makes sense, even (especially ?) if you are buying a dossier from the shelf.

❝

❝ ❝ Possibly the statistician was fired by the man in the Armani suit shouting

❝ ❝ I calculated post hoc power in FARTSSIE and it’s already 81.6%! What’s the point of having 99.4% at the end?«

❝

❝ Come on, the man in the Armani suit has no idea how to use FARTSSIE and has even never heard of it

How could I forget? Only professional in PowerPoint, Twitter, WhatsApp, copypasting from one document to another, networking at conferences, smalltalk at cocktail parties. Preferred communication via LinkedIn.

---

• Fuglsang A. Detection of data manipulation in bioequivalence trials. Eur J Pharm Sci. 2021; 156: 105595. doi:10.1016/j.ejps.2020.105595.

---

Currently only for 2×2×2 crossover designs. If you are not interested in the cumulative PE and CI, delete everything in red. Then you need only the columns Subject, Treatment, Cmax in the file. If you want a version reading from XLS(X), let me know.

library(PowerTOST)
path <- "your path"
file <- "yourdata.csv"
# required case-sensitive column names (in any order):
#   Subject, Sequence, Treatment, Period, Cmax
#   the treatments must be coded as T and R
sep  <- "," # column separator: inspect the file; alternatively "\t" or ";"
dec  <- "." # decimal separator; alternatively ","
com  <- "#" # for eventual comments, sometimes found in the header
na   <- c("NA", "ND", ".", "", "Missing") # will be forced to NA
if (sep == dec) stop("sep and dec must be different!")
data <- read.csv(file = paste0(path, "/", file), sep = sep, dec = dec, quote = "",
                 strip.white = TRUE, comment.char = com, na.strings = na)
# if there are other columns in the file, drop them
keep <- c("Subject", "Sequence", "Treatment", "Period", "Cmax")
data <- data[, keep]
subjects <- unique(data$Subject)
n        <- length(subjects)
Cmax.ratio <- data$Cmax[data$Treatment == "T"] / data$Cmax[data$Treatment == "R"]
Cmax.cum   <- numeric(n)
for (i in 1:n) {
  Cmax.cum[i] <- prod(Cmax.ratio[1:i])^(1/i)
}
res  <- data.frame(PE = rep(NA_real_, n-1), lower = NA_real_,
                   upper = NA_real_, CV = NA_real_, power = NA_real_)
# factorize for the linear model
facs <- c("Subject", "Sequence", "Treatment", "Period")
data[facs] <- lapply(data[facs], factor)
ow   <- options("digits")
options(digits = 12) # more digits for anova
on.exit(ow)          # ensure that options are reset if an error occurs
for (i in 3:n) {
  tmp            <- head(data, i * 2)
  m              <- lm(log(Cmax) ~ Subject + Period + Sequence + Treatment,
                                   data = tmp)
  res$PE[i-1]    <- exp(coef(m)[["TreatmentT"]])
  res[i-1, 2:3]  <- as.numeric(exp(confint(m, "TreatmentT",
                                           level = 1 - 2 * 0.05)))
  res$CV[i-1]    <-mse2CV(anova(m)[["Residuals", "Mean Sq"]])
  res$power[i-1] <- suppressMessages(
                      power.TOST(CV = res$CV[i-1], theta0 = res$PE[i-1],
                                 n = length(unique(tmp$Subject))))
}
# print(round(100 * res, 2)) # uncomment this line for PE, CI, CV, post hoc power
dev.new(width = 4.5 * sqrt(2), height = 4.5)
op   <- par(no.readonly = TRUE)
par(mar = c(3.2, 3.1, 0, 0.1), cex.axis = 0.9,
    mgp = c(2, 0.5, 0), ljoin = "mitre", lend = "square")
plot(1:n, Cmax.ratio, type = "n", log = "y", axes = FALSE, xlab = "analysis",
     ylab = expression(italic(C)[max[T]]*" / "*italic(C)[max[R]]))
x.tick <- pretty(1:n)
x.tick <- c(1, x.tick[2:(length(x.tick)-2)], n)
axis(1, at = x.tick, labels = paste0("#", x.tick))
y.tick <- sort(c(0.8, 1.25, axTicks(side = 2, log = TRUE)))
axis(2, at = y.tick, las = 1)
grid(nx = NA, ny = NULL)
abline(v = x.tick, lty = 3, col = "lightgrey")
abline(h = c(0.8, 1, 1.25), lty = c(2, 1, 2))
box()
legend("topleft", legend = c("individual T/R-ratio", "cumulative T/R-ratio",
                             "point estimate", "confidence limits"),
       pch = c(21, NA, NA, NA), lty = c(0, 1, 1, 1), lwd = c(1, 2, 2),
       col = c(rep("blue", 2), "darkred", "red"), pt.bg = "#AAAAFF80",
       pt.cex = 1, bg = "white", seg.len = 3, cex = 0.9, y.intersp = 1.1)
lines(1:n, Cmax.cum, col = "blue", lwd = 2, type = "s")
lines(2:n, res$PE, col = "darkred", lwd = 2, type = "s")
lines(2:n, res$lower, col = "red", type = "s")
lines(2:n, res$upper, col = "red", type = "s")
points(1:n, Cmax.ratio, pch = 21, col = "blue", bg = "#AAAAFF80")
par(op)