4
$\begingroup$

I’m using the glmer function from the lme4 package in R to model species richness adjacent to aquaculture sites. I have 6 sites: 2 in production, 2 were in production the last years but not anymore at the time of the sampling (fallow), and 2 that were never under production (references). Photographs along transects away from the aquaculture sites were taken each 20-40 m from 0 to 200 m and reference sites were at 1500 m from aquaculture sites. These transects were repeated 7 times over a period of 2 years to determine if the community changed over time.

I’ve followed the steps described in the excellent book from Zuur et al. (2009) Mixed Effects Models and Extensions in Ecology with R and my best model is:

(Note that predictors Distance, Depth and Beggiatoa.sp. have been standardized to remove an lme4 error message.)

glmm.8 <- glmer(sr ~ Distance+Depth+fSubstrate+Beggiatoa.sp.+
                     Distance:Beggiatoa.sp.+(1|fSite),
                glmerControl(optimizer="bobyqa", optCtrl=list(maxfun=100000)),
                family=poisson, data=datsc)

summary(glmm.8)

Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) 
   ['glmerMod']
Family: poisson  ( log )
Formula: sr ~ Distance + Depth + fSubstrate + Beggiatoa.sp. + 
   Distance:Beggiatoa.sp. +      (1 | fSite)
Data: datsc
Control: glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 1e+05))

   AIC      BIC   logLik deviance df.resid 
2279.7   2328.8  -1129.9   2259.7      992 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-1.5171 -0.6376 -0.2008  0.4326  4.9375 

Random effects:
Groups Name        Variance Std.Dev.
fSite  (Intercept) 0.1831   0.4279  
Number of obs: 1002, groups:  fSite, 6

Fixed effects:
                       Estimate Std. Error z value Pr(>|z|)    
(Intercept)             1.49171    0.50388   2.960 0.003072 ** 
Distance                2.05809    0.59940   3.434 0.000596 ***
Depth                  -0.09093    0.02966  -3.066 0.002171 ** 
fSubstrateCoarse       -0.09929    0.08299  -1.196 0.231514    
fSubstrateFine         -0.62376    0.08606  -7.248 4.24e-13 ***
fSubstrateFloc         -1.75314    0.30211  -5.803 6.51e-09 ***
fSubstrateMedium       -0.35201    0.07625  -4.617 3.90e-06 ***
Beggiatoa.sp.           2.42190    1.09521   2.211 0.027011 *  
Distance:Beggiatoa.sp.  3.30995    1.37755   2.403 0.016271 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) Distnc Depth  fSbstC fSbstrtFn fSbstrtFl fSbstM Bggt..
Distance     0.885                                                       
Depth        0.052  0.027                                                
fSubstrtCrs -0.076 -0.028 -0.024                                         
fSubstratFn -0.177 -0.058 -0.145  0.325                                  
fSubstrtFlc  0.047  0.132 -0.022  0.066  0.155                           
fSubstrtMdm -0.088 -0.029 -0.102  0.314  0.380     0.097                 
Beggiat.sp.  0.927  0.947  0.039 -0.024 -0.088     0.080    -0.027       
Dstnc:Bgg..  0.925  0.950  0.037 -0.024 -0.089     0.119    -0.027  0.996

My question is: How do I validate this model to see if it meets the required assumptions?

I did a series of plots but I'm not sure if they are the appropriate ones and if they are, if they violate the assumptions.

EP <- residuals(glmm.8,type="pearson")
plot(EP~fitted(glmm.8))

enter image description here

qqnorm(EP)
qqline(EP)

enter image description here

plot(datsc$Distance, EP, xlab="Distance", ylab="Pearson Residuals")

enter image description here

plot(datsc$Depth, EP, xlab="Depth", ylab="Pearson Residuals")

enter image description here

plot(datsc$fSubstrate, EP, xlab="Substrate", ylab="Pearson Residuals")

enter image description here

plot(datsc$Beggiatoa.sp., EP, xlab="Beggiatoa.sp.", ylab="Pearson Residuals")

enter image description here

plot(fitted(glmm.8)~predict(glmm.8))

enter image description here

I looked on this and other websites and I couldn't find a "perfect" method to validate Poisson GLMM models. I believe a good answer to my question would be relevant to many people. If needed I can provide a subset of my data but this question can probably be answered without it. Still, let me know if need it.

$\endgroup$
  • 1
    $\begingroup$ Related: Diagnostic plots for count regression. $\endgroup$ – gung - Reinstate Monica May 17 '16 at 18:24
  • $\begingroup$ What do you mean by validate? Did you do model building and are you now trying to assess the extent of overfitting and how well the model performs? $\endgroup$ – Björn May 18 '16 at 5:26
  • $\begingroup$ @gung thanks for the edits and pointing to this excellent post! Based on it, it seems that my data would better fit a negative binomial distribution and that I should used Deviance residuals instead of Pearson's for my diagnostic plots. However, some functions used don't work for mixed models (e.g. #3 Overdispersion, #4 Influential and leverage points, and #8: Diagnostic plot codes for log linear models for count data). I've found the gof function from the AER package to verify overdispersion, but can you recommend other resources especially for mixed models for count data? $\endgroup$ – Mud Warrior May 18 '16 at 13:46
  • $\begingroup$ @Björn by validate I mean what procedure to use to assess whether my model meet the underlying assumptions of non-normal GLMM? For linear models using normal distribution you should verify that your residuals vs fitted values are homogeneous, that residuals are normally distributed, and that they are independent vs each of your explanatory variables. If you can answer yes to those 3 assumptions you know that your model is valid. From my readings, I understand that these assumptions don't hold very well for non-normal distribution, so what should I do to say whether my non-normal GLMM is valid? $\endgroup$ – Mud Warrior May 18 '16 at 14:10

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.