# Diagnostics for generalized linear (mixed) models (specifically residuals)

I am currently struggling with finding the right model for difficult count data (dependent variable). I have tried various different models (mixed effects models are necessary for my kind of data) such as lmer and lme4 (with a log transform) as well as generalized linear mixed effects models with various families such as Gaussian or negative binomial.

However, I am quite unsure on how to correctly diagnose the resulting fits. I found a lot of different opinions on that topic on the Web. I think diagnostics on linear (mixed) regression are quite straight-forward. You can go ahead and analyse the residuals (normality) as well as study heteroscedasticity by plotting fitted values compared to residuals.

However, how do you properly do that for the generalized version? Let us focus on a negative binomial (mixed) regression for now. I have seen quite opposing statements regarding the residuals here:

1. In Checking residuals for normality in generalised linear models it is pointed out in the first answer that the plain residuals are not normally distributed for a GLM; I think this is clear. However, then it is pointed out that Pearson and deviance residuals are also not supposed to be normal. Yet, the second answer states that deviance residuals should be normally distributed (combined with a reference).

2. That deviance residuals should be normally distributed is hinted at in the documentation for ?glm.diag.plots (from R's boot package), though.

3. In this blog post, the author first studied normality of what I assume are Pearson residuals for a NB mixed-effects regression model. As expected (in my honest opinion) the residuals did not show to be normal and the author assumed this model to be a bad fit. However, as stated in the comments, the residuals should be distributed according to a negative binomial distribution. In my opinion, this comes closest to the truth as GLM residuals can have other distributions than the normal one. Is this correct? How to check for things like heteroscedasticity here?

4. The last point (plotting residuals against quantiles of the estimated distribution) is emphasized in Ben & Yohai (2004). Currently, this seems the way to go for me.

In a nutshell: How do you properly study the model fits of generalized linear (mixed) regression models specifically with a focus on residuals?

• Residuals for GLMs aren't in general normal (cf here), but note that there are lots of kinds of residuals for GLMs. Eg, glm.diag.plots says it's for jackknifed deviance residual (I suspect that distinction is important). Also, I gather you have count data; you might want to focus on that fact. Eg, counts are supposed (in some sense) to be heteroscedastic. Diagnostic plots for count regression should be helpful for you (although it doesn't address the mixed effects aspect). – gung Dec 9 '15 at 17:18
• I am familiar with the post you mentioned. However, there is also a statement that suggests that (deviance) residuals should be normal "we see very large residuals and a substantial deviance of the deviance residuals from the normal (all speaking against the Poisson)". – fsociety Dec 9 '15 at 18:06

This answer is not based on my knowledge but rather quotes what Bolker et al. (2009) wrote in an influential paper in the journal Trends in Ecology and Evolution. Since the article is not open access (although searching for it on Google scholar may prove successful, I thought I cite important passages that may be helpful to address parts of the questions. So again, it's not what I came up with myself but I think it represents the best condensed information on GLMMs (inlcuding diagnostics) out there in a very straight forward and easy to understand style of writing. If by any means this answer is not suitable for whatever reason, I will simply delete it. Things that I find useful with respect to questions regarding diagnostics are highlighted in bold.

Page 127:

Researchers faced with nonnormal data often try shortcuts such as transforming data to achieve normality and homogeneity of variance, using nonparametric tests or relying on the robustness of classical ANOVA to nonnormality for balanced designs [15]. They might ignore random effects altogether (thus committing pseudoreplication) or treat them as fixed factors [16]. However, such shortcuts can fail (e.g. count data with many zero values cannot be made normal by transformation). Even when they succeed, they might violate statistical assumptions (even nonparametric tests make assumptions, e.g. of homogeneity of variance across groups) or limit the scope of inference (one cannot extrapolate estimates of fixed effects to new groups). Instead of shoehorning their data into classical statistical frameworks, researchers should use statistical approaches that match their data. Generalized linear mixed models (GLMMs) combine the properties of two statistical frameworks that are widely used in ecology and evolution, linear mixed models (which incorporate random effects) and generalized linear models (which handle nonnormal data by using link functions and exponential family [e.g. normal, Poisson or binomial] distributions). GLMMs are the best tool for analyzing nonnormal data that involve random effects: all one has to do, in principle, is specify a distribution, link function and structure of the random effects.

Page 129, Box 1:

The residuals indicated overdispersion, so we refitted the data with a quasi-Poisson model. Despite the large estimated scale parameter (10.8), exploratory graphs found no evidence of outliers at the level of individuals, genotypes or populations. We used quasi-AIC (QAIC), using one degree of freedom for random effects [49], for randomeffect and then for fixed-effect model selection.

Page 133, Box 4:

Here we outline a general framework for constructing a full (most complex) model, the first step in GLMM analysis. Following this process, one can then evaluate parameters and compare submodels as described in the main text and in Figure 1.

1. Specify fixed (treatments or covariates) and random effects (experimental, spatial or temporal blocks, individuals, etc.). Include only important interactions. Restrict the model a priori to a feasible level of complexity, based on rules of thumb (>5–6 random-effect levels per random effect and >10–20 samples per treatment level or experimental unit) and knowledge of adequate sample sizes gained from previous studies [64,65].

2. Choose an error distribution and link function (e.g. Poisson distribution and log link for count data, binomial distribution and logit link for proportion data).

3. Graphical checking: are variances of data (transformed by the link function) homogeneous across categories? Are responses of transformed data linear with respect to continuous predictors? Are there outlier individuals or groups? Do distributions within groups match the assumed distribution?

4. Fit fixed-effect GLMs both to the full (pooled) data set and within each level of the random factors [28,50]. Estimated parameters should be approximately normally distributed across groups (group-level parameters can have large uncertainties, especially for groups with small sample sizes). Adjust model as necessary (e.g. change link function or add covariates).

5. Fit the full GLMM. Insufficient computer memory o r too slow: reduce model complexity. If estimation succeeds on a subset of the data, try a more efficient estimation algorithm (e.g. PQL if appropriate). Failure to converge (warnings or errors): reduce model complexity or change optimization settings (make sure the resulting answers make sense). Try other estimation algorithms. Zero variance components or singularity (warnings or errors): check that the model is properly defined and identifiable (i.e. all components can theoretically be estimated). Reduce model complexity. Adding information to the model (additional covariates, or new groupings for random effects) can alleviate problems, as will centering continuous covariates by subtracting their mean [50]. If necessary, eliminate random effects from the full model, dropping (i) terms of less intrinsic biological interest, (ii) terms with very small estimated variances and/or large uncertainty, or (iii) interaction terms. (Convergence errors or zero variances could indicate insufficient data.)

6. Recheck assumptions for the final model (as in step 3) and check that parameter estimates and confidence intervals are reasonable (gigantic confidence intervals could indicate fitting problems). The magnitude of the standardized residuals should be independent of the fitted values. Assess overdispersion (the sum of the squared Pearson residuals should be $\chi^2$ distributed [66,67]). If necessary, change distributions or estimate a scale parameter. Check that a full model that includes dropped random effects with small standard deviations gives similar results to the final model. If different models lead to substantially different parameter estimates, consider model averaging.

Residuals plots should be used to assess overdispersion and transformed variances should be homogeneous across categories. Nowhere in the article was mentioned that residuals are supposed to be normally distributed.

I think the reason why there are contrasting statements reflects that GLMMs (page 127-128)...

...are surprisingly challenging to use even for statisticians. Although several software packages can handle GLMMs (Table 1), few ecologists and evolutionary biologists are aware of the range of options or of the possible pitfalls. In reviewing papers in ecology and evolution since 2005 found by Google Scholar, 311 out of 537 GLMM analyses (58%) used these tools inappropriately in some way (see online supplementary material).

And here are a few full worked examples using GLMMs including diagnostics.

I realize that this answer is more like a comment and should be treated as such. But the comment section doesn't allow me to add such a long comment. Also since I believe this paper is of value for this discussion (but unfortunately behind a pay-wall), I thought it would be useful to quote important passages here.

Cited papers:

[15] - G.P. Quinn, M.J. Keough (2002): Experimental Design and Data Analysis for Biologists, Cambridge University Press.

[16] - M.J. Crawley (2002): Statistical Computing: An Introduction to Data Analysis Using S-PLUS, John Wiley & Sons.

[28] - J.C. Pinheiro, D.M. Bates (2000): Mixed-Effects Models in S and S-PLUS, Springer.

[49] - F. Vaida, S. Blanchard (2005): Conditional Akaike information for mixed-effects models. Biometrika, 92, pp. 351–370.

[50] - A. Gelman, J. Hill (2006): Data Analysis Using Regression and Multilevel/Hierarchical Models, Cambridge University Press.

[64] - N.J. Gotelli, A.M. Ellison (2004): A Primer of Ecological Statistics, Sinauer Associates.

[65] - F.J. Harrell (2001): Regression Modeling Strategies, Springer.

[66] - J.K. Lindsey (1997): Applying Generalized Linear Models, Springer.

[67] - W. Venables, B.D. Ripley (2002): Modern Applied Statistics with S, Springer.

• Thanks, that is indeed helpful, I knew about the coding examples of Bolker, but not the actual paper somehow. What I still wonder though is how graphical checking applies to very large-scale data with thousands of groups. The few papers (such as that one) that try give some guidelines on how to properly check your models all only apply to very small-scale data. Then, it is much easier to pick e.g., the groups and visualize something. I really think that a good scientific contribution can be made if someone goes through a more complex example in future. – fsociety Dec 16 '15 at 19:42
• I'm glad this was useful! I think the presented examples are already quite complex (at least to me). I guess the bigger problem is that larger datasets and more complex models may become computationally infeasible as is mentioned in the text: "[...] to find ML estimates, one must integrate likelihoods over all possible values of the random effects. For GLMMs this calculation is at best slow, and at worst (e.g. for large numbers of random effects) computationally infeasible." What I find amazing though, and what should be kept in mind, is that we are using tools that are under active research! – Stefan Dec 16 '15 at 20:02

This is an old question, but I thought it would be useful to add that option 4 suggested by the OP is now available in the DHARMa R package (available from CRAN, see here).

The package makes the visual residual checks suggested by the accepted answer a lot more reliable / easy.

From the package description:

The DHARMa package uses a simulation-based approach to create readily interpretable scaled residuals from fitted generalized linear mixed models. Currently supported are all 'merMod' classes from 'lme4' ('lmerMod', 'glmerMod'), 'glm' (including 'negbin' from 'MASS', but excluding quasi-distributions) and 'lm' model classes. Alternatively, externally created simulations, e.g. posterior predictive simulations from Bayesian software such as 'JAGS', 'STAN', or 'BUGS' can be processed as well. The resulting residuals are standardized to values between 0 and 1 and can be interpreted as intuitively as residuals from a linear regression. The package also provides a number of plot and test functions for typical model mispecification problem, such as over/underdispersion, zero-inflation, and spatial / temporal autocorrelation.

• Very good addition to this thread! – Stefan Feb 15 '17 at 21:05