T/R-ratios are useful [GxP / QC / QA]

posted by Helmut Homepage – Vienna, Austria, 2022-05-25 15:35 (81 d 18:29 ago) – Posting: # 23021
Views: 1,124

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. [image] 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.

[image]

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 :-D

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.



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)


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
22,301 posts in 4,667 threads, 1,585 registered users;
online 3 (0 registered, 3 guests [including 3 identified bots]).
Forum time: Monday 10:04 CEST (Europe/Vienna)

Science… never solves a problem
without creating ten more.    George Bernard Shaw

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