Skip to content

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

8 messages · Greg Snow, Douglas Bates, Ben Bolker +2 more

Original
Greg Snow
# 
Martijn,

Absolute truth is the realm of pure mathematics, philosophy, and religion.  R deals with statistics where we are more concerned with usefulness than absolute correctness.

One famous quote among statisticians is: "All models are wrong, some models are useful".

Given that, I would suggest looking at the quality of the predictions from your models.  Get predicted values from the 2 different fits and compare them to each other and the original data.  If one set of predictions is clearly superior, then you have your answer.  If the predictions are not that different, then that can tell you that possibly there are multiple models that are equally useful for this data.  Then you can look at the predictions from models fit with fewer predictors to see how that affects the predictions.

Hope this helps,

--
Gregory (Greg) L. Snow Ph.D.
Statistical Data Center
Intermountain Healthcare
greg.snow at imail.org
801.408.8111
Douglas Bates
# 
Due to some travel and the need to attend to other projects, I haven't
been keeping up as closely with this list as I normally do.  Regarding
the comparison between the PQL and Laplace methods for fitting
generalized linear mixed models, I believe that the estimates produced
by the Laplace method are more reliable than those from the PQL
method.  The objective function optimized by the Laplace method is a
direct approxmation, and generally a very good approximation, to the
log-likelihood for the model being fit.  The PQL method is indirect
(the "QL" part of the name stands for "quasi-likelihood") and, because
it involves alternating conditional optimization, can alternate
back-and-forth between two potential solutions, neither of which is
optimal.  (To be fair, such alternating occurs more frequently in the
analogous method for nonlinear mixed-models, in which I was one of the
co-conspirators, than in the PQL method for GLMMs.)

It may be that the problem you are encountering has more to do with
the use of the quasipoisson family than with the Laplace
approximation.  I am not sure that the derivation of the standard
errors in lmer when using the quasipoisson family is correct, in part
because I don't really understand the quasipoisson and quasibinomial
families.  As far as I know, they don't correspond to probability
distributions so the theory is a bit iffy.

Do you need to use the quasipoisson family or could you use the
poisson family?  Generally the motivation for the quasipoisson familiy
is to accomodate overdispersion.  Often in a generalized linear mixed
model the problem is underdispersion rather than overdispersion.

In one of Ben's replies in this thread he discusses the degrees of
freedom attributed to certain t-statistics.  Regular readers of this
list are aware that degrees of freedom is one of my least favorite
topics.  If one has a reasonably large number of observations and a
reasonably large number of groups then the issue is unimportant.
(Uncertainty in degrees of freedom is important only when the value of
the degrees of freedom is small.  In fact, when I first started
studying statistics we used the standard normal in place of the
t-distribution whenever the degrees of freedom exceeded 30).
Considering that the quasi-Poisson doesn't correspond to a probability
distribution in the first place, (readers should feel free to correct
me if I am wrong about this) I find the issue of the number of degrees
of freedom that should be attributed to a distribution of a quantity
calculated from a non-existent distribution to be somewhat off the
point.

I think the problem is more likely that the standard errors are not
being calculated correctly.  Is that what you concluded from your
simulations, Ben?

On Tue, Oct 7, 2008 at 8:21 AM, Martijn Vandegehuchte
<martijn.vandegehuchte at ugent.be> wrote:
Ken Beath
# 
There is a large difference between the estimated std of the random  
effect, usually a sign that the glmmPQL approximation isn't working.

I would try using adaptive Gauss-Hermite. Set the nAGQ parameter to  
increasing values until the results stabilise.

Ken
On 08/10/2008, at 12:21 AM, Martijn Vandegehuchte wrote:

            

      
      
      
    

        

      
      
    

  
  


  


      

    

    

      
      

        


  

        

      Douglas Bates
    

    

      
      #
    

  


  

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

            

      
      
      
    

        
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.

            

      
      
      
    

        

      
      
    

            

      
      
      
    

        

      
      
    

  
  


  


      

    

    

      
      

        


  

        

      Ben Bolker
    

    

      
      #
    

  


  

      
Douglas Bates wrote:

        
[snip]

            

      
      
      
    

        
In ecological data it's quite common to have dispersion remaining
in a GLMM even after accounting for known grouping factors.

            

      
      
      
    

        
Fair enough, but backing up to wanting p-values associated with
(sensible???) hypothesis tests -- we don't have mcmcsamp for GLMMs,
so our options are (1) no p-values at all, (2) Z tests (i.e. don't
worry about uncertainty in the scale parameter estimate, (3) ??
simulation from the null hypothesis.  See

  http://www.zoo.ufl.edu/bolker/glmm/glmersim.pdf

            

      
      
      
    

        
Yes.

  http://www.zoo.ufl.edu/bolker/glmm/quasitest.pdf

  By the way, I'd be more than happy for any input on the
above-referenced URLs -- if anyone thinks (and can argue
reasonably convincingly) that they're wrong and/or misguided
I'll take them down so as not to mislead the public.
The Sweave files are up in the same place (substitute .Rnw
for .pdf)

  cheers
    Ben

  
  


  


      

    

    

      
      

        


  

        

      Martijn Vandegehuchte
    

    

      
      #
    

  


  

      
Thank you for this explanation. My dataset consists of 120 observations, 
with 6 levels of the random effect (20 observations per level). I think in 
my case the df, like you say, are not the point, because the number is 
rather large.

Considering the quasipoisson family, I do not really know how it works or 
what it does (like I said, I am not a statistician at all). I only learned 
that it is a way to deal with overdispersion, and yes, my data are highly 
overdispersed. If I make the same models with lmer with a poisson family, 
all the effects become highly significant, so the correction seems 
necessary.

Another remark: by now I found out that when I use poisson or quasipoisson, 
the results with glmmPQL are exactly the same, and still quite close to the 
quasipoisson results of lmer (as opposed to the poisson results of lmer). I 
have been searching for information about the quasi families, but I can't 
seem to find any. I'm a bit puzzled by this at the moment. Why is there a 
huge difference between quasipoisson and poisson with lmer and not with 
glmmPQL? It seems that glmmPQL already accounts for the overdispersion, also 
because when I compare the results with the output of SAS proc glimmix, the 
results of glmmPQL are exactly the same only if I put the overdispersion 
correction "random _residual_" in SAS.

Also, the overdispersion in my data is to a large extent due to the amount 
of zeroes. I have been searching for ways to build ZIP models, but in SAS 
with proc nlmixed, it is quite complex, and with glmm.admb() in R, I also 
get strange results and warning messages, and I have problems with obtaining 
the output I need. Now a colleague from my lab just pointed out that a glmm 
with a correction for overdispersion gave very similar results compared with 
ZIP models in a study of his. And besides that, the overdispersion in my 
case is not only due to the zeroes, there are also some large values in my 
dependent variables, so I'm not even sure if a ZIP model is the right way to 
deal with it. But maybe it is a way to avoid the quasifamilies, if they 
should be avoided at all. It's just a thought.

I would appreciate your ideas on these matters,

Thank you for your time and effort,

Martijn.

  
    

      
      
    

  


  


      

    

    
2 days later

    

      
      

        


  

        

      Ken Beath
    

    

      
      #
    

  


  

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

        

            

      
      
      
    

        
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