Estimate negative binomial dispersion parameter $k$ using mean and proportion of zeros I came across supplemental methods of a paper estimating the mean ($R$) and dispersion ($k$) of a negative binomial distribution that stated:

Page 8: "Given estimates of the mean ($\hat{R}$) and proportion of zeros ($\hat{p_0}$) of a negative binomial distribution, the parameter $k$ can be estimated by solving the equation $\hat{p_0} = (1+\frac{\hat{R}}{k})^{-k}$ numerically."

This "zero-class estimator" approach is also used here.
I would like to compare the accuracy of this method (with regards to inference of $k$) with my results using MLE. This initially seemed simple, but I have been unable to successfully estimate $k$ using this approach in R programming. I tried solving for $k$ algebraically but my algebra may be wrong (happy to post if requested, but omitting equations for readability/brevity).
Any advice on how to use this approach to estimate $k$ (in R or other statistical software) would be much appreciated.
 A: Since $(1+\frac{\hat{R}}{k})^{-k}$ is a decreasing function of $k$, the function $\hat{p}_0 -(1+\frac{\hat{R}}{k})^{-k}$ will have at most one zero.  However, when the lower asymptote of $(1+\frac{\hat{R}}{k})^{-k}$ is larger than $\hat{p}_0$, no zero will be found.  Here's some simple R code
set.seed(555)
x = rnbinom(100,size=.5,prob=.4)
p0 = mean(x==0)
mu = mean(x)
f = function(k,p0,mu){
    return(p0 - (1+mu/k)^(-k))
}
uniroot(f,c(.001,1000),p0=p0,mu=mu)

```

A: The exact answer can be given by the Lambert $W$ function as follows:
$$k=-\frac{\hat{R}\ln(\hat{p}_0)}{\ln(\hat{p}_0)-\hat{R}\,W\!\left(\frac{\hat{p}_0^{1/{\hat{R}}}\,\ln(\hat{p}_0)}{\hat{R}}\right)},$$
where $W$ is the Lambert $W$ function, which is the inverse of the function $f(x)=x\,e^{x}.$ In R, you can install the lamW package, which gives access to the lambertWm1 function (a better branch of the $W$ function to use; see Henry's comment), and then code up the equation above as follows:
install.packages('lamW')
require(lamW)
k = - R * log(p) / (log(p) - R * lambertWm1(p^(1/R) * log(p) / R))

As mentioned in the comments, this does require $p>0$ to satisfy the domain requirements of the logarithm function.
A: Given that both $\hat{p}_0$ and $(1+\frac{\hat{R}}{k})^{-k}$ are by definition non-negative, you can just minimize $(\hat{p}_0 - (1+\frac{\hat{R}}{k})^{-k})^2$. Almost any software will have some generic minimization tools (e.g. ucminf in R). The only difficulty will be the scenarios where the estimate of the dispersion parameter should be 0 or $\infty$ (depending on your parametrization), which is a genuine potential issue.
