Bioequivalence and Bioavailability Forum

Hi BEQool,

❝ 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 n^th 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.

For SASians and users of Statistica: 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)