[R] glmmADMB: Generalized Linear Mixed Models using AD Model Builder

[Hans Julius Skaug]

Douglas Bates wrote:

>
>The "Laplace" method in lmer and the default method in glmm.admb,
>which according to the documentation is the Laplace approximation,
>produce essentially the same model fit.  One difference is the
>reported value of the log-likelihood, which we should cross-check, and
>another difference is in the execution time
>

Yes, glmmADMB has sqrt(2*pi) constants missing. Thanks for pointing that out.

Execution time: As pointed out by Roel de Jong, the underlying software AD Model Builder
does not use hand-coded derivatives for the Hessian involved in the Laplace approximation, 
but calculates these by automatic differentiation. There is a cost in terms of execution
speed, but on the other hand it is very quick to develop new models, as you do not
have to worry about derivatives. I hope to exploit this beyond standard GLMMs, and provide
other R packages.

Comparison of glmmADMB with lmer: I find that the two packages do not give the same
result on one of the standard datasets in the literature (Lesaffre et. al., Appl. Statist. (2001) 50, Part3, pp 325-335).
The full set of R commands used to download data and fit the model is given at the end of this email.

> fit_lmer_lapl <- lmer(y~ treat + time  + (1|subject),data=lesaffre,family=binomial,method="Laplace")
Warning message:
optim or nlminb returned message ERROR: ABNORMAL_TERMINATION_IN_LNSRCH 
 in: LMEopt(x = mer, value = cv) 

PART OF OUTPUT:

Fixed effects:
             Estimate Std. Error  z value Pr(>|z|)    
(Intercept) -0.626214   0.264996  -2.3631  0.01812 *  
treat       -0.304660   0.360866  -0.8442  0.39853    
time        -0.346605   0.026666 -12.9979  < 2e-16 ***

The corresponding estimates with glmmADMB is:

> fit_glmmADMB <- glmm.admb(y~ treat + time,random=~1,group="subject",data=lesaffre,family="binomial",link="logit")

PART OF OUTPUT:

Fixed effects:
  Log-likelihood: -359.649 
  Formula: y ~ treat + time 
(Intercept)       treat        time 
   -2.33210    -0.68795    -0.46134 


So, the estimates of the fixed effects differ. lmer() does infact produces a warning, and it appears that
it method="Laplace" and method="PQL" produce the same results.


Best regards,

hans


# Load data
source("http://www.mi.uib.no/~skaug/cash/lesaffre_dat.s")

# Run lmer
library(lme4)
fit_lmer <- lmer(y~ treat + time + (1|subject),data=lesaffre,family=binomial)
fit_lmer_lapl <- lmer(y~ treat + time  + (1|subject),data=lesaffre,family=binomial,method="Laplace")


# Run glmmADMB
library(glmmADMB)
example(glmm.admb)	# Must be run once in each new directory (this feature will be removed in future version of glmmADMB).
fit_glmmADMB <- glmm.admb(y~ treat + time,random=~1,group="subject",data=lesaffre,family="binomial",link="logit")