[R-sig-ME] generalized linear mixed models: large differences when using glmmPQL or lmer with laplace approximation

[Ken Beath]

On 08/10/2008, at 8:18 AM, Douglas Bates wrote:

> 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