Dealing with heteroskedasticity in negative binomial GLM

Question

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:

EDIT

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 :)

Answer 1

This answer (Negative Binomial Regression and Heteroskedasticity) on the same forum explains very nicely that models such as yours are predicated on a certain type of relationship between the (conditional) mean and variance of the response variable. The author of the answer advises the following:

"So in place of investigating violation of homoskedasticity, you should investigate violation of the assumed mean-variance relationship."

In your context, one way to achieve this would be to use the gamlss package in R, which allows for distributions other than negative binomial for a count response, thereby allowing you to explore different types of mean-variance relationships and determine which one is most suitable for your data.

This package includes distributions such as Poisson (PO), Poisson Inverse Gaussian (PIG), Negative Binomial Type I (NBI), Negative Binomial Type II (NBII), Sichel (SICHEL) and Delaporte (DEL). See Table 2 of https://docs.ufpr.br/~taconeli/CE06218/Artigo.pdf.

For example, fit your models like this:

library("gamlss")
library("gamlss.dist")

mPO <- gamlss(harvest ~ YAratio + scale(inventory) + 
                        offset(log(hours)) + 
                        random(year),
              data = daily_harvest, 
              family = PO, 
              trace = FALSE)

mPIG <- gamlss(harvest ~ YAratio + scale(inventory) + 
                         offset(log(hours)) +
                         random(year),
              data = daily_harvest, 
              family = PIG, 
              trace = FALSE)

etc.

and then compare them on the basis of their (generalized) AIC:

GAIC(mPO, mPIG, ...)

The model with the smallest AIC value is best supported by your data.

See also slide 17 of the presentation on Flexible Regression and Smoothing: Discrete Distributions by Mikis Stasinopoulos et al.(https://gamlss.com/wp-content/uploads/2019/10/discrete-distributions.pdf) for a nice overview of the expression of the (conditional) mean and variance of the count response variable corresponding to different distributions supported by gamlss:

This article on Analyzing Highly Dispersed Crash Data Using the Sichel Generalized Additive Models for Location, Scale and Shape might come in handy as well:

https://ceprofs.civil.tamu.edu/dlord/Papers/Zou_et_al_Sichel_GAMLSS.pdf

Addendum 1

For your question on year having both a fixed linear effect and a random effect on the variance parameter, see Would including "year" as a categorical random effect remove a long-term trend in a mixed effects model? and especially Thierry Onkelinx's post on Same variable in fixed and random effects at https://rpubs.com/INBOstats/both_fixed_random.

While Thierry Onkelinx's post refers to a mean parameter, I would suspect even in your case the random effect of year included in the variance parameter would pick up nonlinearities in temporal trend.

But if nonlinearities are really present, it seems to me a cleaner way would be to allow for that possibility in the fixed effects portion of your model rather than use a statistical artifice such as the inclusion of a random year effect - especially since you mentioned you prefer simpler interpretations.

Capturing nonlinearity of the year effect depends on how many years you have in your study and/or what the nonlinearity looks like.

If you only have a few years, you can just convert your year variable from numeric to factor and include it in the model as a factor - that would be sufficient to capture nonlinearities that may be present in the fixed year effect.

If you have lots of years, you can allow your model to include a smooth, possibly nonlinear effect of year. See https://www.gamlss.com/additive-terms/ for various ways to include such effects in a gamlss model.

Addendum 2

One thing I like about the gamlss framework is that it comes with nice diagnostic plots for its models. If you want to keep things simple and go back to using your original model for your data, why not fit it to your data using the gamlss framework and then do the following:

1. Compare its GAIC value against that of the "best" (but more complex) model you identified via the gamlss framework - how much of a difference do you see in these values (e.g., more than 3 points?)? Which one is smaller?

2. Compare its diagnostic plots produced by the plot() function against that of the "best" (but more complex) model? Are there diagnostic plots where your simpler model seems to do a poorer job? How much poorer? Is that something you can live with (i.e., not much poorer) or is it something of concern?

3. Compare its (overall) worm plot produced by the wp() function against that of the "best" (but more complex) model? Does your simpler model seems to do a poorer job? How much poorer? Is that something you can live with (i.e., not much poorer) or is it something of concern?

4. Compare its worm plots produced by the wp() function separately for each predictor variable in your fixed effects portion of the model against that of the "best" (but more complex) model? Does your simpler model seems to do a poorer job for a specific predictor variable? How much poorer? Is that something you can live with (i.e., not much poorer) or is it something of concern? The wp() function has an xvar = option you can use to specify the predictor variable of interest. According to the help function for this function: van Buuren and Fredriks M. (2001) used these residuals to identify regions (intervals) of the explanatory variable within which the model does not fit adequately the data (called "model violation").

Only after doing something like this you would feel fully comfortable about choosing to stick with your simpler model (even if you decide to fit that simpler model in a framework other than gamlss). Simpler is not always better - but if it isn't, at least you can quantify precisely where simpler falls short and qualify the findings of your statistical analyses appropriately.

Personally, I find gamlss valuable exactly for this reason - I use it a lot behind the scenes to build a better intuition for what is going on when I model my data. If you don't think it's valuable, you don't have to use it. At least you saw it once in action and decided it is not for you.