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.

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)

```

