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

posted by Helmut Homepage – Vienna, Austria, 2022-12-25 13:17 (429 d 14:28 ago) – Posting: # 23411
Views: 1,630

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,912 posts in 4,806 threads, 1,636 registered users;
26 visitors (0 registered, 26 guests [including 7 identified bots]).
Forum time: 03:45 CET (Europe/Vienna)

Exploratory analysis: The art of finding a Rembrandt
in a Jackson Pollock.    Stephen Senn

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