Compare glm.nb vs glm(..., negative.binomial(k), ..) models

Question

My $Y$ is a count variable and I am estimating it using a negative binomial function. I am not sure how I can decide which model is better. In the first model I am selecting a value for $k$, and the second the MASS package is estimating a value of $k$ for me.

Will

anova(mod1, mod2, test="Chisq") 

work?

If so, which model is the null referring to?

Otherwise, is there a better way to decide which model is better?

library(MASS)
mod1<- glm(Y ~ X1 + X2 + X3, negative.binomial(1), myData)
mod2<- glm.nb(Y ~ X1 + X2 + X3, myData)

Answer 1

The chi-squared test is valid here, it is testing the hypothesis that

$$H_0: \phi = 1$$ with $\phi = \frac{1}{\theta}$ the overdispersion parameter. It basically performs a likelihood-ratio test to compare both models (see ?anova). This is valid since your model with fixed overdispersion parameter is nested in the one where the overdispersion is freely estimated. The test statistic will follow a $\chi_1^2$ distribution, since your extended model has one additional estimated parameter. If you are unsure, you can always run some simulations under $H_0$ to prove this null-distribution.

Note that the null-distribution would be an even mixture between a point mass at 0 and $\chi_1^2$ if you were to test $$H_0: \phi = 0$$, since this would put the overdispersion at the edge of the parameter space. However, this does not apply here since $\phi$ can shrink to values smaller than 1.

Answer 2

ANOVA is used to compare nested models, i.e., models $M_1$ and $M_2$ where all predictors that appear in $M_1$ also appear in $M_2$, but $M_2$ contains additional ones. ANOVA then answers the question whether the additional predictors explain more variance than we would expect by chance alone.

Of course, one could also see this as an example of constrained model fitting: we could see $M_1$ as containing the same predictors as $M_2$, but with the parameter estimates of the additional predictors constrained to be zero.

Your use case is different, even if it is also a case of constrained parameter estimates. I don't know whether ANOVA would be valid here. (You could run a few simulations.)

I would recommend an approach using cross-validation (CV) and scoring-rules as follows. For each of the $k$ folds of your CV:

• Fit both negative binomial (NB) regression models to the in-sample data.
• Predict conditional NB distributions for your out-of-sample data.
• Assess the quality of these conditional distributions for the holdout using proper scoring rules. If you want to use the Brier or the spherical score, you will need to calculate the sum of squared NB probability masses; see here on how to do this.

Repeat this process $k$ times and average the scores. Go with the model that yields the smallest score.

A few points to keep in mind:

• I'd recommend running the entire CV multiple times with different seeds for the random number generator. See how consistently one model outperforms the other in terms of the score.
• Similarly, use multiple different scores to see how robust your results are. The Brier, spherical and logarithmic scores are probably the most common ones. Merkle & Steyvers (2013, Decision Analysis) may be helpful.
• The approach above does not consider parameter estimation uncertainty - strictly speaking, the predictions should not use the NB distribution, but an (even) more overdispersed one to account for this. Compare how we use the $t$ distribution for predicting in an OLS context, even if we assume normally distributed residuals. If you want to get fancy, you could bootstrap within each CV fold to account for this uncertainty, but I don't know whether this is worth the additional effort.