What is theta in a negative binomial regression fitted with R?

Question

I've got a question concerning a negative binomial regression: Suppose that you have the following commands:

require(MASS)
attach(cars)
mod.NB<-glm.nb(dist~speed)
summary(mod.NB)
detach(cars)

(Note that cars is a dataset which is available in R, and I don't really care if this model makes sense.)

What I'd like to know is: How can I interpret the variable theta (as returned at the bottom of a call to summary). Is this the shape parameter of the negbin distribution and is it possible to interpret it as a measure of skewness?

Answer 1

Yes, theta is the shape parameter of the negative binomial distribution, and no, you cannot really interpret it as a measure of skewness. More precisely:

• skewness will depend on the value of theta, but also on the mean
• there is no value of theta that will guarantee you lack of skew

If I did not mess it up, in the mu/theta parametrization used in negative binomial regression, the skewness is

$$ {\rm Skew}(NB) = \frac{\theta+2\mu}{\sqrt{\theta\mu(\theta+\mu)}} = \frac{1 + 2\frac{\mu}{\theta}}{\sqrt{\mu(1+\frac{\mu}{\theta})}} $$

In this context, $\theta$ is usually interpreted as a measure of overdispersion with respect to the Poisson distribution. The variance of the negative binomial is $\mu + \mu^2/\theta$, so $\theta$ really controls the excess variability compared to Poisson (which would be $\mu$), and not the skew.

Answer 2

I was referred to this site by one of my students in my Modeling Count Data course. There seems to be a lot of misinformation about the negative binomial model, and especially with respect to the dispersion statistic and dispersion parameter.

The dispersion statistic, which gives an indication of count model extra-dispersion, is the Pearson statistic divided by the residual DOF. $\mu$ is the location or shape parameter. For count models, the scale parameter is set at 1. The R glm and glm.nb $\theta$ is a dispersion parameter, or ancillary parameter. I called it the heterogeneity parameter in the first edition of my book, Negative Binomial Regression (2007, Cambridge University Press), but call it the dispersion parameter in my 2011 second edition. I give a complete rationale for the various terms in the NB model in my forthcoming book, Modeling Count Data (Cambridge) which is going to press today. It should be for sale (paperback) by July 15.

glm.nb and glm are unusual in how they define the dispersion parameter. The variance is given as $\mu+\frac{\mu^2}{\theta}$ rather than $\mu+\alpha\mu^2$, which is the direct parameterization. It is the way NB is modeled in SAS, Stata, Limdep, SPSS, Matlab, Genstat, Xplore, and most all software. When you compare glm.nb results with other software results, remember this. The author of glm (which came from S-plus) and glm.nb apparently took the indirect relationship from McCullagh & Nelder, but Nelder (who was the co-founder of GLM in 1972) wrote his kk system add-on to Genstat in 1993 in which he argued that the direct relationship is preferred. He and his wife used to visit me and my family about every other year in Arizona starting in early 1993 until the year before he died. We discussed this pretty thoroughly, since I had put a direct relationship into the glm program I wrote in late 1992 for Stata and Xplore software, and for a SAS macro in 1994.

The nbinomial function in the msme package on CRAN allows the user to employ the direct (default) or indirect (as an option, to duplicate glm.nb) parameterization, and provides the Pearson statistic and residuals to output. Output also displays the dispersion statistic, and allow the user to parameterize $\alpha$ (or $\theta$), giving parameter estimates for the dispersion. This allows you to assess which predictors add to the extra-dispersion of the model. This type of model is generally referred to as heterogeneous negative binomial. I'll put the nbinomial function into the COUNT package before the new book comes out, plus a number of new functions and scripts for graphics.

Answer 3

glm reference negative binomial :

Wikipedia negative binomial 'r' is glm's 'theta' which implies glm 'theta' is shape parameter. In Simple terms, glm's 'theta' is number of failures.

Answer 4

Here comes an attempt to clean up the "Negative Binomial Parametrisation Confusion". In the following, $\mathrm{NB}(x; \ldots)$ stands for the PMF of the Negative Binomial (N.B.).

1. The $(r,p)$ parametrisation: $\mathrm{NB}(x; r,p)$ is the probability of observing $x$ "failures" in a series of independent Bernoulli trials with "success" probability $p$ until $r$ "successes" are observed. (Note that because "success" and "failure" are arbitrary labels of the two outcomes of a Bernoulli trial, you may see definitions talking about observing $r$ "failures" and $p$ being the probability of "failure".) In this parametrisation the PMF is: $$\mathrm{NB}(x; r,p) = \binom{r+x-1}{r-1} p^r (1-p)^x$$ its mean is $$\mu = r\frac{1-p}{p}$$ and the variance is $$\sigma^2 = r\frac{(1-p)}{p^2} = \frac{\mu}{p}$$ If $p=1$ then $\sigma^2 = \mu$ and the N.B. becomes a Poisson. Otherwise, when $0 < p < 1$, then $\sigma^2 > \mu$, the N.B. is overdispersed.
2. The $(\mu,\theta)$ parametrisation: To emphasise the relationship between the N.B. and the Poisson, let's parametrise the N.B. with its mean $\mu$. From the formula for the N.B. mean we can see that $$r = \frac{p}{1-p} \mu$$ where $\frac{p}{1-p}$ is the odds of "success". Although in the original definition of the N.B. $r$ was an integer ("number of successes"), from this formula it's clear that it could be a non-negative real number. To emphasise this insight, let's rename $r \rightarrow \theta$. From the equation above we can work out the value of $p$ in the $(\mu, \theta)$ parametrisation: $$p = \frac{\theta}{\theta + \mu}$$ After replacing the factorials in the PMF with Gamma functions we get: $$\mathrm{NB}(x; \mu, \theta) = \frac{\Gamma(\theta + x)}{\Gamma(\theta) x!} \left(\frac{\theta}{\theta + \mu}\right)^\theta \left(\frac{\mu}{\theta + \mu}\right)^x$$ The variance in this formulation is: $$\sigma^2 = \mu + \frac{\mu^2}{\theta}$$

We can see that if $\theta \rightarrow \infty$ then $\sigma^2 \rightarrow \mu$, the N.B. becomes the Poisson.

To call $\theta$ the "dispersion parameter" is somewhat confusing because the smaller $\theta$ gets, the larger is the overdispersion. Further confusion arises from the fact that $r$ and $\theta$ are essentially the same thing. This is why in the R functions dnbinom, pnbinom etc. both are set with the size= parameter, just to make sure you can't recognise them :-). Moreover, I think size= is a misnomer as $\theta$ is actually a kind of "shape" parameter.

It would have been more intuitive to define an overdispersion parameter $\alpha = 1/r$, which would give $$\sigma^2 = \mu + \alpha \mu^2$$ as explained by Joseph Hilbe above. In this formulation $\alpha=0$ indicates no overdispersion. But it's probably too late for R to change, I guess.