[R-sig-ME] Poisson mixed models

[Ken Beath]

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