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

Question

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).

Answer 1

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.