[R] GLMM (lme4) vs. glmmPQL output (summary with lme4 revised)
Dieter Menne
dieter.menne at menne-biomed.de
Fri Jan 30 10:20:13 CET 2004
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)
}
More information about the R-help
mailing list