using rnegbin (plot negative binomial distribution based on real data)

Question

Say you have data with mean $\mu$ and standard deviation $\sigma$. You think they came from a negative binomial ($\sigma > \mu$), and you want to simulate a negative binomial distribution based on those parameters.

The MASS package in R has the handy function rnegbin(n, mu, theta) that simulates $n$ samples from a negative binomial distribution with mean $\mu$, Gamma-shape parameter $\theta$ and variance $\mu + \frac{\mu^2}{\theta}$.

Cool. The mean of the observed sample fits right in. But what about $\sigma$? I imagine it should play with $\theta$, but I don't see how.

Answer 1

You said you already know $\mu$ and $\sigma$, and you're trying to find $\theta$, so you're solving the equation

$$\sigma^2 = \mu + \mu^2/\theta$$ for $\theta$ with everything else given.

This is elementary algebraic manipulation. I'm not sure I see the difficulty in isolating $\theta$ here.

However, moment-matching is not necessarily a particularly efficient way to estimate parameters in the negative binomial. If you have all the data (not just the mean and standard deviation), why would you not use maximum likelihood to get $\theta$?

This is straightforward in R in several ways. One is to use

 library(MASS)  #this comes with R but is not loaded by default
 summary(glm.nb(y~1))

Here's an example on a set of data:

   7  3  7  7  2  6  3 12 12  8  4  6 26 21  8 15 12 13  2 11  0 27  
   7  8  5 13  2 12  4  1  5  6 17  2  3  6 12 12  2 12

Using the moment-matching approach you get $\hat{\mu}=8.525$ and $\hat{\theta}=2.268387$. The ML estimate of the mean parameter will be the same but the estimate of $\theta$ is not:

> nbfity = glm.nb(y~1)
> summary(nbfity)

[snip a few lines]

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)   2.1430     0.1168   18.35   <2e-16 ***

(Dispersion parameter for Negative Binomial(2.3352) family taken to be 1)

    Null deviance: 43.291  on 39  degrees of freedom
Residual deviance: 43.291  on 39  degrees of freedom
AIC: 251.26

Number of Fisher Scoring iterations: 1


              Theta:  2.335 
          Std. Err.:  0.672 

 2 x log-likelihood:  -247.257 

So $\hat{\theta}=2.335$ and you get a standard error for free. The estimate of $\mu$ is obtained by exponentiating the coefficient of the intercept:

  exp(nbfity$coefficients[1])
  [1] 8.525

One can get a standard error on this scale if needed using predict.

---

Note that besides MASS::rnegbin there's stats::rnbinom which uses another parameterization (which the Wikipedia article on the negative binomial helpfully gives the mean and variance in terms of).