[R-sig-ME] Over-dispersed models

Jarrod Hadfield j.hadfield at ed.ac.uk
Wed Aug 5 14:10:33 CEST 2009


Hi,

I am having trouble obtaining the variance component estimates I  
expect when using quasi GLMM models with lmer.   My first thought is  
that I am misinterpreting how lmer handles over-dispersion. For  
example, I thought a quasi Poisson model with log-link would  have the  
form:

mu = exp(Xb+Zu+e)

where mu is the Poisson rate parameter and e is a vector of residuals,  
the variance of which is estimated (the "Residual" term in the model  
summary).  However, if I use simulated data (see below) the estimates  
do not correspond. I f I simulate data  with Xb=0, and var(u)=var(e)=1  
(with a well replicated balanced design), the mean estimate is  
var(u)=3.4 and var(e)= 2.62, instead of 1.

If I fit a quasibinomial model with the linear predictor of the same  
form

Pr = inverse.logit(Xb+Zu+e)

where Pr is the probability of success, then I get var(u)=0.038 and  
var(e)= 0.042,  instead of 1.

I guess the most likely explanation is that lmer is using some other  
model for over-dispersion (multiplicative?), but if someone could  
verify that it would be great. The other odd thing is that the  
likelihood does not seem to change between the standard and quasi  
models on the same data.

Any help would be appreciated,

Jarrod

I am using version 0.999375-31 on a MacBook Pro and Linux Fedora Core8.

library(gtools)

quasipoisson=FALSE
nsim<-100
VGhat<-1:nsim
VRhat<-1:nsim
fac<-gl(250,4)

for(i in 1:nsim){
	
   fac<-gl(250,4)
   l<-rnorm(250, 0, sqrt(1))[fac]+rnorm(1000, 0, sqrt(1))

   if(quasipoisson){
     y<-rpois(1000, exp(l))
     m1<-glmer(y~1+(1|fac), family="quasipoisson")
   }else{
     y2succ<-rbinom(1000, 10, inv.logit(l))
     y2fail<-10-y2succ
     m1<-glmer(cbind(y2succ, y2fail)~1+(1|fac), family="quasibinomial")
   }

   VGhat[i]<-VarCorr(m1)$fac[1]
   VRhat[i]<-attr(VarCorr(m1), "sc")^2
}



-------------- next part --------------
An embedded and charset-unspecified text was scrubbed...
Name: not available
URL: <https://stat.ethz.ch/pipermail/r-sig-mixed-models/attachments/20090805/6a7f1384/attachment.pl>


More information about the R-sig-mixed-models mailing list