Paraphrasing Hilbe (2014), the standard Poisson model assumes equidispersion or mean and variance sharing a single parameter. In practice this is not usually the case, with most data exhibiting overdispersion or more variance than you would expect from the mean. As already mentioned in a comment dispersion is also conditional, it may be caused by failure to include relevant or appropriately scaled predictors or interactions, or missing-not-at-random observations. In either case, an extradispersed model (also underdispersed) will not provide correct estimates of standard errors and may render your inference invalid.
Pearson dispersion
The most straightforward indication for extradispersion is the Pearson dispersion statistic, which is the sum of squared Pearson residuals divided by the residual degrees of freedom - see also the worked example below. Note that you should not use the deviance residuals or its dispersion statistic as these are biased (see also here). Equidispersion implies a value of 1 for the Pearson dispersion statistic. Hilbe proposes a threshold of 1.25 for 'moderate' or as low as 1.05 for 'large' numbers of observations to indicate problematic overdispersion (and presumably 0.8-0.95 would mean the same for underdispersion). I would put $n\approx3e7$ in the 'large' category.
Tests for overdispersion
There are also tests you can conduct, one such is a score test that directly compares the Poisson and negative binomial distributions from $z_i=\big((y_i-\mu_i)^2-y_i\big)/(\mu_i\sqrt2)$ with $y_i$ each observed response and $\mu_i$ its predicted mean. The Poisson distribution implies $z\sim N(0,1)$ so a one-sample $t$ test can provide a P-value for testing Poisson vs. negative binomial.
Another test for equidispersion is the Lagrange Multiplier $\big(\sum(\mu_i^2)-n\bar y\big)^2/(2\sum\mu_i^2)$ which follows a one-degree $\chi^2$ distribution under the null. Finally, you could conduct a contingency table test (Fisher's, chi-square) on the observed counts versus those predicted from the model, the argument being that as long as both quantities match sufficiently the model fit is adequate.
I'm going to make a brief aside here on comparing models, feel free to skip this paragraph. AIC is only valid if the (log-)likelihood is calculated consistently (see e.g. here), so it would depend on the full likelihood being calculated in the first place and the exact implementation next. Skimming the MASS::glm.nb source code it seems to call stats::glm.fit() by default, so it should be appropriate to compare your models via AIC. Another option may be the Vuong or closely related Clarke test, though I cannot say for sure whether Poisson and negative binomial models are in fact partially non-nested (which would make such test invalid).
Handling overdispersion
By now you've almost certainly established that your data are not likely to be equidispersed. A possible solution to address extradispersion is to fit a negative binomial distribution, which adds an extra parameter that uncouples the variance from the mean. However, that is not the only option.
A very straightforward alternative is to fit a quasi-Poisson model. The underlying idea is to fit your model assuming the observed Pearson dispersion instead of equidispersion. R provides the quasipoisson GLM family for this (which might in fact be fitting 'quasi-likelihood'? Not sure). You can also simply scale your Poisson GLM's covariance (standard errors) by the the Pearson dispersion (its square root).
A second option is to apply a robust (aka. sandwich, empirical) variance estimator. These are usually intended for correlated data, but they do a decent job at handling overdispersed data as well. R's sandwich is your friend here.
Third, bootstrapping your standard errors can also be used to relax distributional assumptions in your model, but is probably computationally less attractive than any of the other options if fitting even a single model takes that much time.
Finally, while all of these are alternatives to the standard Poisson GLM for extradispersed data, I consider it not too likely that your actual conclusion will meaningfully change in the case where both the Poisson and negative binomial GLMs give you $P<2e{-}16$. The Poisson model may estimate P too low, but I assume there's quite a few orders of magnitude to go before the interpretation changes (without considering any multiplicity or other possible issues such as the zero-inflation).
Worked example
Checking extradispersion in a Poisson GLM:
data("quine", package="MASS")
poi_glm <- glm(Days ~ Sex/(Age + Eth*Lrn), data = quine, family="poisson")
## Pearson dispersion statistic
pearson_dispersion <- sum(residuals(poi_glm, type="pearson")**2) / poi_glm$df.residual
## Score test
y <- quine$Days
mu <- predict(poi_glm, quine, type="response")
z <- ((y - mu)**2 - y) / (mu * sqrt(2))
t.test(z)
## Lagrange Multiplier test
lagrange <- (sum(mu**2) - length(mu)*mean(mu))**2 / (2*sum(mu**2))
pchisq(lagrange, 1, lower.tail = FALSE)
## Test of observed vs. predicted counts
y_range <- seq(0, max(y))
observed <- vapply(y_values, \(yi) sum(y == yi), numeric(1))
predicted <- vapply(y_values, \(yi) {
round(sum(dpois(yi, fitted(poi_glm))))
}, numeric(1))
fisher.test(observed, predicted, simulate.p.value = TRUE)
Some alternative standard error estimators:
## Quasi-count model
qpoi_glm <- glm(Days ~ Sex/(Age + Eth*Lrn), data = quine, family="quasipoisson")
## Could also have updated Poisson SEs manually [multiply by sqrt(dispersion)]
all.equal(vcov(qpoi_glm), vcov(poi_glm) * pearson_dispersion, tolerance = 1E-5)
## Robust variance-covariance estimator
robust_se <- sqrt(diag(sandwich::vcovHC(poi_glm)))
## Bootstrap
bootfun <- function(d, i) {
sqrt(diag(vcov(glm(Days ~ Sex/(Age + Eth*Lrn), data = d[i,], family="poisson"))))
}
boot_samples <- boot::boot(quine, bootfun, R=1000)
boot_se <- colMeans(boot_samples$t)
Reference
Modeling Count Data (2014), Joseph M. Hilbe, Cambridge University Press
