[R] identical results with PQL and Laplace options in lmer function (package lme4)

[Emmanuel Tillard]

Dear R users

I encounter a problem when i perform a generalized linear mixed model (binary data) with the lmer function (package lme4)
with R 2.1.0 on windows XP and the latest version of package "lme4" (0.96-1) and "matrix" (0.96-2)
both options "PQL" and "Laplace" for the method argument in lmer function gave me the same results (random and fixed effects estimates, standard error and p.values). However, Loglikelihood and deviance are different.

here is an example reproduced with the bacteria data set available in the MASS package:

library(lme4)
library(MASS)

data(bacteria)
bacteria$week2 <- as.factor(ifelse(bacteria$week <=2, 0, 1))
model.PQL <- lmer(y ~ trt + week2 + (1 | ID), family = binomial, data  = bacteria, method ="PQL")
model.Laplace <- lmer(y ~ trt + week2 + (1 | ID), family = binomial, data = bacteria, method ="Laplace")

model.PQL
model.Laplace

> model.PQL
Generalized linear mixed model fit using PQL 
Formula: y ~ trt + week2 + (1 | ID) 
   Data: bacteria 
 Family: binomial(logit link)
     AIC      BIC   logLik deviance
 152.443 138.8685 -80.2215  160.443
Random effects:
     Groups        Name    Variance    Std.Dev. 
         ID (Intercept)      3.2721      1.8089 
# of obs: 220, groups: ID, 50

Estimated scale (compare to 1)  0.7800484 

Fixed effects:
            Estimate Std. Error z value  Pr(>|z|)    
(Intercept)  3.41227    0.65884  5.1792 2.228e-07 ***
trtdrug     -1.24743    0.81841 -1.5242 0.1274555    
trtdrug+    -0.75440    0.82009 -0.9199 0.3576229    
week21      -1.60737    0.45527 -3.5306 0.0004146 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

> model.Laplace
Generalized linear mixed model fit using Laplace 
Formula: y ~ trt + week2 + (1 | ID) 
   Data: bacteria 
 Family: binomial(logit link)
      AIC      BIC    logLik deviance
 304.9488 291.3743 -156.4744 312.9488
Random effects:
     Groups        Name    Variance    Std.Dev. 
         ID (Intercept)      3.2721      1.8089 
# of obs: 220, groups: ID, 50

Estimated scale (compare to 1)  0.7800484 

Fixed effects:
            Estimate Std. Error z value  Pr(>|z|)    
(Intercept)  3.41227    0.65884  5.1792 2.228e-07 ***
trtdrug     -1.24743    0.81841 -1.5242 0.1274555    
trtdrug+    -0.75440    0.82009 -0.9199 0.3576229    
week21      -1.60737    0.45527 -3.5306 0.0004146 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

is anybody else aware of this? or did I forget something important ?

Many thanks for your help.

-- 
Emmanuel Tillard
Veterinaire
CIRAD-EMVT
Unite de recherche 18

UMR868 Elevage des Ruminants en Regions Chaudes (ERRC)
Campus ENSA-INRA
2 place Viala
34060 Montpellier cedex 1

tel:	0499612265 (fixe)
	0633850598 (gsm)
fax:	0467545694
e-mail: tillard at cirad.fr