[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