generalized linear mixed models: large differences when using glmmPQL or lmer with laplace approximation

On Tue, Oct 7, 2008 at 4:07 PM, Ken Beath <ken at kjbeath.com.au> wrote:

There is a large difference between the estimated std of the random  
effect,
usually a sign that the glmmPQL approximation isn't working.

Or that there is a mistake in the calculation of the standard errors
for the random effects, which is more likely in this case.

The actual optimization is with respect to the relative standard
deviation of the random effects (relative to the scale parameter in
the conditional standard deviation of the response).  For the Poisson
family or the binomial family that scale parameter is fixed at 1 (you
could also consider the situation to be that there isn't a scale
parameter in those cases).  For the quasipoisson and quasibinomial
families you maybe estimate a value there or maybe not.  I don't know.
I believe Ben's simulations showed that I was doing the wrong thing
there

Definitely something wrong. I did some simulations of my own using  
Poisson distributed data. The standard error of the fixed effects also  
seems rather large.

 > 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
  3329 3341  -1663     3325
Random effects:
  Groups  Name        Variance Std.Dev.
  subject (Intercept) 0.9102   0.95405
Number of obs: 2000, groups: subject, 100

Fixed effects:
             Estimate Std. Error z value Pr(>|z|)
(Intercept)   9.9734     0.0954   104.5   <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
  3331 3348  -1663     3325
Random effects:
  Groups   Name        Variance Std.Dev.
  subject  (Intercept) 11794    108.60
  Residual             12957    113.83
Number of obs: 2000, groups: subject, 100

Fixed effects:
             Estimate Std. Error t value
(Intercept)    9.973     10.860  0.9184
 >


Ken