Poisson mixed models

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