How to statistically test if the dispersion parameter in a negative binomial distribution is different? I have a negative binomial distribution for two species. I have been able to fit (dummy) data to estimate mean (mu) and dispersion parameter (size). However, I want to statistically test if the two distributions are statistically different specifically in their dispersion parameter. Any ideas how to go about testing this?
Here is the dummy data, code in R, and graphical output to aid in understanding the question.
library(ggplot2)
library(MASS)
set.seed(111)

df <- data.frame(
  count = rnbinom(500, rep(c(5, 10), each  = 250), 0.5),
  species = rep(c("A", 'B'), each = 250)
)

ests <- sapply(split(df$count, df$species), function(x) {
  est <- fitdistr(x, "negative binomial", method = "SANN")$estimate
  formatted <- paste0(names(est)\[1\], " = ", format(est, digits = 2)\[1\], ",",
                      names(est)\[2\], " = ", format(est, digits = 2)\[2\])
  formatted
})

mybinwidth <- 1

spec_A = df\[df$species == "A",\]
spec_B = df\[df$species == "B",\]

ggplot(df, aes(count)) +
  geom_histogram(binwidth = mybinwidth,
                 aes(fill = species), alpha = 0.5,
                 position = "identity") +
  stat_function(aes(color = "A"), 
                data = data.frame(species = "A"),
                fun = function(x, size, mu) {
                  mybinwidth * nrow(spec_A) * dnbinom(x,size = size, mu = mu)
                },
                args = fitdistr(spec_A$count, "negative binomial", method="SANN")$estimate, 
                xlim = c(0, max(df$count)), n= max(df$count) + 1, inherit.aes = FALSE) +
  stat_function(aes(color = "B"), 
                data = data.frame(species = "B"),
                fun = function(x, size, mu) {
                  mybinwidth * nrow(spec_B) * dnbinom(x,size = size, mu = mu)][1]][1]


 A: gamlss model with NBI parametrization
Can you try using the gamlss() function in the gamlss package? Something like this:
library(gamlss)

mI <- gamlss(count ~ species, 
             sigma.formula = ~ species, 
             family = NBI, 
             mu.link = "log", 
             sigma.link = "log",
             data = df) # fits the model 

plot(mI) # plots the model diagnostics

summary(mI) # summarizes the model 

Note that both the mean mu (aka $\mu$) and the dispersion sigma (aka $\sigma$) are modelled on the log scale and are allowed to differ across the two species.
Mu Coefficients
When you look at the model summary, you will see a portion called Mu Coefficients.  This portion provides you estimates for the parameters $\beta_0$ and $\beta_1$ in the equations:
$log(\mu) = \beta_0$ for species A
$log(\mu) = \beta_0 + \beta_1$ for species B
You can find the estimated value of $\beta_0$ at the intersection of the Intercept row and the Estimate column and the estimated value of $\beta_1$ at the intersection of the speciesB row and the Estimate column under the Mu Coefficients.
Note that the value of $\mu$ for each species can be obtained by exponentiating the above equations:
$\mu = exp(\beta_0)$ for species A
$\mu = exp(\beta_0 + \beta_1)$ for species B
To get the estimated value of $\mu$ for each species, simply replace $\beta_0$ and $\beta_1$ with their estimated values in these last equations.
In the portion of model summary output for Mu Coefficients, you will find the p-value for the test comparing the following competing hypotheses:
$Ho: \beta_1 = 0$ versus $Ha: \beta_1 \neq 0$
Just look for it at the intersection of the speciesB row and Pr(>|t|) column of the Mu Coefficients output.
This test will help you determine whether the log-transformed value of the mean $\mu$ is different between species B and species A.
Sigma Coefficients
The portion called Sigma Coefficients of the model summary provides you with estimates for the parameters $\gamma_0$ and $\gamma_1$ in the equations:
$log(\sigma) = \gamma_0$ for species A
$log(\sigma) = \gamma_0 + \gamma_1$ for species B
You can find the estimated value of $\gamma_0$ at the intersection of the Intercept row and the Estimate column and the estimated value of $\gamma_1$ at the intersection of the speciesB row and the Estimate column under the Sigma Coefficients.
The value of $\sigma$ for each species can be obtained by exponentiating the last two equations above:
$\sigma = exp(\gamma_0)$ for species A
$\sigma = exp(\gamma_0 + \gamma_1)$ for species B
To get the estimated value of $\sigma$ for each species, simply replace $\gamma_0$ and $\gamma_1$ with their estimated values in these last two equations.
In the portion of model summary output for Sigma Coefficients, you will find the p-value for the test comparing the following competing hypotheses:
$Ho: \gamma_1 = 0$ versus $Ha: \gamma_1 \neq 0$
The p-value can be found at the intersection of the speciesB row and Pr(>|t|) column of the Sigma Coefficients output.
I believe this last test is the test you want, as $\sigma$ is a dispersion parameter. The test compares the log-transformed value of the dispersion parameter $\sigma$ between species B (non-reference species) and species A (reference species).
Note
The negative binomial distribution has two different parametrizations in gamlss: NBI and NBII.  Both of these parametrizations use a mu parameter and a sigma parameter and use log links for modelling mu and sigma, respectively.  However, the results produced for sigma by the two parametrizations will be different.  For this reason, if you choose to work with one particular parametrization among NBI and NBII, you must specify which parametrization you chose.
