# Is there any R packages allow direct MLE estimation of dispersion in negative binomial distribution?

Using the built-in function, I can get

data=rnbinom(1000,mu=0.1,size=0.5)

MASS::fitdistr(test, dnbinom, list(mu=0.1, size=0.5),
method = "L-BFGS-B",lower = c(0,0),upper = c(Inf,Inf))

mu           size
0.085004342   0.867696303
(0.009659642) (0.618157097)


Is there any way to directly estimate the dispersion parameter rather than size parameter?

A follow-up question:

We can observe that the standard error of the MLE for size is large. This is quite often seen when the mean is small, so it's hard for the algorithm to find the global maximum. Are there any methods to overcome this limitation especially when mean is extremely small (e.g., < 0.01).

• The size parameter is the same as the dispersion parameter. Oct 27, 2021 at 3:07
• I suspect the imprecision of the estimates of the dispersion parameter when the mean is very small is largely unavoidable. The variance of the NB = $\mu + \mu^2/k$, where $k$ is the dispersion parameter; at $\mu = 0.01$ we have $0.01 + 0.0001/k$. Choosing between, say, $k=1$ and $k=2$ is going to make only a tiny difference in the distribution, so naturally $k$ will be hard to estimate precisely. Oct 27, 2021 at 3:20
• Sorry if it causes any misunderstanding. The definition for dispersion varies in different resources. I am referring to the dispersion (alpha) where variance=mean+mean^2*alpha. Oct 28, 2021 at 2:51
• Then the MLE of $\alpha$ = $1/$ the MLE of size. If you really care, you can just add an extra line of code to invert the estimate and confidence intervals yourself. Oct 28, 2021 at 3:03
• Thanks. I understand that when the mean is small, it is naturally hard to estimate the k precisely. However, do you happen to know if any method to provide better estimate in this case? I heard some quasi-likelihood methods can do this job. Oct 30, 2021 at 4:39

It appears that at least some of your difficulty arises not because the MLE of the dispersion parameter is relatively inefficient in smallish samples, but because the distribution of the MLE is not well-approximated by the Gaussian distribution when the mean of the data's distribution is small, even with sample sizes that appear quite large.

For example, we simulate 1000 draws from a Negative Binomial distribution with $$\mu=0.05$$ and $$\alpha = 2$$, calculate $$\hat{\alpha}$$, and replicate 10,000 times. The code and histogram follow:

de <- rep(NA, 10000)
true_mu <- 0.05
for (i in seq_along(de)) {
test <- rnbinom(1000, mu=true_mu, size=0.5)

tryCatch({
tmp <- MASS::fitdistr(test, dnbinom, list(mu=true_mu, size=0.5),
method = "L-BFGS-B",lower = c(0.001,0.01),upper = c(1,100))
de[i] <- 1/tmp$estimate[2] }, error = function(e) {}) } hist(de, xlab="MLE of dispersion parameter (true value = 2)", main="Histogram of MLE values")  In a case like this, calculating confidence intervals for the dispersion parameter based on the asymptotic distribution will produce misleading results: test <- rnbinom(1000, mu=true_mu, size=0.5) tmp <- MASS::fitdistr(test, dnbinom, list(mu=true_mu, size=0.5), + method = "L-BFGS-B",lower = c(0.001,0.001), upper = c(1,100)) > tmp mu size 0.057006143 0.425602884 (0.008037687) (0.288389381) >  As we can see, using the asymptotic Gaussian distribution to construct, say, 95% confidence intervals would lead to a zero lower bound for size (assuming we were paying attention to the fact that size and $$\alpha > 0$$.) The actual 95% bounds for this example are $$(1.01, \infty)$$, which is less helpful than it might be. Bootstrapping the sample gives us much better confidence intervals: # Returns estimate of dispersion; if the estimate fails, resamples the data dispersion_estimate <- function(df, indices) { res <- NA while (is.na(res)) { res <- tryCatch( 1 / MASS::fitdistr(df[indices], dnbinom, list(mu=true_mu, size=0.5), method = "L-BFGS-B",lower = c(0.001,0.001), upper = c(1,100))$estimate[2]
, error = function(e) {
indices <<- sample(1:length(df), replace=TRUE)
NA
}
)
}
res
}

library(boot)
#  test <- rnbinom(1000, mu=true_mu, size=0.5) (use same data as before)
res <- boot(test, dispersion_estimate, 10000)
boot.ci(res)

Intervals :
Level      Normal              Basic
95%   (-0.635,  5.202 )   (-1.138,  4.662 )

Level     Percentile            BCa
95%   ( 0.046,  5.846 )   ( 0.167,  6.693 )
Calculations and Intervals on Original Scale


The percentile and BCa estimates are the ones we would prefer.

As we can see, the confidence intervals are still quite wide relative to the true value of 2. This is not because fitdistr is having a hard time finding the global maximum, it's because the likelihood function is quite flat in the vicinity of the maximum, so the spread of reasonable values really is large. We show this by plotting $$2*\log \mathcal{L(\alpha | \mu)}$$, the profile log-likelihood for the dispersion parameter given the true $$\mu$$, and calculating a 95% confidence interval based on the asymptotic $$\chi^2(1)$$ distribution of this statistic:

llf <- function(dispersion) {
sum(dnbinom(test, mu=true_mu, size=1/dispersion, log=TRUE))
}

llf_values <- rep(0,1000)

dispersion_values <- seq(0.1,10,length.out=1000)

for (i in seq_along(dispersion_values)) {
llf_values[i] <- llf(dispersion_values[i])
}

plot(2*llf_values~dispersion_values, main="Log likelihood function",
xlab="Dispersion parameter", ylab="Log likelihood")
abline(h = max(2*llf_values) - qchisq(0.95,1), lwd=2, col=2, lty=2)


Although the plot doesn't look particularly flat, the y-axis scale reveals the truth - it just doesn't make much difference to the log-likelihood whether the dispersion parameter equals, say, 1 or 2 or 5, even with a sample size of 1000.

• Hi jbowman, thanks for your valuable input. I also tested some examples and found the likelihood function for the dispersion parameter given the true mean is flat when the mean is small. Nov 4, 2021 at 4:19