I'm analyzing a harvest dataset and I'm trying to figure out which parameters influence hunting success. My data is a daily number of hunted birds and I have multiple covariates, effort (number of hours spent hunting every day; continuous variable), inventories of birds in the area (count data ranging from 200-100 000), ratio of young:adult in the fall population which is an important predictor of hunting success, etc. Since my data is analogous to count data and because it is overdispersed, I'm using a mixed-effects negative binomial glm to analyze this in R. Here is what the dataset looks like:
year day date harvest inventory YAratio hours <dbl> <dbl> <date> <dbl> <dbl> <dbl> <dbl> 1 2000 276 2000-10-02 96 23000 26 76.5 2 2000 277 2000-10-03 95 21500 26 139. 3 2000 278 2000-10-04 323 26000 26 143 4 2000 279 2000-10-05 356 16500 26 135. 5 2000 280 2000-10-06 314 19000 26 131. 6 2000 281 2000-10-07 147 30000 26 66.8 7 2000 284 2000-10-10 87 35000 26 80 8 2000 285 2000-10-11 223 27500 26 156. 9 2000 286 2000-10-12 151 17500 26 155 10 2000 287 2000-10-13 86 19000 26 148.
Here is the model I'm trying to fit:
mod<-glmer.nb(data=daily_harvest, formula = harvest ~ offset(log(hours)) + YAratio + scale(inventory) + (1|year))
The inventories are usually around 5 000 - 60 000 birds, but there are a few inventories that are wild (> 100 000 birds). I believe this is creating heteroskedasticity the residuals. Here is the plot for residuals vs. fitted values from my model:
I know this is caused by the inventory variable since this doesn't happen when I don't include it in the model and heteroskedasticity is then pretty ok. Is there any way that I can deal with this? Is it necessarily a problem? (i.e. maybe the fact that I have only 4 points with such high fitted values is simply not enough to show variability at this end of the scale?)
There is still some heteroskedasticity in the 100-300 range as shown by this plot zoomed in on the left part:
After Isabella's answer I read up on GAMs and
gamlss and while I learned a lot of cool stuff, I'm not sure a non-linear function is what I necessarily need (I'm no stats genius so I might still be very wrong). But looking at my data, I don't see why I should suspect a non-linear pattern (I could see a situation where having more birds in an area could mean more vigilant individuals warning each other and that could affect hunting success, resulting in an increase in hunting success with inventory sizes until a certain threshold after which hunting success would decrease for example; but there is no previous evidence of such a process in the literature and there is no evidence of this in my data either). So I went back to the dataset and tried to make sense of the pattern I was seeing in the residuals vs. fitted values plot. I think the fact that the residuals vary less as the fitted values get bigger is normal with the data I have. Here is a scatterplot of the harvest vs. the number of birds present on the reserve.
As the inventories get larger, there are more opportunities for increasing success, but it does not automatically mean that hunters will be successful (they can either get unlucky, or may be they have terrible aim or something). So you get this pattern where you can harvest a lot (or not) when there's a lot of birds in the area, but you usually don't get a very high harvest if there are only a few birds around. This explains the pattern in residuals vs. fitted values very well : lot of residual variation at low fitted values (i.e. inventories can be very high or very low at low harvest values so you can be far from the expected mean based on inventory size), and less variation as you go further along since you are very unlikely to have a high hunting success when there are fewer birds present. So this and my results and the residual plots all make sense (I think?).
The remaining question is: is it a problem for interpretation of the coefficients and resulting standard errors? Can I use the output of a negative binomial GLM with this type of pattern in the residuals or should I be trying to fit this with
gamlss nonetheless? I'm always more prone to try and use the less complex model that will fit the data (something about not using a jackhammer to crack a nut) but in this case maybe I should just go with the more general models?
Thanks for bearing with me all the way here :)