[R] glm(family=binomial) and lmer

[Prof Brian Ripley]

On Tue, 14 Aug 2007, Chris O'Brien wrote:

> 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):

and in the places he got it from (it is not his original work).

These are not the only two ways, and they are not the same analyses as the 
saturated models differ.  The usual way to use weights is

y <- dead/batch
model3 <- glm(y ~ log(dose), binomial, weights=batch)
summary(model3)

and internally glm converts models with a two-column response to this 
form, for it is in this form the binomial fits into the GLM framework.

See the White Book or MASS (even the 1994 edition).


> >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.

This form has ten obeservations of groups with weights 2,98,10,80 ....

> 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
>
> ______________________________________________
> R-help at stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>

-- 
Brian D. Ripley,                  ripley at stats.ox.ac.uk
Professor of Applied Statistics,  http://www.stats.ox.ac.uk/~ripley/
University of Oxford,             Tel:  +44 1865 272861 (self)
1 South Parks Road,                     +44 1865 272866 (PA)
Oxford OX1 3TG, UK                Fax:  +44 1865 272595