# How to validate a Poisson GLMM model?

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))


qqnorm(EP)
qqline(EP)


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


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


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


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.

• – gung - Reinstate Monica May 17 '16 at 18:24
• 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? – Björn May 18 '16 at 5:26
• @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? – Mud Warrior May 18 '16 at 13:46
• @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? – Mud Warrior May 18 '16 at 14:10