An impossible estimation problem? Question
The variance of a negative binomial (NB) distribution is always greater than its mean. When the mean of a sample is greater than its variance, trying to fit the parameters of a NB with maximum likelihood or with moment estimation will fail (there is no solution with finite parameters).
However, it is possible that a sample taken from a NB distribution has mean greater than variance. Here is a reproducible example in R.
set.seed(167)
x = rnbinom(100, size=3.2, prob=.8);
mean(x) # 0.82
var(x) # 0.8157576

There is a non-zero probability that the NB will produce a sample for which parameters cannot be estimated (by maximum likelihood and moment methods).

*

*Can decent estimates be given for this sample?

*What does estimation theory say when estimators are not defined for all samples?

About the answer
The answers of @MarkRobinson and @Yves made me realize that parametrization is the main issue. The probability density of the NB is usually written as
$$P(X = k) = \frac{\Gamma(r+k)}{\Gamma(r)k!}(1-p)^rp^k$$
or as
$$P(X = k) = \frac{\Gamma(r+k)}{\Gamma(r)k!} \left(\frac{r}{r+m}\right)^r \left(\frac{m}{r+m}\right)^k.$$
Under the first parametrization, the maximum likelihood estimate is $(\infty, 0)$ whenever the variance of the sample is smaller than the mean, so nothing useful can be said about $p$. Under the second, it is $(\infty, \bar{x})$, so we can give a reasonable estimate of $m$. Finally,  @MarkRobinson shows that we can solve the problem of infinite values
by using $\frac{r}{1+r}$ instead of $r$.
In conclusion, there is nothing fundamentally wrong with this estimation problem, except that you cannot always give meaningful interpretations of $r$ and $p$ for every sample. To be fair, the ideas are present in both answers. I chose that of @MarkRobinson as the correct one for the complements he gives.
 A: In the Negative Binomial (NB) example, the likelihood may have its
maximum at an infinite distance for $p \to 0$ and $r \to \infty$, on
the boundary of the domain $\Theta := (0,\,1)\times(0,\,\infty)$. If it
turns out that the Poisson distribution leads for some mean $\lambda >0$
to a likelihood which is greater than the NB, then the likelihood can
increase when $[p,\,r] \in \Theta$ moves along a path with $p \to 0$,
$r \to \infty$ and $rp/(1-p) \to \lambda$. The probability that
the maximum likelihood is found on the boundary is not zero.
A similar problem but with simpler diagnostic is for the
Lomax distribution:
it is known that the ML estimate of the shape is infinite when the
sample has coefficient of variation $\text{CV} < 1$. Yet the
probability of this event is positive for any sample size, and is
for example $>0.3$ for $\alpha = 20$ and $n = 200$.
ML properties are for a large sample size: under
regularity conditions, a ML estimate is shown to exist, to be unique
and to tend to the true parameter. Yet for a given finite
sample size, the ML estimate can fail to exist in the domain,
e.g. because the maximum is reached on the boundary.  It can also
exist in a domain which is larger than the one used for maximisation.
In the Lomax example, some people would chose to use the exponential
distribution, which is the limit for $\alpha \to \infty$ and $\lambda
/ \alpha \to \theta >0$. This boils down to accepting an infinite ML
estimate. Since the Lomax is a special re-parameterisation of the
two-parameter Generalised Pareto Distribution $\text{GPD}(\sigma,\,\xi)$
with shape $\xi >0$, we could as
well fit a GPD, then finding $\widehat{\xi} < 0$ instead of the exponential
$\widehat{\xi} = 0$. For the NB example, we can chose
to fit a Poisson distribution thus accepting a boundary value of the
NB parameter.
For the sake of invariance by re-parameterization, I believe that
infinite parameters can make sense in some cases. 
A: Basically, for your sample, the estimate of the size parameter is on the boundary of the parameter space. One could also consider a reparameterization such as d = size / (size+1); when size=0, d=0, when size tends to infinity, d approaches 1.  It turns out that, for the parameter settings you have given, size estimates of infinity (d close to 1) happen about 13% of the time for Cox-Reid adjusted profile likelihood (APL) estimates, which is an alternative to MLE estimates for NB (example shown here).  The estimates of the mean parameter (or 'prob') seem to be ok (see figure, blue lines are the true values, red dot is the estimate for your seed=167 sample).  More details on the APL theory are here.
So, I would say to 1.: Decent parameter estimates can be had .. size=infinity or dispersion=0 is a reasonable estimate given the sample.  Consider a different parameter space and the estimates will be finite.
