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

[Ben Bolker]

  [It's helpful to keep the full content of the thread included ...]

The SAS documentation says:

You can also use TYPE=FA0(2) here to request a estimate that is
constrained to be nonnegative definite.

  This suggests to me that the model specified is probably singular
(variances equal to zero or correlation equal to +/- 1).  (What do SAS
etc. say if you ask them to report random effect var-cov matrices?)

  It is possible to hack lmer to get term-specific residual variances,
*if* you know in advance which residual variance is the smallest (since
this is done by adding an additional variance term, which must be positive).

  glmmTMB can do term-specific residual variances as well, but isn't as
happy with singular fits.

  nlme takes a lot of iterations to think it's converged because (like
glmmTMB) it's fitting the variance-covariance matrix on a log-Cholesky
scale, which means that it is basically trying to optimize to an
infinite parameter value -- it stops when things look "flat enough".

My results (still had to hack the sign of the effect ...)

         estimate  std.error       PE
glmmTMB 0.1454574 0.04609280 1.156568
nlme    0.1466579 0.04673533 1.157958
lmer    0.1454644 0.04650146 1.156577
SAS     0.1454643 0.04650124 1.156576

  I agree that nlme does worst, but I'm a little alarmed that someone at
the FDA is quibbling over a 0.2% difference [all.equal(1.158,1.156)] in
a numerical computation that has a SE that's about 33% of its mean ...

====
library(RCurl)
library(lme4)
library(nlme)
library(glmmTMB)
library(emmeans)
library(broom.mixed)
url  <- getURL("https://bebac.at/downloads/ds01.csv")
data <- read.csv(text = url, colClasses=c(rep("factor", 4), "numeric"))
options(contrasts=c("contr.SAS","contr.poly"))
##     ‘contr.SAS’ is a wrapper for contr.treatment that sets the base
##     level to be the last level of the factor.  The coefficients
##     produced when using these contrasts should be equivalent to those
##     produced by many (but not all) SAS procedures.
mod1 <- glmmTMB(log(PK) ~ period + sequence + treatment +
(treatment|subject),
              disp=~treatment-1,
              data = data)

mod2 <- lme(log(PK) ~ period + sequence + treatment ,
            random = ~ treatment | subject, data = data,
            weights = varIdent(~treatment),
            control=lmeControl(maxIter=10000, msMaxIter=10000,
                               msMaxEval=10000))

data$obs <- factor(seq(nrow(data)))
mod3 <- lmer(log(PK) ~ period + sequence + treatment + (treatment|subject) +
                 (dummy(treatment,"R")+0|obs),
             data=data,
             control=lmerControl(check.nobs.vs.nlev="ignore"))

modList <- list(glmmTMB=mod1, nlme=mod2, lmer=mod3)
ff <- function(x) {
    res <- subset(tidy(x, effect="fixed"), term=="treatmentR",
                  select=c(estimate,std.error))
    ## hack sign
    res$estimate <- -res$estimate
    return(as.data.frame(res))
}
res <- do.call(rbind,c(lapply(modList, ff),
                       list(SAS=data.frame(estimate=0.1454643,
                                           std.error=0.04650124))))
res$PE <- exp(res$estimate)
lapply(modList,logLik)


On 2019-08-19 8:41 a.m., Helmut Schütz wrote:
> Dear Thierry,
> 
> Thierry Onkelinx wrote on 2019-08-19 14:13:
>> Please do check all parameters. I recall that SAS and R use a
>> different style of dummy coding factor variable.
> 
> I’m not a SASian. Hence, I can’t tell. In R factors specified as
> characters are handled in lexical order and internally coded as integers
> starting with 1.
> 
>> If that is not the case, please write the full equation for the SAS
>> model. The language of mathematics is the best way to clearly describe
>> a model and thus to compare models.
> 
> Agree.
> Unfortunately the SAS code is all the FDA gives. I think everything is
> clear except the covariance structure given by FA0(2).
> See
> https://support.sas.com/documentation/cdl/en/statug/63033/HTML/default/viewer.htm#statug_mixed_sect019.htm#statug.mixed.mixedrandomtypefa0
> Table 56.13
> Sorry, beyond me.

> Best regards,
> Helmut
>