In the R function, rnbinom, one of the parameters is the dispersion or shape parameter. This can be parameterized as theta or alpha, depending on how the model is written. I can't tell from ?rnbinom what its asking for. Anyone have an idea?

EDIT: I've run a simple negative binomial regression model, and want to use the model parameters to produce the theoretical distribution for simulation work. I'm not exactly sure how to use the dispersion parameter. Here's the output from R:


glm.nb(formula = exit ~ 1 + offset(log(stock)), data = dt, init.theta = 5.855047422, 
    link = log)

Deviance Residuals: 
     Min        1Q    Median        3Q       Max  
-2.83778  -0.86369   0.00863   0.62604   1.80784  

            Estimate Std. Error z value Pr(>|z|)    
(Intercept)   -3.689      0.029  -127.2   <2e-16 ***
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

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

    Null deviance: 218.61  on 211  degrees of freedom
Residual deviance: 218.61  on 211  degrees of freedom
AIC: 2297.5

Number of Fisher Scoring iterations: 1

              Theta:  5.855 
          Std. Err.:  0.582 

 2 x log-likelihood:  -2293.500 

I will use rnbinom to model the distribution, taking as parameters:


My question is if I'm parametrizing the size parameter appropriately. Should it be 5.855, or 1/5.855? I more or less understand the different parametrizations of the model, as either $\theta$ (or $r$) or $\alpha$, and from here I know glm.nb is reporting $\theta$. I'm not exactly sure what rnbinom is looking for with its size parameter - am I correct in assuming it is $\theta$, and my code here correct (size=5.855).


The documentation calls it "size":

size   target for number of successful trials, or dispersion parameter (the shape parameter of the gamma mixing distribution). Must be strictly positive, need not be integer.

That is the simplest way to understand it. The negative binomial distribution is typically understood as:

a discrete probability distribution of the number of successes in a sequence of independent and identically distributed Bernoulli trials before a specified (non-random) number of failures (denoted r) occurs.

In other words, when performing a series of coin flips, you could count how many tails you got before you got $r$ heads, with a coin that has a $p$ probability of heads.

  • 1
    $\begingroup$ also useful, possibly worth linking in your answer: stats.stackexchange.com/questions/6728/… ; my answer there explains one interpretation of non-integer shape parameters ... this question might be a dupe, it depends how narrowly you interpret it ... $\endgroup$ – Ben Bolker Nov 3 '16 at 0:02
  • $\begingroup$ To clarify my Q, I've got a regression model using glm.nb, and want to use the parameters from the model to draw from the negative binomial distribution. glm.nb produces theta - my understanding is the reciprocal of this is alpha (e.g. that's what Stata's nbreg outputs). I'm wondering if I want to use rnbinom to simulate what the theoretical distribution would be based on the empirical parameters, am I using size=theta or size=1/theta? From the wikipedia page it looks like they are saying that r=1/alpha, which would be r=theta. $\endgroup$ – robin.datadrivers Nov 7 '16 at 21:56
  • $\begingroup$ That's a long way from what your wrote in the question, @robin.datadrivers (& I'm not sure if it's an on topic statistical question or an off topic R coding question). Why don't you edit your question to state what you're working with & what you want to do. Eg, you could paste in the output from your NB regression model & show the rnbinom() call you are trying. $\endgroup$ – gung - Reinstate Monica Nov 7 '16 at 22:04

I have recently just had trouble with this question, and as I found no concrete answer on the edit of the original post I will share my findings.

The glm.nb function estimates a dispersion parameter, labelled $\theta$. This is often called the size parameter and notated by $r$, as shown in the Negative Binomial Wikipedia article.

The rnbinom function (with documentation) takes prob and size as parameters, or alternatively mu and size, where mu is calculated by prob = size/(size+mu).

Therefore the post is correct in its assumption that you want:

x <- rnbinom(nrow(dt), size=5.855, mu=1/exp(-3.689))

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.