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

From: Hans Julius Skaug <Hans.Skaug_at_mi.uib.no>
Date: Tue 20 Dec 2005 - 05:19:33 EST


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")



R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html Received on Tue Dec 20 05:25:48 2005

This archive was generated by hypermail 2.1.8 : Fri 03 Mar 2006 - 03:41:40 EST