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