# How to develop a negative binomial model where the overdispersion parameter varies as a function of one of the independent variables/covariates?

I am trying to develop a negative binomial model where the dependent variable is crash count, and the independent variables are traffic count and roadway length. Currently, with the below code, I get only one value of overdispersion parameter. But I want to get the overdispersion parameter that varies as a function of roadway length. Can you please tell me what should I change in the code? Thank you very much!

Crash_count<-
glm.nb(Total_crashes~LN(traffic_count) +

• Something tells me that this might go beyond what the package provides and you might need to code this by yourself... But I might be wrong. Mar 16, 2022 at 21:04
• Well, I was wondering, if I need to vary the overdispersion parameter relative to the segment length, then is it as simple as dividing the model overdispersion parameter by the segment length? But it cannot be that easy, can it? I am really confused!
– Meg
Mar 16, 2022 at 21:11
• The GAMLSS package for R should allow this. Mar 16, 2022 at 21:56

I don't think you will be able to do that with glm.nb() but it is possible to do this is a Bayesian framework. I think there are multiple options on how to do this, but the most straightforward for someone without experience in Bayesian modelling will be to use the brms package, which allows the user to fit models using the same model syntax as R packages like lme4.

Specifically, the model formula in the brms package would be the following, assuming you want to model Total_crashes by both traffic_count and road_length and just the overdispersion parameter by road_length:

Crash_count <- brm(bf(Total_crashes~LN(traffic_count) + road_length,
data=mydata,
family = negbinomial())


See this post for an example of someone predicting the shape/phi parameter, which is the overdispersion, using variables in a model:

https://discourse.mc-stan.org/t/negative-binomial-shape-and-phi-parameters/8120/3

That post, in combination with the vignette I link below about fitting distributional models in brms, should give you a good start on fitting the model you'd like to fit

https://cran.r-project.org/web/packages/brms/vignettes/brms_distreg.html

EDIT: As others have pointed out in the comments, there are other packages that allow for this to be modelled.

• The VGAM::negbinomial function will also fit classical versions of the models -- the point of VGAM is to allow linear predictors on multiple parameters of a base distribution. Mar 17, 2022 at 3:22