[R] GLMM (lme4) vs. glmmPQL output

[Dieter Menne]

Goran,

from my reply to a message from Douglas Bates; ">" is quoted from a mail by
DG.

> I believe the distinction is explained in the lme4 documentation but,
> in any case, the standard errors and the approximate log-likelihood
> for glmmPQL are from the lme model that is the last step in the
> optimization.  The corresponding quantities from GLMM are from another
> approximation that should be more reliable.

I have compared glmmPQL, glmmML, geese and GLMM, results and code see below.
I am aware that glmmPQL uses another method to handle the problem, and
geese (geepack) has considerable different assumptions, but the
results are very similar. On the other hand, I had expected that glmmML
results, if reasonable at all, should be close to GLMM. Yet they are not,
but rather come close to the other three.

See for example the last factor (I(week>2): significant effects in all
cases except GLMM.

-----------------------------
(glmmPQL)
Fixed effects: y ~ trt + I(week > 2)
                Value Std.Error  DF t-value p-value
(Intercept)      3.41     0.519 169    6.58  0.0000
trtdrug         -1.25     0.644  47   -1.94  0.0588
trtdrug+        -0.75     0.645  47   -1.17  0.2484
I(week > 2)TRUE -1.61     0.358 169   -4.49  0.0000

(glmmML)
                  coef se(coef)     z Pr(>|z|)
(Intercept)      3.579    0.701  5.10  3.3e-07
trtdrug         -1.369    0.694 -1.97  4.8e-02
trtdrug+        -0.789    0.700 -1.13  2.6e-01
I(week > 2)TRUE -1.627    0.482 -3.38  7.3e-04

(geese from geepack)
                estimate san.se  wald        p
(Intercept)        2.844  0.529 28.92 7.56e-08
trtdrug           -1.113  0.586  3.61 5.76e-02
trtdrug+          -0.634  0.544  1.36 2.44e-01
I(week > 2)TRUE   -1.325  0.368 12.96 3.18e-04

(GLMM)
                 Estimate Std. Error  DF z value Pr(>|z|)
(Intercept)       3.41202    3.93293 169  0.8676   0.3856
trtdrug          -1.24736    1.52156  47 -0.8198   0.4123
trtdrug+         -0.75433    1.21963  47 -0.6185   0.5363
I(week > 2)TRUE  -1.60726    2.19660 169 -0.7317   0.4644

---------------
data(bacteria,package="MASS")
UseMASS<-F# must restart R after changing because of nlme/lme4 clash
if (UseMASS){
  library(MASS) # required for bacteria
  options(digits=3)
  print(summary(glmmPQL(y ~ trt + I(week > 2), random = ~ 1 | ID,
                  family = binomial, data = bacteria)))
  library(glmmML)
  print(glmmML(y=="y"~trt+I(week>2), data=bacteria,cluster=bacteria$ID))
  library(geepack)
  summary(geese(y == "y" ~ trt + I(week > 2), family = binomial,
     id = ID, corstr = "exchangeable",data=bacteria))
} else
{
  library(lme4)
  # try a well documented "old" lme to compare lme4/lme3 for normal things
#  data(Orthodont,package=nlme)
#  fm <- lme(distance ~ Sex*I(age-11), data =
Orthodont,random=~I(age-11)|Subject)
#  summary(fm)
  summary(GLMM(y ~ trt + I(week > 2), random = ~ 1 | ID,
                  family = binomial, data = bacteria,method="PQL"))
 }