[R-sig-ME] Poisson mixed models
Ken Beath
ken at kjbeath.com.au
Wed Oct 22 01:54:53 CEST 2008
On 21/10/2008, at 11:37 PM, Douglas Bates wrote:
> On Tue, Oct 21, 2008 at 5:24 AM, Renwick, A. R.
> <a.renwick at abdn.ac.uk> wrote:
>>
>>> Dear All
>>> There has been a lot of talk recently on this forum regarding (over)
>>> dispersion and quasi models. I am running a GLMM with a poisson
>>> family for some tick burden data I have and I wanted to check if I
>>> had
>>> overdispersion in my model (and thus a poisson family would be
>>> inappropriate). The only method I have found to do this is to run
>>> the
>>> model with a quasipoisson family and then ask for the scale
>>> parameter
>>> using:
>>>
>>> lme4:::sigma(model)
>>>
>>> However, when I do this my model appears severely UNDER dispersed:
>>> sigmaML
>>> 3.779694e-06
>>>
>>> Without the random effect in the model (i.e a GLM) the scale
>>> parameter
>>> is 1.07 - almost perfect for a poisson family. Is the method I am
>>> trying not appropriate to determine the dispersion in the mixed
>>> model?
>>> Does anyone know a better method?
>>>
>>> Many thanks,
>>> Anna
>
> That seems to be an unusually low value for the dispersion.
>
> I would have to check the code to see exactly what the sigma function
> returns in the case of the quasipoisson family. It is quite possible
> that it is an inappropriate value.
>
> I think it is more straightforward to look at the penalized, weighted
> residual sum of squares, model at deviance["pwrss"], divided by the
> number of observations, model at dims["n"]. You can do that for a model
> fit with the Poisson family.
>
> It also looks as if you will want to reduce the number of
> fixed-effects terms in the model.
>
I posted this before but it seems to have been missed.
A quick simulation shows that something is wrong. What I did was to
generate random effects Poisson data and fit both as Poisson and quasi
Poisson, and should give similar results except with slightly larger
SE for quasi. Note the change in the subject random effect std dev.
This seems to be close to the square of the mean of the Poisson data.
Standard errors for fixed effects are totally wrong.
Ken
> nsubj <- 100
> npersubj <- 20
>
> subject <- factor(rep(1:nsubj,each=npersubj))
>
> means <- exp(rep(10+rnorm(nsubj),each=npersubj))
>
> y <- rpois(nsubj*npersubj,means)
>
> simdata <- data.frame(y,subject)
>
> lmer1 <- lmer(y~(1|subject),data=simdata,family=poisson)
> summary(lmer1)
Generalized linear mixed model fit by the Laplace approximation
Formula: y ~ (1 | subject)
Data: simdata
AIC BIC logLik deviance
3305 3316 -1650 3301
Random effects:
Groups Name Variance Std.Dev.
subject (Intercept) 0.9891 0.99453
Number of obs: 2000, groups: subject, 100
Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 10.12714 0.09945 101.8 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
>
> lmer2 <- lmer(y~(1|subject),data=simdata,family=quasipoisson)
> summary(lmer2)
Generalized linear mixed model fit by the Laplace approximation
Formula: y ~ (1 | subject)
Data: simdata
AIC BIC logLik deviance
3307 3324 -1650 3301
Random effects:
Groups Name Variance Std.Dev.
subject (Intercept) 15137 123.03
Residual 15304 123.71
Number of obs: 2000, groups: subject, 100
Fixed effects:
Estimate Std. Error t value
(Intercept) 10.13 12.30 0.8231
More information about the R-sig-mixed-models
mailing list