[R-sig-ME] Selecting random effects in lme4: ML vs. QAICc

[Richard Feldman]

Thank you very much for your informative response.

I guess I'm confused about how to interpret overdispersion. Using the 
Zuur et al. Owls data that you also present in your worked example, I 
ran the following models, keeping the fixed and random effects just as 
you present it:

glmer, family = gaussian
glmer, family = poisson
glmer, family = quasipoisson
glmmPQL, family = quasipoisson
glmmADMB, family = poisson ##zeroInflation=FALSE
glmmADMB, family = negative binomial ##zeroInflation=FALSE
glmmADMB, family = poisson ##zeroInflation=TRUE
glmmADMB, family = negative binomial ##zeroInflation=TRUE


For all but the quasipoisson glmer, I calculated dispersion following 
your example:

 >rdev<-sum(residuals(model)^2)
 >mdf<-length(fixef(model))
 >rdf<-nrow(Data)-mdf
 >rdev/rdf

For the quasipoisson glmer, I extracted dispersion as:

  >lme4:::sigma(model)^2

The results are as follows:

#                          model dispersion
#1                glmer.gaussian  35.534989
#2                    glmer.pois   5.630751
#3 glmer.quasi -- sigma(model)^2  29.830221
#4                 glmmPQL.quasi   1.076906
#5                 glmmADMB.pois   7.585654
#6               glmmADMB.nbinom   1.085072
#7            glmmADMB.pois.zero   8.255389
#8          glmmADMB.nbinom.zero   1.516587

Am I right to interpret this as saying 1) the sigma(model)^2 method is 
inaccurate and 2) the glmmPQL.quasi and glmmADMB.nbinom are sufficiently 
correcting for the overdispersion? I had thought I had read you advising 
using the quasipoisson model to calculate the overdispersion parameter 
(assuming it did so correctly) needed to adjust AIC for QAIC. It would 
seem the dispersion parameter should come from the Poisson regression. 
More likely I just misread you. Also, is it a concern that the 
dispersion is higher in the zero-inflated models? Does this mean 
zero-inflation is not an issue, at least when wanting to calculate AIC 
values?

Again, thank you for all the help!

Richard


Ben Bolker wrote:
>   There are a number of issues here, I will comment in-line ...
> 
> Richard Feldman wrote:
>> I am trying to select between two models that differ in their random
>> effects. Running a likelihood test gives me a different result than
>> using information criteria.
> 
>   In general you shouldn't be applying both methods in the same
> analysis, except for pedagogical/self-teaching purposes. The two tests
> answer different questions, and trying both leads to the temptation to
> cherry-pick.
>> My models:
>>
>> model.1 <- glmer(Y ~ V1 + V2 + V1:V2 + (V2|SITE), data=Data, family =
>> quasipoisson, REML = TRUE)
>>
>> model.2 <- glmer(Y ~ V1 + V2 + V1:V2 + (1|SITE), data=Data, family =
>> quasipoisson, REML = TRUE)
> 
>   Specifying REML=TRUE has no effect in glmer models.  (I don't know
> offhand whether you got a warning, arguably you should have.)  Note that
> in ?glmer there is no REML argument specified for glmer -- there is no
> warning because it gets swallowed by the optional "..." argument.
> 
>> I use quasipoisson because they are highly overdispersed:
>>
>> lme4:::sigma(model.1) #5.886659 lme4:::sigma(model.2) #101.6434
>  
>   If you are going to fit a quasi-likelihood model then generally you
> should find sigma for the most complex model (model 1 in your case) and
> then apply the same estimated value of sigma for all models.
> 
>> The results of the likelihood test: (I know that technically one
>> should not use (RE)ML for quasi-distributions. However the results
>> are nearly identical whether I use quasipoisson or poisson as the
>> family)
> 
>   Clarification: generally one should not use likelihood ratio tests (ML
> vs REML is a separate issue) at all on quasi-likelihood fits, although
> Venables and Ripley (MASS, p. 210) do suggest using an F test on the
> scaled difference in deviance (this is implemented for GLMs with the
> 'test="F"' option to anova.glm()).  (In contrast, they say in the same
> section that one *cannot* use AIC in this case, because one doesn't have
> a real likelihood -- presumably they don't agree with the line of
> reasoning leading to QAIC.)
> 
>    I'm also assuming that your sample size is large enough that the LRT applies (which is often a problem with GLMMs 
> if the number of random-effects levels is small)
> 
> 
>> anova(model.2, model.1)
>>
>> #         Df    AIC    BIC   logLik  Chisq Chi Df Pr(>Chisq) #model.2
>> 9 4648.3 4665.1 -2315.13 #model.1 14  346.8  373.0  -159.39 4311.5
>> 5  < 2.2e-16 ***
> 
>   Besides the fact that you shouldn't do this, it's unlikely that
> anova() is correcting for overdispersion.
>> Now, I run the same models with a poisson distribution and then
>> adjust the AIC by the overdispersion and the number of parameters to
>> obtain QAICc. With all these penalties taken into account, model.2
>> has the lowest QAICc. My gut instinct is to go with model.1, however,
>> because the overdispersion of model.2 is so high that perhaps it
>> shouldn't even be a candidate model. On the other hand, perhaps
>> adjusting the AIC really does put the two models on a more level
>> playing field.
> 
>   How did you calculate overdispersion?  If you're going to do this (and
> there are still some questions about whether this works) you should use
> sigma^2, not sigma, as your estimate of overdispersion.  If you are
> (heaven forbid) following the supplementary material of the Bolker et al
> 2009 TREE article, please look in the worked examples section of
> glmm.wikidot.com for an updated and corrected version.
> 
>  
> 
>   
> 
>   
> 

-- 
Richard Feldman, PhD Candidate
Dept. of Biological Sciences, McGill University
W3/5 Stewart Biology Building
1205 Docteur Penfield
Montreal, QC H3A 1B1
514-212-3466
richard.feldman at mail.mcgill.ca