Residuals of GAM models not improving with poisson or ziplss, but better with negative binomial (help with high values) I am running GAM models on species counts with lots of zeros and high values or high counts. Residuals under poisson family have a s-like curve on qq line with models not predicting lower and higher values well (see c. bottom figure). Residuals under negative binomial are better and appear more normal for the lower and middle values, while extreme values are off the qq line (see see a. top figure. I thought ziplss() might help normalize the residuals on the qq line, but low and high values are mapping off the qq line (s-like curve and wiggly, see b. middle figure) similar to the residuals for poisson. There are a few predictors that are factors in the dataset assigned as bs="re". The other predictors are continuous and assigned a bs="cr". Species counts are not log transformed, kept as counts. What could be the culprit? How to address this?




I would use a rootogram to look at how well you are modelling the 0s before moving to the ziplss() family from the nb() one. Also, turn on reference band for the plot (assuming via qq.gam()?) with rep = 50 say. That will help confirm whether deviations from the line are within expectations from draws from the model or not.
Rootograms are in package {countreg} on RForge (not CRAN), or a less sophisticated version is available in my {gratia} package.
If you are modelling the 0s pretty well with the NB model, then it might be that you need more flexibility in the extra dispersion part of the NB. This parameter can't be modelled via a linear predictor in {mgcv} (it is modelled as a constant over all observations), though other software might allow it - {brms} will but you will be fitting using MCMC (HMC) so it will likely be much slower that with {mgcv}.
You don't want to transform the counts - leave them as counts! The whole point of GLMs and GLM-like models with link functions is that you can model the mean (expectation) and possibly other parameters of the conditional distribution of the response on a scale where things are approximately linear (like the log scale) without actually transforming the data; we're just modelling the mean (which we didn't observe) on the log scale in the case of the Poisson or NB).
You could also be missing one or more important covariates and if you haven't got those data then there may be nothing that you can do and you should proceed but with a little bit of caution when you start to interpret the model as something may be missing.
