[R] glm(family=binomial) and lmer

Chris O'Brien obrienc at email.arizona.edu
Tue Aug 14 17:47:41 CEST 2007


Dear R users,

I've notice that there are two ways to conduct a binomial GLM with binomial 
counts using R.  The first way is outlined by Michael Crawley in his 
"Statistical Computing book" (p 520-521):

 >dose=c(1,3,10,30,100)
 >dead = c(2,10,40,96,98)
 >batch=c(100,90,98,100,100)
 >response = cbind(dead,batch-dead)
 >model1=glm(y~log(dose),binomial)
 >summary(model1)

Which returns (in part):
Coefficients:
             Estimate Std. Error z value Pr(>|z|)
(Intercept)  -4.5318     0.4381  -10.35   <2e-16 ***
log(dose)     1.9644     0.1750   11.22   <2e-16 ***
Null deviance: 408.353  on 4  degrees of freedom
Residual deviance:  10.828  on 3  degrees of freedom
AIC: 32.287

Another way to do the same analysis is to reformulate the data, and use GLM 
with weights:

 >y1=c(rep(0,5),rep(1,5))
 >dose1=rep(dose,2)
 >number = c(batch-dead,dead)
 >data1=as.data.frame(cbind (y1,dose,number))
 >model2=glm(y1~log(dose1),binomial,weights=number,data=data1)
 >summary(model2)

Which returns:

Coefficients:
             Estimate Std. Error z value Pr(>|z|)
(Intercept)  -4.5318     0.4381  -10.35   <2e-16 ***
log(dose1)    1.9644     0.1750   11.22   <2e-16 ***
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 676.48  on 9  degrees of freedom
Residual deviance: 278.95  on 8  degrees of freedom
AIC: 282.95

Number of Fisher Scoring iterations: 6

These two methods are similar in the parameter estimates and standard 
errors, however the deviances, their d.f., and AIC differ.  I take the 
first method to be the correct one.
However, I'm really interested in conducting a GLM binomial mixed model, 
and I am unable to figure out how to use the first method with the lmer 
function from the lme4 library, e.g.

 >model3=lmer(y~log(dose)+time|ID)    # the above example data doesn't have 
the random effect, but my own data set does.

   Does anyone have any suggestions?

thanks,
chris

Thanks,
Chris O'Brien



More information about the R-help mailing list