Understanding the parameters inside the Negative Binomial Distribution

Question

I was trying to fit my data into various models and figured out that the fitdistr function from library MASS of R gives me Negative Binomial as the best-fit. Now from the wiki page, the definition is given as:

NegBin(r,p) distribution describes the probability of k failures and r successes in k+r Bernoulli(p) trials with success on the last trial.

Using R to perform model fitting gives me two parameters mean and dispersion parameter. I am not understanding how to interpret these because I cannot see these parameters on the wiki page. All I can see is the following formula:

where k is the number of observations and r=0...n. Now how do I relate these with the parameters given by R? The help file does not provide much information either.

Also, just to say a few words about my experiment: In a social experiment that I was conducting, I was trying to count the number of people each user contacted in a period of 10 days. The population size was 100 for the experiment.

Now, if the model fits the Negative Binomial, I can blindly say that it follows that distribution but I really want to understand the intuitive meaning behind this. What does it mean to say that the number of people contacted by my test subjects follows a negative binomial distribution? Can someone please help clarify this?

Answer 1

You should look further down the wikipedia page, to where it says "gamma-Poisson mixture". While the definition you cite (which I call the "coin-flipping" definition since I usually define it for classes as "suppose you want to flip a coin until you get k heads") is easier to derive and makes more sense in an introductory probability or mathematical statistics context, the gamma-Poisson mixture is (in my experience) a much more generally useful way to think about the distribution in applied contexts. (In particular, this definition allows non-integer values of the dispersion/size parameter.) In this context, your dispersion parameter describes the distribution of a hypothetical Gamma distribution that underlies your data and describes unobserved variation among individuals in their intrinsic level of contact. In particular, it is the shape parameter of the Gamma, and it may be helpful in thinking about this to know that the coefficient of variation of a Gamma distribution with shape parameter $\theta$ is $1/\sqrt{\theta}$; as $\theta$ becomes large the latent variability disappears and the distribution becomes Poisson.

Answer 2

As I mentioned in my earlier post to you, I'm working on getting my head around fitting a distribution to count data also. Here's among what I've learned:

When the variance is greater than the mean, overdispersion is evident and thus the negative binomial distribution is likely appropriate. If the variance and mean are the same, the Poisson distribution is suggested, and when the variance is less than the mean, it's the binomial distribution that's recommended.

With the count data you're working on, you're using the "ecological" parameterization of the Negative Binomial function in R. Section 4.5.1.3 (Page 165) of the following freely-available book speaks to this specifically (in the context of R, no less!) and, I hope, might address some of your questions:

http://www.math.mcmaster.ca/~bolker/emdbook/book.pdf

If you come to conclude that your data are zero-truncated (i.e., the probability of 0 observations is 0), then you might want to check out the zero-truncated flavor of the NBD that's in the R VGAM package.

Here's an example of its application:

library(VGAM)

someCounts = data.frame(n = c(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16),
                     freq = c(182479,76986,44859,24315,16487,15308,5736,
                              2843,1370,1115,1127,49,100,490,106,2))

fit = vglm(n ~ 1, posnegbinomial, control = vglm.control(maxit = 1000), weights=freq,
           data=someCounts)

Coef(fit)

pdf2 = dposnegbin(x=with(someCounts, n), munb=0.8344248, size=0.4086801)

print( with(someCounts, cbind(n, freq, fitted=pdf2*sum(freq))), dig=9)

I hope this is helpful.