Helmut
★★★
avatar
Homepage
Vienna, Austria,
2019-08-19 13:07

Posting: # 20503
Views: 608
 

 FDA’s mixed model for ABE: another round [R for BE/BA]

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.

Cheers,
Helmut Schütz
[image]

The quality of responses received is directly proportional to the quality of the question asked. 🚮
Science Quotes
Activity
 Admin contact
20,138 posts in 4,246 threads, 1,387 registered users;
online 10 (0 registered, 10 guests [including 2 identified bots]).
Forum time (Europe/Vienna): 15:33 CET

The fundamental cause of trouble in the world today is
that the stupid are cocksure
while the intelligent are full of doubt.    Bertrand Russell

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