What statistical test should i use for my count data with extreme overdispersion?

Question

I want to test the number of varroa mites in individual apiaries is higher in the south of England vs the north. I'm working with count data (mites) and categorical data (location). I used a binominal GLM, however the output produces a EXTREMELY high overdispersion (around 1600) is there any way to handle this? Or should I be using a different statistical test? I used this code:

glm(formula = Varroa ~ location, family = poisson)

This is the subsequent output the two locations are south and north:

Coefficients: 
(Intercept) locationsouth 7.73165 0.09428 Degrees of Freedom: 1999 Total (i.e. Null); 1998 Residual Null Deviance:  3387000 Residual Deviance: 3380000 AIC: 3394000 –

Answer 1

There are two sorts of issues here: The treatment of your independent variable and the model chosen.

For your IV (location) as people noted in the comments, using it as a categorical variable is not usually great. If you only have two locations, then it's fine, but if you have more locations, you will want to look into other ways to treat it. These might be based on why you think there are different numbers in the north and south. E.g. if you think it is due to hours of daylight, then you could use latitude; if you think it is due to temperature, you could use average temperature in a location and so on. Even if you don't have a particular reason, you will want to do something other than a purely categorical variable. One choice is to use latitude and longitude.

Then there is your model. If you have a count dependent variable, you want a count regression. The usual starting place is Poisson regression, but overdispersion is very common (I've never had a data set that didn't have overdispersion). The usual solution there is a negative binomial regression. There are also zero-inflated versions of these models, if you have a lot of sites with no mites.

Answer 2

Peter is correct that negative binomial regressions are typically used to address dispersion (in fact it is probably the default option for most), but this is usually in the case where overdispersion is not very extreme such as yours. My guess is that you have a response which has a low mean (around 1 or 2), high variance, and a heavy right skew. This typically contributes to this level of dispersion (barring no issues with zeroes like Peter noted), as the variance of the counts can only meaningfully fluctuate near the base of the distribution, whereafter it varies to a very minor degree for higher values. Given your concern, I think this would actually be a better candidate for a Poisson inverse Gaussian (PIG) model, which is tailored to handling these sorts of problems. That is not a guarantee, only an option that may potentially be helpful.

Here is how that is accomplished. The negative binomial (NB) regression is technically a mixture model, in that it uses a Poisson (discrete) distribution along with a Gamma (continuous) distribution to model the response. This gives the model it's flexibility in dealing with dispersion. The PIG model is similar in that it uses a Poisson and an inverse Gaussian distribution (the IG distribution specifically because of it's highly right-skewed nature), and models a variance of $\mu + \alpha \mu^3$ rather than $\mu + \alpha \mu^2$. If zero-inflation is indeed an issue, one can just use a zero-inflated PIG (ZIPIG) to deal with both concerns. The gamlss package can fit both models with relative ease.

An example can be found below, where the PIG model lowers the dispersion by a little but not by a crazy amount. In some models this may vary more:

#### Load Libraries and Data ####
library(gamlss) # for model 
library(COUNT) # for data and vcov 
library(msme) # for dispersion checks
data(rwm5yr) # data for model
rwm1984 <- subset(rwm5yr, year==1984) # only use this year 

#### Poisson Model ####
pois.mod <- glm(docvis ~ outwork + age, data=rwm1984, family = poisson)
P__disp(pois.mod) # dispersion = 11.343

#### NB2 Model ####
nb.mod <- glm.nb(docvis ~ outwork + age, data=rwm1984)
summary(nb.mod)
P__disp(nb.mod) # dispersion near 1.410

#### PIG Model ####
pig.mod <- gamlss(docvis ~ outwork + age, data=rwm1984, family=PIG) 
summary(pig.mod) # sigma dispersion = 1.344

PIGs can also utilize zero-truncated and hurdle models. Modeling With Count Data, while not a perfect book, does a good job of at least explaining some of the differences between PIG models.