Bioequivalence and Bioavailability Forum

To all RUsers,

had a nice conversion with Russell Lenth (maintainer of emmeans and author of the VBA-code within FARTSSIE for the non-central t).
What he suggested was actually what one gets in replicateBE by method.B(..., option = 1) – though coded with library lmerTest in our package – i.e., the EMA’s but with Satterthwaite’s degrees of freedom.

The FDA is happy with R but being close is not close enough.

If you have a lot of spare time and any ideas, great. Be aware that many tried before to no avail. 
Below the code I posted to R-help ready to be treated like a child (“who cares about SAS, R rulez, etc.”).

# 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) ~ sequence + period + 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(...) of base-R
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"))


Gives:

EMA example data set 1
Analysis of log-transformed data
Satterthwaite’s degrees of freedom, 90% CI

  SAS 9.4, Phoenix/WinNonlin 8.1
                   mean         SE         df    p.value
    R      :  7.6704296 0.10396421  74.762420
    T      :  7.8158939 0.09860609  74.926384
    T vs. R:  0.1454643 0.04650124 207.734958 0.00201129
                     PE  lower.CL  upper.CL
    antilog:  1.1565764 1.0710440 1.2489393

  lmer (lme 1.1-21), emmeans 1.4
                   mean         SE         df    p.value
    R      :  7.6700137 0.10129485  83.037216
    T      :  7.8161019 0.10139525  83.354517
    T vs. R:  0.1460882 0.04651301 216.938614 0.04951878
                     PE  lower.CL  upper.CL
    antilog:  1.1572982 1.0717073 1.2497248

Slightly different PE and CI little bit conservative.

---

Edit: Wasn’t aware that there is a mixed models list for R. THX to Thierry Onkelinx for pointing that out and forwarding my question. I tweaked his suggestion (changed the default optimizer from "nlminb" to the old one "optim" since the former failed to converge).

library(nlme)
mod  <- lme(log(PK) ~ sequence + period + treatment,
                      random = ~ treatment | subject,
                      data = data, weights = varIdent(~ treatment),
                      method = "REML", na.action = na.exclude,
                      control = list(opt = "optim"))

Gives:

  lme (nlme 3.1-140), emmeans 1.4
                   mean         SE         df    p.value
    R      :  7.6698013 0.10396673  74.970750
    T      :  7.8164670 0.09876994  74.613143
    T vs. R:  0.1466657 0.04707054 127.728002 0.05334051
                     PE  lower.CL  upper.CL
    antilog:  1.1579668 1.0710890 1.2518914

A small step for a man but a giant leap for mankind? Whereas lmer() assumes identical within-subject variances of T and R, varIdent(~ treatment) in lme() estimates them separately. I could live with something which is more conservative than the FDA’s. Of course, this needs a lot of testing to check whether this is always the case.