fixed or mixed effects model [BE/BA News]
❝ I am confused, if I understand well, we should all treat subject as a fixed effect?
In a 2×2×2 crossover with subject(s) missing in one period you get (apart from a slight difference in the nth figure, which is not relevant) the same result for a mixed effects model with all subject(s) and a fixed effects model excluding the subject(s) with missings.
If you exclude the subject(s) with missings, a mixed effects model will given you essentially the same result than a fixed effects model. Quite often the differences are at the numerical resolution of the software. Since the EMA implemented ICH M13A already (effective 25 January 2025; see this post), you could use any.
However, I would still use a mixed effects model for the FDA and Health Canada and a fixed effects model for all other agencies. We don’t want to confuse assessors who are used to what they required for ages.
There is one advantage of a mixed effects model: Apart from the within-subject subject CV you get additionally the between-subject subject CV. At least nice to know.
A simulation (mixed effects models in are a beast). I give the point estimate and confidence interval in percent with five decimal places for comparison. If rounded to two decimal places, differences for the balanced and imbalanced cases will disappear (note that ICH M13A is silent about rounding).
CV : 0.225
T/R-ratio : 0.95
lower limit : 0.8000
upper limit : 1.2500
power : at least 0.8
alpha : 0.0500
n : 24 (complete data, balanced sequences)
dropout(s) : 1 (both periods)
missing(s) : 1 (second period)
data package function subject df method PE lower upper width
balanced stats lm fixed 22.0 Residual 93.71542 85.34248 102.90983 17.56735
balanced nlme lme random 22.0 Residual 93.71542 85.34247 102.90983 17.56736
balanced lmerTest lmer random 22.0 Satterthwaite 93.71542 85.34248 102.90983 17.56735
balanced lmerTest lmer random 22.0 Kenward-Roger 93.71542 85.34248 102.90983 17.56735
imbalanced stats lm fixed 21.0 Residual 92.96790 84.36634 102.44643 18.08008
imbalanced nlme lme random 21.0 Residual 92.96790 84.36634 102.44644 18.08010
imbalanced lmerTest lmer random 21.0 Satterthwaite 92.96790 84.36634 102.44643 18.08008
imbalanced lmerTest lmer random 21.0 Kenward-Roger 92.96790 84.36634 102.44643 18.08008
incomplete stats lm fixed 20.0 Residual 92.19586 83.31751 102.02028 18.70277
incomplete nlme lme random 20.0 Residual 92.03122 83.19633 101.80431 18.60798
incomplete lmerTest lmer random 20.2 Satterthwaite 92.03122 83.20136 101.79815 18.59679
incomplete lmerTest lmer random 20.1 Kenward-Roger 92.03122 83.19399 101.80718 18.61319
balanced sequences (nRT = nRT)
imbalanced sequences (one subject missing)
incomplete data (as above and the second period of the last subject missing)
df : degrees of freedom =
n – 2 for Residual, approximated by Satterthwaite and Kenward-Roger
method: df method applied
PE : Point Estimate
lower : lower limit of the 90% CI
upper : upper limit of the 90% CI
- For both balanced and imbalanced data the point estimate is identical regardless the model applied.
- For balanced data we get practically the same results for the fixed and mixed effects model (slight difference in the fifth decimal place of the lower confidence limit).
- For imbalanced data, there are slight differences in the fifth decimal place of the upper confidence limit.
- There are differences for incomplete data (period missing) – not only for the CI but also the PE. However, in such a case the subject should be excluded according to M13A and its Q&A anyway.
Residual
(in R
and Phoenix/WinNonlin) is equivalent to DDFM = CONTAIN
and containment
, respectively.Lengthy -script for simulating a 2×2×2 crossover. Defaults: alpha 0.05 (90% CI), BE-limits 0.8–1.25, T/R-ratio 0.95. power ≥ 0.8. Optionally specification of period- and/or carryover-effects.
library(PowerTOST) # for sample size estimation
library(nlme) # for mixed effects model by lm()
suppressMessages(library(lmerTest)) # requires also the packages lme4 and Matrix
# for mixed effects model by lmer() of lme4,
# Satterwaite and Kenward-Roger df
sim.study <- function(alpha = 0.05, theta0, theta1, theta2, targetpower, CV, n, CVb,
per.effect = c(0, 0), carryover = c(0, 0), setseed = TRUE) {
# simulate subject data
if (setseed) set.seed(1234569)
if (missing(theta0)) theta0 <- 0.95
if (missing(theta1) & missing(theta2)) theta1 <- 0.80
if (missing(theta1)) theta1 <- 1 / theta2
if (missing(theta2)) theta2 <- 1 / theta1
if (missing(targetpower)) targetpower <- 0.80
if (missing(CV)) stop("CV must be given.")
if (missing(n)) {
n <- sampleN.TOST(alpha = alpha, CV = CV, theta0 = theta0, theta1 = theta1,
theta2 = theta2, targetpower = targetpower,
print = FALSE)[["Sample size"]]
}
CVb <- CV * 1.5 # arbitrary (cosmetics)
sd <- CV2se(CV)
sd.b <- CV2se(CVb)
# within subjects
T <- rnorm(n = n, mean = log(theta0), sd = sd)
R <- rnorm(n = n, mean = 0, sd = sd)
# between subjects
between <- rnorm(n = n, mean = 0, sd = sd.b)
T <- T + between
R <- R + between
data <- data.frame(subject = rep(1:n, each = 2), period = 1:2L,
sequence = c(rep("RT", n), rep("TR", n)),
treatment = c(rep(c("R", "T"), n / 2),
rep(c("T", "R"), n / 2)), logPK = NA_real_)
subj.T <- subj.R <- 0L # subject counters
for (j in 1:nrow(data)) { # clumsy but transparent
if (data$treatment[j] == "T") {
subj.T <- subj.T + 1L
if (data$period[j] == 1L) {
data$logPK[j] <- T[subj.T] + per.effect[1]
} else {
data$logPK[j] <- T[subj.T] + per.effect[2] + carryover[1]
}
} else {
subj.R <- subj.R + 1L
if (data$period[j] == 1L) {
data$logPK[j] <- R[subj.R] + per.effect[1]
} else {
data$logPK[j] <- R[subj.T] + per.effect[2] + carryover[2]
}
}
}
facs <- c("subject", "period", "sequence", "treatment")
data[facs] <- lapply(data[facs], factor)
return(data)
} # eof sim.study()
fixed.mixed <- function(data) {
# analyse data by fixed and three mixed effects models
od <- options("digits")
oc <- options("contrasts")
options(digits = 12)
options(contrasts = c("contr.treatment", "contr.poly"))
on.exit(od)
on.exit(oc)
res <- data.frame(package = c("stats", "nlme", rep("lmerTest", 2)),
fun = c("lm", "lme", rep("lmer", 2)),
subjects = c("fixed", rep("random", 3)), df = NA_real_,
method = c(rep("Residual", 2), "Satterthwaite", "Kenward-Roger"),
PE = NA_real_, lower = NA_real_, upper = NA_real_)
model <- c("lm", "lme", "lmer.Satt", "lmer.kr") # for switching
for (j in 1:4L) {
switch (j,
lm = {
# all effects fixed: evaluation by lm() of package stats / ANOVA
# bogus nested model acc. to the GL
m <- lm(logPK ~ sequence + subject %in% sequence +
period + treatment, na.action = na.omit, data = data)
PE <- 100 * exp(coef(m)[["treatmentT"]])
ci <- 100 * exp(confint(m, "treatmentT", level = 1 - 2 * alpha))
df <- anova(m)["Residuals", "Df"]
res[j, c(4, 6:8)] <- c(df, PE, ci)
},
lme = {
# subjects random: evaluation by lm() of package nlme
m <- lme(logPK ~ sequence + period + treatment,
random = ~ 1 | subject, na.action = na.omit, data = data)
s <- summary(m)
PE <- 100 * exp(s$tTable["treatmentT", "Value"])
ci <- 100 * exp(intervals(m, which = "fixed",
level= 1 - 2 * alpha)[[1]]["treatmentT", c(1, 3)])
df <- s$tTable["treatmentT", "DF"]
res[j, c(4, 6:8)] <- c(df, PE, ci)
},
lmer.Satt = {
# subjects random: evaluation by lmer() of package lme4
# Satterthwaite df by package lmerTest
m <- lmer(logPK ~ sequence + period + treatment +
(1 | subject), na.action = na.omit, data = data)
s <- summary(m, ddf = "Satterthwaite")
PE <- s$coefficients["treatmentT", "Estimate"]
ci <- 100 * exp(PE + c(-1, +1) *
qt(1 - alpha, s$coef["treatmentT", "df"]) *
s$coef["treatmentT", "Std. Error"])
PE <- 100 * exp(PE)
df <- s$coefficients["treatmentT", "df"]
res[j, c(4, 6:8)] <- c(df, PE, ci)
},
lmer.kr = {
# subjects random: evaluation by lmer() of package lme4
# Kenward-Roger df by package lmerTest
m <- lmer(logPK ~ sequence + period + treatment +
(1 | subject),
na.action = na.omit, data = data)
s <- summary(m, ddf = "Kenward-Roger")
PE <- s$coefficients["treatmentT", "Estimate"]
ci <- 100 * exp(PE + c(-1, +1) *
qt(1 - alpha, s$coef["treatmentT", "df"]) *
s$coef["treatmentT", "Std. Error"])
PE <- 100 * exp(PE)
df <- s$coefficients["treatmentT", "df"]
res[j, c(4, 6:8)] <- c(df, PE, ci)
}
)
}
return(result = res)
} # eof fixed.mixed()
####################
# simulation setup #
####################
CV <- 0.225
theta0 <- 0.95
theta1 <- 0.80
theta2 <- 1.25
target <- 0.80
alpha <- 0.05
per.effect <- 0
carryover <- c(0, 0)
dropouts <- 1
missings <- 1
n <- sampleN.TOST(alpha = alpha, CV = CV, theta0 = theta0, theta1 = theta1,
theta2 = theta2, targetpower = target,
print = FALSE)[["Sample size"]]
per.effect <- c(0, per.effect)
# simulate balanced (nRT = nTR) and complete (nR = nT) data
data0 <- sim.study(alpha = alpha, CV = CV, theta0 = theta0, theta1 = theta1,
targetpower = target, n = n, per.effect = per.effect,
carryover = carryover)
tmp0 <- fixed.mixed(data0)
# remove droput(s): gives imbalanced sequences
remove <- tail(unique(data0$subject), dropouts)
if (length(remove) >= 1) {
data1 <- data0[!data0$subject == remove, ]
} else {
data1 <- data0
}
tmp1 <- fixed.mixed(data1)
# set logPK in 2nd period of last subject(s) to NA: gives incomplete data
data2 <- data1
remove <- tail(unique(data2$subject), missings)
if (length(remove) >= 1) {
data2$logPK[data2$subject == remove & data2$period == 2] <- NA
}
tmp2 <- fixed.mixed(data2)
# aggregate results
res <- data.frame(data = c(rep("balanced", 4),
rep("imbalanced", 4),
rep("incomplete", 4)),
package = "", FUN = "", subject = "", df = NA_real_,
method = "", PE = NA_real_, lower = NA_real_,
upper = NA_real_, width = NA_real_)
res[1:4, 2:9] <- tmp0 # balanced, complete
res[5:8, 2:9] <- tmp1 # imbalanced, complete
res[9:12, 2:9] <- tmp2 # imbalanced, incomplete
res[, 5] <- round(res[, 5], 1)
res[, 10] <- res[, 9] - res[, 8]
res[, 7:10] <- round(res[, 7:10], 5)
names(res)[3] <- "function"
t1 <- paste("CV :", sprintf("%.4g", CV),
"\nT/R-ratio :", sprintf("%.4g", theta0),
"\nlower limit :", sprintf("%.4f", theta1),
"\nupper limit :", sprintf("%.4f", theta2),
"\npower : at least", sprintf("%.4g", target),
"\nalpha :", sprintf("%.4f", alpha),
"\nn :", n, "(complete data, balanced sequences)")
if (dropouts > 0)
t1 <- paste(t1, "\ndropout(s) :", sprintf("%2i", dropouts), "(both periods)")
if (missings > 0)
t1 <- paste(t1, "\nmissing(s) :", sprintf("%2i", missings), "(second period)")
if (!per.effect[2] == 0)
t1 <- paste(t1, "\nperiod effect:", per.effect[2])
if (!carryover[1] == 0 & !carryover[2] == 0)
t1 <- paste(t1, "\ncarryover :", paste(carryover, collapse = ", "))
t1 <- paste(t1, "\n\n")
t2 <- paste("\nbalanced sequences (nRT = nRT)",
"\nimbalanced sequences (one subject missing)",
"\nincomplete data (as above and the second period",
"of the last subject missing)",
"\ndf : degrees of freedom =",
"\n n – 2 for Residual,",
"approximated by Satterthwaite and Kenward-Roger",
"\nmethod: df method applied",
"\nPE : Point Estimate",
"\nlower :",
sprintf("lower limit of the %.4g%% CI", 100*(1-2*alpha)),
"\nupper :",
sprintf("upper limit of the %.4g%% CI\n", 100*(1-2*alpha)))
cat(t1); print(res, row.names = FALSE, right = FALSE); cat(t2)
Dif-tor heh smusma 🖖🏼 Довге життя Україна!
Helmut Schütz
The quality of responses received is directly proportional to the quality of the question asked. 🚮
Science Quotes
Complete thread:
- ICH M13A: Changes to Step 2
Helmut 2024-07-31 13:16
- ICH M13A: Changes to Step 2 Helmut 2024-07-31 14:19
- AUCres mittyri 2024-08-05 21:08
- ‘Percentage covered’ Helmut 2024-08-05 22:13
- ICH M13A: Changes to Step 2 BEQool 2024-09-09 05:48
- fixed or mixed effects modelHelmut 2024-09-09 07:04
- fixed or mixed effects model BEQool 2024-09-09 10:18
- fixed or mixed effects model Helmut 2024-09-09 11:13
- fixed or mixed effects model BEQool 2024-09-10 07:35
- fixed or mixed effects model Helmut 2024-09-09 11:13
- fixed or mixed effects model Helmut 2024-09-10 08:12
- fixed or mixed effects model BEQool 2024-09-09 10:18
- fixed or mixed effects modelHelmut 2024-09-09 07:04
- AUCres mittyri 2024-08-05 21:08
- ICH M13A: Changes to Step 2 Helmut 2024-08-05 12:53
- period within group and formulation mittyri 2024-08-05 20:42
- period within group and formulation Helmut 2024-08-05 22:29
- ICH M13A: Step 4 → 5 Helmut 2024-08-08 11:57
- Formal ICH Procedure Helmut 2024-08-09 09:45
- ICH M13A: Changes to Step 2 Helmut 2024-09-06 08:04
- ICH M13A: Changes to Step 2 Helmut 2024-07-31 14:19