[R-sig-ME] Test and account for overdispersion in GLMM - lme4 & glmmML packages.

[Lucas Marie-Orleach]

Dear list members,

I am stuck in statistical problems related to overdispersion in GLMM  
with a Poisson distribution. I have vainly explored the archives of  
the mailing list, books and forums and now count on the kindness and  
knowledge of the list members.

Following the guidelines of the TREE?s review (Bolker et al., 2008), I  
would like to test fixed effects in a GLMM fit by Laplace  
approximation by using the package lme4. I encounter 2 problems:

1- Since overdispersion (indicated by the scale value) is a crucial  
feature in GLMM, I need to test it in my models. The scale value is  
not available in the summary of lmer models, but I have read that it  
is given by the command:

> lme4:::sigma(model)

As far as I understood, when the family used in the model is poisson;  
the function ?lme4:::sigma(model)? gives 1, by definition, which is  
meaningless in my case.
Therefore, it was advised to used ?quasi-? families, then run the  
command ?lme4:::sigma(model)? which give the scale parameter of the  
model. And keep on using quasi-families in case of overdispersion.
However, since ?quasi-families? are no longer available in lme4, how  
does one know if a lmer model is overddispersed?

I have tried to avoid this problem by using the glmmML package and fit  
the same model:

> model<-glmmML(var1~var2+var3+var4,cluster=pair, family=poisson(link="log"))
> summary(model)
> Call:  glmmML(formula = var1 ~ var2 + var3 + var4, family =  
> poisson(link = >"log"), cluster = pair)                   coef   
> se(coef)        z    Pr(>|z|)
> (Intercept)     -0.6039   0.2851   -2.1184    0.03410
> var1            -0.2602   0.2728   -0.9536    0.34000
> var2            -0.7628   0.2753   -2.7710    0.00559
> var3            -0.3075   0.2524   -1.2182    0.22300
>
> Scale parameter in mixing distribution:  0.6403 gaussian Std. Error:  
>                              0.2277      LR p-value for H_0: sigma =  
> 0:  0.04954 Residual deviance: 156.7 on 181 degrees of freedom       
> AIC: 166.7

In the first hand, I have seen that overdispersion of glmmML may be  
indicated by the commands:

> sum.model<-summary(model)       sum.model$dispersion NULL

Such commands give me NULL, which I guess means that the  
under/over-dispersion are null. Is it right?
In another hand, the summary indicates that Residual deviance is not  
equal to the df, which I thought meant overdispersion. Moreover, is  
the given value as ?scale parameter in mixing distribution? (here  
0.6403) stands for the scale value? Should it be equal to 0 or 1 when  
no overdispersion?


2- My second problem is the following step. In case of overdispersion,  
after reading various sources, the solution seems to include the  
individual-level as an additional random effect.
First I do not fully understand how one can test other effects while  
the smallest level of observation is taken into account by this  
additional random effect.
Second, if I test the overdispersion of the new model, should the  
outcome no longer indicate overdispersion?
Third, practical question, shall I just include the individual-level  
in the model and then perform a summary which indicates the outcome of  
a Wald test, as below?

> indiv <- 1:nrow(subdata)
> model<-lmer(var1~var2+var3+var4+(1|indiv)+(1|pair),family=poisson(link="log"))
> summary(model)
> Generalized linear mixed model fit by the Laplace approximation  
> Formula: var1 ~ var2 + var3 + var4 + (1 | indiv) + (1 | pair)      
> AIC   BIC  logLik  deviance
>   168.7 188.0 -78.34    156.7
> Random effects:
> Groups Name        Variance   Std.Dev.  indiv  (Intercept)  
> 1.8331e-09 4.2814e-05
> pair   (Intercept) 4.0996e-01 6.4028e-01
> Number of obs: 186, groups: indiv, 186; pair, 93
>
> Fixed effects:
>                Estimate   Std. Error z value  Pr(>|z|)   (Intercept)  
>     -0.6039     0.2699    -2.237     0.0253 * var1             
> -0.2602     0.2842    -0.916     0.3599   var2            -0.7629     
>  0.2885    -2.644     0.0082 **
> var3            -0.3075     0.2671    -1.151     0.2496   ---
> Signif. codes:  0 ?***? 0.001 ?**? 0.01 ?*? 0.05 ?.? 0.1 ? ? 1  
> Correlation of Fixed Effects:
>             (Intr)       var1       var2
> var1    -0.545              var2    -0.495     0.065         var3     
> -0.420    -0.060     0.059



Any help would be greatly appreciated,

Best wishes,

Lucas Marie-Orleach



-- 
Lucas Marie-Orleach


Zoological Institute, University of Basel
Vesalgasse 1, 4051 Basel, Switzerland

Tel:      0041 612 670 375
Fax:      0041 612 670 362
E-mail:   l.marie-orleach at unibas.ch
Web:      http://evolution.unibas.ch/scharer/