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.

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.
