[R] GLMM (lme4) vs. glmmPQL output (summary with lme4 revised)

[Dieter Menne]

This is a summary and extension of the thread
"GLMM (lme4) vs. glmmPQL output"

http://maths.newcastle.edu.au/~rking/R/help/04/01/0180.html

In the new revision (#Version: 0.4-7) of lme4 the standard
errors are close to those of the 4 other methods. Thanks to Douglas Bates,
Saikat DebRoy for the revision, and to Göran Broström who run a
simulation.

In response to my first posting, Prof. B. Ripley wrote:
___
Although it has not been stated nor credited, this is very close to an
example in MASS4 (there seems a difference in coding).  Both the dataset
and much of the alternative analyses are from the work of my student James
McBroom (and other students have contributed).
____

Well, I thought having repeated "MASS" four times in the header of
my submitted test programm was enough credit for the r-help audience,
but he certainly is right: the base example glmmPQL is from MASS,
http://www.stats.ox.ac.uk/pub/MASS4/



#Package: lme4
#Version: 0.4-7
#Date: 2004/01/26  !!!!!! Revised

--- GLMM/lme4/Pinheiro/Bates
                 Estimate Std. Error  DF z value  Pr(>|z|)
(Intercept)       3.41202    0.65874 169  5.1796 2.223e-07
trtdrug          -1.24736    0.81824  47 -1.5244 0.1273995
trtdrug+         -0.75433    0.81993  47 -0.9200 0.3575758
I(week > 2)TRUE  -1.60726    0.45525 169 -3.5305 0.0004148

--- glmmPQL/MASS/Venables&Ripley
                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/glmmML/Broström
                  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/geepack/Jun Yan
                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/repeated/J.K.Lindsey
            Estimate Std. Error z value Pr(>|z|)
(Intercept)   -3.569      0.549   -6.51  7.8e-11 ***
trtdrug        1.367      0.486    2.81  0.00490 **
trtdrug+       0.786      0.498    1.58  0.11408
week2TRUE      1.623      0.459    3.53  0.00041 ***
sd             1.294      0.250    5.17  2.3e-07 ***
---

---------------------------------------------------------------

data(bacteria,package="MASS")
UseMASS<-T # must restart R after changing because of nlme/lme4 clash
if (UseMASS){
  library(MASS) # required for bacteria
  options(digits=3)
  cat("--- glmmPQL/MASS/Venables&Ripley\n")
  print(summary(glmmPQL(y ~ trt + I(week > 2), random = ~ 1 | ID,
                  family = binomial, data = bacteria)))
  cat("--- glmmML/glmmML/Broström\n")
  library(glmmML)
  print(glmmML(y=="y"~trt+I(week>2), data=bacteria,cluster=bacteria$ID))
  cat("--- geese/geepack/Jun Yan\n")
  library(geepack)
  print(summary(geese(y == "y" ~ trt + I(week > 2), family = binomial,
     id = ID, corstr = "exchangeable",data=bacteria)))
  library(repeated)
  cat("--- glmm/repeated/J.K.Lindsey\n")
  # reformat data as required by glmm. reshape should be safer...

  bac<-as.data.frame(xtabs(~ID+trt+I(week>2)+y,data=bacteria))
  names(bac,)[3]<-"week2"
  nr<-nrow(bac)
  bac<-cbind(bac[1:(nr/2),-4],bac$Freq[(nr/2+1):nr])
  names(bac,)[4]<-"Freq.n"
  names(bac,)[5]<-"Freq.y"
  attach(bac) # couldn't find the right sytax without this
  print(summary(glmm(cbind(Freq.n,Freq.y)~trt+week2, data=bac, nest=ID,
          family=binomial)))

} else
{
  cat("--- GLMM/lme4/Pinheiro/Bates\n")
  library(lme4)
  bac.GLMM<-GLMM(y ~ trt + I(week > 2), random = ~ 1 | ID,
                  family = binomial, data = bacteria,method="PQL")
  print(bac.GLMM)
 }