[R] linear mixed model required for the U.S. FDA

[Helmut Schütz]

Dear all,

I’m struggling to set up a model required for the FDA (haha, and the 
Chinese agency). The closest I could get given at the end (which matches 
the one preferred by other regulatory agencies worldwide). The FDA is 
happy with R but "close" is not close /enough/.

Don't hit me. I'm well aware of the community's attitudes towards SAS. 
I'm not a SASian myself (software agnostic) but that's not related to 
SAS; one could set up this model in other (commercial...) software as well.

The FDA’s model allows different subject effects for each treatment 
(i.e., a subject-by-treatment interaction), and therefore, has 5 
variance terms:
   1. within subject for T
   2. within subject for R
   3. between subject for T
   4. between subject for R
   5. covariance for between subject Test and Reference
The last three are combined to give the subject × formulation 
interaction variance component.

The code provides a lot of significant digits only for comparison.

# FDA 2001 (APPENDIX E)
# https://www.fda.gov/media/70958/download
# FDA 2011 (p. 8)
# 
https://www.accessdata.fda.gov/drugsatfda_docs/psg/Progesterone_caps_19781_RC02-11.pdf
###############################################
# PROC MIXED;                                 #
# CLASSES SEQ SUBJ PER TRT;                   #
# MODEL Y = SEQ PER TRT/ DDFM = SATTERTH;     #
# RANDOM TRT/TYPE = FA0(2) SUB = SUBJ G;      #
# REPEATED/GRP=TRT SUB = SUBJ;                #
# ESTIMATE 'T vs. R' TRT 1 -1/CL ALPHA = 0.1; #
###############################################
# Example data set (EMA)
# 
https://www.ema.europa.eu/en/documents/other/31-annex-ii-statistical-analysis-bioequivalence-study-example-data-set_en.pdf
library(RCurl)
library(lme4)
library(emmeans)
url  <- getURL("https://bebac.at/downloads/ds01.csv")
data <- read.csv(text = url, colClasses=c(rep("factor", 4), "numeric"))
mod  <- lmer(log(PK) ~ period + sequence + treatment + (1|subject),
                        data = data)
res1 <- test(pairs(emmeans(mod, ~ treatment, mode = "satterth"),
                    reverse = TRUE), delta = log(0.8))
res2 <- confint(emmeans(mod, pairwise ~ treatment, mode = "satterth"),
                 level = 0.9)
# Workaround at the end because of lexical order
#   I tried relevel(data$treatment, ref = "R") /before/ the model
#   However, is not observed by confint(...)
cat(paste0("\nEMA Example data set 1",
            "\nAnalysis of log-transformed data",
            "\nSatterthwaite\u2019s degrees of freedom, 90% CI",
            "\n\n  SAS 9.4, Phoenix/WinNonlin 8.1",
            "\n                   mean         SE       df p.value",
            "\n    R      :  7.6704296 0.10396421  74.762420",
            "\n    T      :  7.8158939 0.09860609  74.926384",
            "\n    T vs. R:  0.1454643 0.04650124 207.734958 0.00201129",
            "\n                     PE  lower.CL  upper.CL",
            "\n    antilog:  1.1565764 1.0710440 1.2489393",
            "\n\n  lmer (lme 1.1-21), emmeans 1.4",
            "\n                   mean         SE       df p.value",
            "\n    R      :  ", sprintf("%.7f %.8f  %3.6f",
                                        res2$emmeans$emmean[1],
                                        res2$emmeans$SE[1],
                                        res2$emmeans$df[1]),
            "\n    T      :  ", sprintf("%.7f %.8f  %3.6f",
                                        res2$emmeans$emmean[2],
                                        res2$emmeans$SE[2],
                                        res2$emmeans$df[2]),
            "\n    T vs. R:  ", sprintf("%.7f %.8f %3.6f %.8f",
                                        res1$estimate, res1$SE, res1$df,
                                        res1$p.value),
            "\n                     PE  lower.CL  upper.CL",
            "\n    antilog:  ", sprintf("%.7f %.7f %.7f",
exp(-res2$contrasts$estimate),
exp(-res2$contrasts$upper.CL),
exp(-res2$contrasts$lower.CL)), "\n"))

Cheers,
Helmut

-- 
Ing. Helmut Schütz
BEBAC – Consultancy Services for
W https://bebac.at/
C https://bebac.at/Contact.htm
F https://forum.bebac.at/