# Rule of thumb for deciding between Poisson and negative binomal models

I am analyzing the number of a specific class of mutations in cancer genome sequencing data, as a function of: (1) overall mutation rate per sample, (2) length of the genomic segment (gene) being analyzed, and (3) a particular covariate of interest that applies to a small number of genes. Not surprisingly, the vast majority of genes, in any sample, have zero mutations - but I don't think modeling this as a "zero-inflated" process makes sense, because there's no reason to think that there is a separate zero-generating process at work, it's just that most genes by chance don't end up getting hit. The most straightforward models would seem to be Poisson and negative binomial. The latter, overall, seems to be favored in the literature, but it is much more computationally intensive: if I run a simple model looking at all samples in my dataset (~3e7 observations), considering mutation count as a function only of sample-specific mutation rate and gene length, it takes 2 min to return results with Poissson (glm(..., family='poisson')), but 2 hours to return results with negative binomial (glm.nb(...)). The results are similar: in both cases the two explanatory variables are positively correlated with mutation count, at p<2e-16. The NB model has a modestly (~10%) lower AIC value, but my understanding is that one can't compare AIC between classes of model.

My understanding is that NB models are preferred when the variance of the measured variable exceeds its mean. In my case, these values are 0.110 and 0.023, respectively. Is there any sort of rule of thumb as to how much the variance should exceed the mean, before going to the effort of running a negative binomial model, or better yet some test that can decide how much more inaccurate a Poisson model would be, given a set of data, compared to NB?

• re: variance exceeding the mean, this is conditional variance/mean, which you can't just calculate via "mean(y)" or "var(y)". you may have better luck with pseudo-poisson regression in terms of execution time. Commented Nov 6, 2023 at 16:58
• Thank you for this correction, John. Is there a simple way to calculate the conditional variance/mean? Also, is pseudo-poisson the same as quasi-poisson, discussed by PBulls below? Commented Nov 10, 2023 at 20:14
• One option is to use the tests that Pbull recommends to assess overdispersion. Personally though, whenever possible, I like to evaluate assumptions graphically. In this case, we could plot predictive mean versus (absolute) residual (from a poisson regression), and add lines to the plot indicating the one, two and three theoretical standard devs. And yes, theyre the same, and quasi is correct, not pseudo (misremembered). Commented Nov 10, 2023 at 20:21

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

• Thanks for a detailed explanation and worked example. It does seem like the Pearson dispersion statistic makes the most sense as a simple way to analyze the data. I'm not sure how one would compare the Poisson vs NB models using the score test, however. In the worked model, the t test of poi_glm gives p=4.6e-9. If I use glm.nb to analyze the same data, I get a t test p=2.5e-8. I'm not sure how to compare these. If I run a two-sample t test on the NB vs Poisson z values, I get p=0.90.) Whether meaningful or not, the Pearson statistic of the NB model is 0.94, vs 10.5 for Poisson. Commented Nov 10, 2023 at 22:22
• These tests are only for Poisson fits, to check whether mean and variance are indeed the same. In a NB model the variance is a second parameter so they do not apply. I think you can still use the Pearson statistic for NB as you do for Poisson, there are other models (e.g. NB with parametrized exponent, Poisson-inverse-Gaussian) that can handle even more overdispersion but 0.94 seems quite OK. FYI, the data were taken from the glm.nb example so I would assume a NB model is appropriate here. I only just learned about the DHARMa package which may also be useful for goodness of fit. Commented Nov 10, 2023 at 22:34