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
road_length and just the overdispersion parameter by
Crash_count <- brm(bf(Total_crashes~LN(traffic_count) + road_length,
shape ~ road_length),
family = negbinomial())
See this post for an example of someone predicting the shape/phi parameter, which is the overdispersion, using variables in a model:
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
EDIT: As others have pointed out in the comments, there are other packages that allow for this to be modelled.