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.