Bob,

I messed around with this a while ago.  I can't remember the details,  
but with simulations I convinced myself that using a quasi-binomial  
distribution for binary data with lmer could give very wrong results  
(underestimated standard errors, etc.).  I didn't find this too  
surprising, given (as you said) it doesn't really make sense to do  
this.  A simple simulation might be the easiest way to answer your  
problem.

Cheers, Tony


On Sep 28, 2010, at 7:46 AM, Bob Farmer wrote:

> Hi again.
> I've been working with a Bernoulli dependent variable and modeling its
> values as a mixed model using glmer().  My understanding is that it's
> not possible to model overdispersion in a Bernoulli dataset (for
> instance, see Gelman and Hill 2007, p302, or
> http://finzi.psych.upenn.edu/R/Rhelp02a/archive/91242.html ), however
> I can find some suggestions that while it's technically not possible
> to do this, estimating some analog of this parameter may nonetheless
> improve model fit (see Skrondal and Rabe-Hesketh 2007, and a paragraph
> in Venables and Ripley 2002, p297-8 (table 10.4, the "bacteria"
> data)), even if that parameter doesn't actually represent
> "overdispersion" per se.  I don't really understand how this is
> possible for Bernoulli datasets -- I think the mechanics of this
> process are well beyond me -- so I'd have to take it on faith that
> it's working properly, and that makes me uncomfortable.
>
> What I find interesting is that in the following self-contained
> example, glmer() happily models with a quasibinomial family and gives
> tighter SEs about the fixed effects compared to binomial(), but
> glmmPQL (and incidentally, a glm() without the arbitrary grouping
> factor) doesn't (or maybe it does? -- it has broader SEs, but also a
> residual SD term).  In general, my suspicion is that something's amiss
> -- I'm somehow "abusing" the glmer -- and I'm better just sticking
> with the binomial family for Bernoulli datasets.  Any thoughts?
>
> Thanks very much.
> --Bob Farmer
> Dalhousie University, Halifax NS
>
> rm(list=ls())
> table(bacteria$y) #Bernoulli dependent
> m1<-glm(y ~ trt * week, data = bacteria, family=binomial)
> #above example from ?bacteria
> m2<-update(m1, family=quasibinomial)
> summary(m1)
> summary(m2)  #standard errors are similar enough
>
> #use ID as a grouping factor
> library(lme4)
> m1.mix<-glmer(y ~ trt * week + (1 | ID), data =  
> bacteria,family=binomial)
> m2.mix<-update(m1.mix, family=quasibinomial)
> summary(m1.mix)
> summary(m2.mix)  #vastly different standard errors, plus residual
> sigma (the overdisp?)
>
> library(MASS)
> m1.mix.PQL<-glmmPQL(y ~ trt * week, random = ~ 1 | ID, data =  
> bacteria,
>   family=binomial)
> m2.mix.PQL<-update(m1.mix.PQL, family=quasibinomial)
> summary(m1.mix.PQL) #also with residual sigma, but m1-style SEs
> summary(m2.mix.PQL)
>
> Note:  As far as I can tell, in these non-BUGS-based models, this is
> multiplicative overdispersion (which estimates an overdispersion
> parameter based on how the model deviance differs from its expected
> value (which is distributed as a chi-square when overdispersion == 1)
> and then scales down the standard deviation of fixed-effects
> parameters accordingly.  My question does not concern the alternative
> of adding a sigma[i] parameter to the main model formula directly.
>
> _______________________________________________
> R-sig-mixed-models@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models

Anthony Ragnar Ives
Department of Zoology
UW-Madison
(608) 262-1519


	[[alternative HTML version deleted]]


