Models, models! [General Sta­tis­tics]

posted by Helmut Homepage – Vienna, Austria, 2022-12-25 13:17 (44 d 12:22 ago) – Posting: # 23411
Views: 607

Hi Brus,

❝ Can anyone confirm that it is correct?


That’s an approximation (like a model with subject as the only effect).
According to the EMA’s Q&A document we have to fit a linear model of R-data with fixed effects sequence, subject(sequence), period and calculate CVwR from the models’ residual MSE. As usual, you could drop sequence (keep only subject and period). :-D

library(replicateBE)
library(PowerTOST)
df       <- data.frame(method = c("approximate", "simple", "EMA",
                                  "no sequence", "FDA"),
                       CVwR = NA_real_)
# EMA reference data set I, only untransformed PK after R
dataR    <- rds01[rds01$treatment == "R", c(1:3, 5)]
subj     <- levels(dataR$subject)
varR     <- rep(NA_real_, length(subj))
dlat     <- data.frame(sequence = NA, cont = rep(NA_real_, length(subj)))
for (j in seq_along(subj)) {
  varR[j] <- var(log(dataR[dataR$subject == subj[j], "PK"]))
  if (length(dataR[dataR$subject == subj[j], "PK"]) == 2) {
    # intra-subject contrast
    dlat$cont[j]     <- diff(log(dataR[dataR$subject == subj[j], "PK"]))
    dlat$sequence[j] <- as.character(unique(dataR[dataR$subject == subj[j], "sequence"]))
  }
}
dlat$sequence <- factor(dlat$sequence)
# All analysis on complete cases only!
# Approximation

df[1, 2] <- 100 * mse2CV(mean(varR, na.rm = TRUE))
# Simple model
mod1     <- lm(log(PK) ~ subject, data = dataR, na.action = na.omit)
df[2, 2] <- 100 * mse2CV(anova(mod1)["Residuals", "Mean Sq"])
# Model acc. to the EMA’s Q&A document
mod2     <- lm(log(PK) ~ sequence + subject %in% sequence + period,
                         data = dataR, na.action = na.omit)
df[3, 2] <- 100 * mse2CV(anova(mod2)["Residuals", "Mean Sq"])
# Sequence effects dropped
mod3     <- lm(log(PK) ~ subject + period, data = dataR, na.action = na.omit)
df[4, 2] <- 100 * mse2CV(anova(mod3)["Residuals", "Mean Sq"])
# FDA’s model
mod4     <- lm(cont ~ sequence, data = dlat, na.action = na.omit)
df[5, 2] <- 100 * mse2CV(anova(mod4)["Residuals", "Mean Sq"] / 2)
print(df, digits = 5, row.names = FALSE, right = FALSE)

 method      CVwR 
 approximate 47.627
 simple      47.627
 EMA         46.964
 no sequence 46.964
 FDA         46.964

BTW, from the FDA’s mixed-effects model we could directly get CVwR with 47.328. I love lacking harmonization.

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,485 posts in 4,710 threads, 1,603 registered users;
15 visitors (0 registered, 15 guests [including 10 identified bots]).
Forum time: 01:40 CET (Europe/Vienna)

No matter what side of the argument you are on,
you always find people on your side
that you wish were on the other.    Thomas Berger

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