# How can I resolve overdispersion in a negative binomial model?

I started what seemed like a straightforward analysis, but I've gotten stuck with overdispersion in my negative binomial model.

I would like to know which sites are different from each other in terms of number of calls. Can anyone tell me please how they would approach this? (we don't have any covariates).

Calls have been counted at each site for several nights over four years. I have aggregated counts of bird calls so that there is one mean per Site and Year; 'n' is the sample size of each mean.

Here is some example data:

site <- as.factor(rep(letters[1:11], each=4))
year <- as.factor(rep(c("2017","2018","2019","2020"),11))
calls <- c(222, 3778,11472,3642,2251,3008,41924,1718,284,29,2508,1610,
16,5,128,8,130,108,75,78,32,54,40,23,4,13,67,11,60,20,26,3,99,26,82,13,
2325,3487,12696,2849,48929,18309,34645,34625)
n <- c(10,8,7,8,12,8,7,8,4,6,7,7,9,6,7,7,9,5,7,8,8,5,7,8,6,
7,7,8,8,7,7,7,8,7,7,7,9,9,7,7,8,10,7,9)
birds <- data.frame(site,year,calls,n)

And the nb model:

require(MASS)
m1 <- glm.nb(calls ~ site, weights=n,link='log',data=birds)
summary(m1)

Which is overdispersed:

df_resid <- nrow(model.frame(m_nb1)) - length(coef(m_nb1)+1)
pearson_resid <- residuals(m_nb1, type = "pearson")
pearson_sq <- sum(pearson_resid^2)
pearson_sq / df_resid

[1] 12.59071

Any thoughts will be greatly appeciated!

First, it would be better to work with the original data, not aggregated per year. Do you have access to that? If so, how many observations per site do you have? You should of course plot the time courses per site, just to check that some sites have not had changes over time.

Then you would be able to use simple ANOVA type approaches, which are very robust to overdispersion. (You could hack up a permutation test to see whether results are appreciably different. I did so in the past and very rarely found a difference.)

Alternatively, you can look directly at contrasts for your negbin model.

However, if you want to compare all sites to each other, you need to account for a multiplicity of $${11\choose 2}=55$$ tests in either case.

Two useful functions would be MASS::confint.glm() and multcomp::glht(). For instance, you can compare the means for siteb and sitec as follows:

> require(multcomp)
> summary(glht(m1,"siteb-sitec=0"))

Simultaneous Tests for General Linear Hypotheses

Fit: glm.nb(formula = calls ~ site, data = birds, weights = n, link = "log",
init.theta = 1.257857112)

Linear Hypotheses:
Estimate Std. Error z value Pr(>|z|)
siteb - sitec == 0   2.0983     0.2364   8.877   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Adjusted p values reported -- single-step method)

If you really want to compare all sites in a pairwise manner, you can create the full combinatorics for the linfct argument to glht() as follows:

> linfct <- c(paste0("site",letters[2:11],"=0"),
+ unlist(sapply(2:10,function(ii)paste0("site",letters[ii],"-site",letters[(ii+1):11],"=0"))))
> summary(glht(m1,linfct))

Simultaneous Tests for General Linear Hypotheses

Fit: glm.nb(formula = calls ~ site, data = birds, weights = n, link = "log",
init.theta = 1.257857112)

Linear Hypotheses:
Estimate Std. Error z value Pr(>|z|)
siteb == 0          0.86749    0.21637   4.009  0.00296 **
sitec == 0         -1.23084    0.23928  -5.144  < 0.001 ***
sited == 0         -4.70713    0.22891 -20.563  < 0.001 ***
sitee == 0         -3.77533    0.22773 -16.578  < 0.001 ***
sitef == 0         -4.80076    0.23130 -20.755  < 0.001 ***
siteg == 0         -5.18820    0.23233 -22.331  < 0.001 ***
siteh == 0         -5.02060    0.22962 -21.864  < 0.001 ***
sitei == 0         -4.33171    0.22830 -18.974  < 0.001 ***
sitej == 0          0.15793    0.22124   0.714  0.99977
sitek == 0          2.04393    0.21790   9.380  < 0.001 ***
siteb - sitec == 0  2.09832    0.23638   8.877  < 0.001 ***
siteb - sited == 0  5.57461    0.22587  24.680  < 0.001 ***
...

The advantage is that glht() will account for the multiple comparisons automatically (as long as you put all your comparisons into a single call to glht(), of course - if you call it 55 times with different contrasts, you will need to do the corrections yourself).

• Thank you so much for your quick and informative response. Just to clear up my confusion though, ..so even though the nb model is overdispersed, I can still go ahead and compare site means?
– Lisa
Sep 23 '20 at 7:53
• Yes indeed. The parameter estimates are the group means in your model, and especially if you work on the original data, the asymptotics should be valid. Sep 24 '20 at 12:05