Hello,

I'm a bit confusing with overdispersion estimation in glmer(). I prepared a reproduzible example to address my doubt.

First, I create a toy data adding more variation that expected by chance from a Poisson model.

# toy data, poisson with overdispersion
da <- expand.grid(trat=factor(LETTERS[1:3]),
                  block=factor(c("I","II","III","IV")))
da$y <- rpois(da$trat, lambda=18)+trunc(rnorm(da$trat,15,10))

Than, I fit a glm poisson model and I, as expected, noticed overdispersion by comparing residual deviance with residual degrees of freedom. I calculate the overdispersion parameter by pearson residuals.

# fixed poisson model
g0 <- glm(y~trat+block, family=poisson, data=da)
summary(g0)

# overdispersion estimated by pearson residuals
sum(residuals(g0, type="pearson")^2)/df.residual(g0)

So, I fit a quasi-poisson model and standard errors were multiplied by the pearson overdipersion estimate above.

# notice that overdispersion estimated by quasi-poisson
# is the pearson estimate
g0q <- glm(y~trat+block, family=quasipoisson, data=da)
summary(g0q)

Now I used glmer fitting block as random.

require(lme4)

# random block poisson model
gm0 <- glmer(y~trat+(1|block), family=poisson, data=da)
summary(gm0)

# random block quasi-poisson
gm0q <- glmer(y~trat+(1|block), family=quasipoisson, data=da)
summary(gm0q)

# pearson overdispersion
sum(residuals(gm0, type="pearson")^2)/(nrow(da)-length(fixef(gm0))-1)

# glmer overdispersion?
summary(gm0q)@sigma

I noticed that the overdispersion estimate (@sigma) is not based on pearson residuals, and that the SE were correct by the value estimated by glmer (@sigma). How is glmer overdispersion parameter estimated? Can I simple use glmer poisson and correct the SE by pearson residuals overdispersion?

Thanks.
Walmes.



      
	[[alternative HTML version deleted]]


