6
$\begingroup$

I am testing for the effect of treatment and fish length on the school size of fish. Treatment is a categorical variable with two levels (e.g., treatment A & treatment B). Fish length is continuous. There are five species of fish in this dataset and school size varies among fish species. Fish species is therefore treated as a random intercept. There are 4,622 data used to fit this model.

That is: glmer.nb(School.Size~Treatment+Fish.Length+(1|Fish.Species),data=dat).

The diagnostics of this model do not indicate a normal distribution of residuals. I think this is because the school size of fish can be vastly greater than predicted values, but not vastly smaller. This dataset is also zero truncated (i.e., there cannot be observations of a school size of zero), which I think is the cause for a small curve without residuals in the bottom left corner of the plot. I also think that this data may be peculiar because it's not "true" count data. Distributions of fish schools are almost certainly going to be over-dispersed relative to typical count data, and observers had to estimate (i.e., round) counts when fish schools were large.

My question is, given the non-normality (but not heteroskedasticity) of the residuals, can I move forward with this model? Or, how can I improve this model? For reference, I have tried to fit a GLMM with a gamma distribution on log-transforming school sizes and the residuals are similarly non-normal.

Thank you very much for any suggestions.

Pearson residuals vs. log-transformed fitted values of negative binomial GLMM. Numbers indicate values of school size (the response variable) and text color indicates species (levels of the random effect). Transparency is set to 25% such that overlaid points appear darker. A loess smoother (not shown) fit to these points shows homogenous variance around a mean of zero.

Pearson residuals vs. log-transformed fitted values of negative binomial GLMM. Numbers indicate values of school size (the response variable) and text color indicates species (levels of the random effect). Transparency is set to 25% such that overlaid points appear darker. A loess smoother (not shown) fit to these points shows homogenous variance around a mean of zero.

Improve this question
$\endgroup$
5
  • 4
    $\begingroup$ Why would you expect the residuals to be normal? The model explicitly says that the conditional response values aren't normal (they're negative binomial, it's right in your title); why would a diagnostic that also suggests they're not normal be a problem? $\endgroup$
    – Glen_b
    Nov 15, 2015 at 3:23
  • 1
    $\begingroup$ Hi Glen_b, thanks for your response. I was confused because in Zuur et al. 2009 (Fig. 9.8 on page 237) they show diagnostic plots for a negative binomial model that indicate a normal distribution of residuals. I wonder if this is because these plots do not show Pearson residuals (I think they are deviance residuals?). $\endgroup$
    – Reginald Fairfield
    Nov 15, 2015 at 18:22
  • 2
    $\begingroup$ Deviance residuals will be nearer to normal than Pearson residuals, but they still may be distinctly non-normal without it meaning there's any issue with the model. If the negative binomial mean is large but your residuals look very non-normal you may have an issue. If the negative binomial mean is small, it may not mean much of anything. $\endgroup$
    – Glen_b
    Nov 16, 2015 at 0:11
  • $\begingroup$ Thanks very much for the clarification. I'll move forward with the model. $\endgroup$
    – Reginald Fairfield
    Nov 16, 2015 at 0:41
  • $\begingroup$ You also may find this answer helpful: stats.stackexchange.com/questions/185491/… $\endgroup$
    – Stefan
    Dec 23, 2015 at 16:11

1 Answer 1

Reset to default
5
$\begingroup$

Reginald, you have probably moved on by now, but for people that stumble across this post I would like to note that the DHARMa package (available from CRAN, see here) that I have created solves this problem, i.e. it will allow you to test if the residuals are compatible with the assumptions of nb (or any other distributions for that matter).

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.

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.