Appendix C: 

Negative Binomial Regression Models and 

Estimation Methods 

Dominique Lord 

Zachry Dept. of 

Civil Engineering 

Texas A & M University 

College Station, TX 

Byung-Jung Park 

Korea Transport Institute 

Goyang, South Korea 

This appendix presents the characteristics of Negative Binomial regression models and 

discusses their estimating methods. 

Probability Density and Likelihood Functions 

The properties of the negative binomial models with and without spatial intersection are 

described in the next two sections. 

Poisson-Gamma Model 

The Poisson-Gamma model has properties that are very similar to the Poisson model 

discussed in Appendix B, in which the dependent variable yi is modeled as a Poisson variable 

with a mean i where the model error is assumed to follow a Gamma distribution. As its name 

implies, the Poisson-Gamma is a mixture of two distributions and was first derived by 

Greenwood and Yule (1920). This mixture distribution was developed to account for overdispersion that is commonly observed in discrete or count data (Lord et al., 2005). It became 

very popular because the conjugate distribution (same family of functions) has a closed form and 

leads to the negative binomial distribution. As discussed by Cook (2009), “the name of this 

distribution comes from applying the binomial theorem with a negative exponent.” There are two 

major parameterizations that have been proposed and they are known as the NB1 and NB2, the 

latter one being the most commonly known and utilized. NB2 is therefore described first. Other 

parameterizations exist, but are not discussed here (see Maher and Summersgill, 1996; Hilbe, 

2007). 

C.1
 

 

 

 

NB2 Model 

Suppose that we have a series of random counts that follows the Poisson distribution: 

(C.1)

i 1, 2,  n ; and 

distribution. If the Poisson mean is assumed to have a random intercept term and this term enters 

the conditional mean function in a multiplicative manner, we get the following relationship 

(Cameron and Trivedi, 1998): 

where yi is the observed number of counts for    i is the mean of the Poisson 

(C.2) 

where is defined as a random intercept; is the log-link

i

between the Poisson mean and the covariates or independent variables xs ; and the s are the 

regression coefficients. As discussed in Appendix B, the relationship can also be formulated 

using vectors, such that

The marginal distribution of yican be obtained by integrating the error term, , 

(C.3)

where h i is a mixing distribution. In the case of the Poisson-Gamma mixture, is 

 is the Gamma distribution. This distribution has a closed 

form and leads to the NB distribution. 

C.2

the Poisson distribution and
 

 

 

 

Assume that the variable  i follows a two-parameter Gamma distribution: 

(C.4)

gives us the one-parameter Gamma

where E  i  1 and VAR  i   1



. We can transform the Gamma distribution as a function of

the Poisson mean, which gives the following probability density function (PDF; Cameron and 

Trivedi, 1998): 

(C.5)

Combining equations C-1 and C-5 with equation C-3 yields the marginal distribution of 

yi : 

(C.6)

Using the properties of the Gamma function, it can be shown that equation C-6 can be 

expressed as: 

( C-7)

The PDF of the NB2 model is therefore the last part of Equation C-7: 

(C.8)

C.3

where, Setting
 

 

 

 

 (C.13)

 

 (C.14)

 

Note that the PDF has also been defined in the literature as: 

(C.9)

The first two moments of the NB2 are the following: 

(C.10)

(C.11)

The next step consists of defining the log-likelihood function of the NB2. It can be 

shown that: 

(C.12)

By substituting equation C-12 into C-8, the log-likelihood can be computed using the 

following equation: 

Note also that the log-likelihood has also been expressed as: 

In the statistical literature, the Poisson-Gamma model has also been defined as: 

(C.15)

where the mean of the Poisson is structured as: 

C.4

Recall that
 

 

 

 

 

 (C.19) 

 

 

 

(C.16)

and where, f (.) is a function of the covariates, X (Miaou and Lord, 2003). As before, β is a 

vector of coefficients and  i is the model error independent of all the covariates with mean equal 

to 1 and a variance equal to1  . 

NB1 Model 

The NB1 is very similar to the NB2, but the parameterization of the variance (the second 

moment) is slightly different than in equation C-11. 

(C.17) 

(C.18)

The log-likelihood of the NB1 is given by: 

The NB1 is usually less flexible in capturing the variance and is not used very often by 

analysts and statisticians. Interested readers are referred to Cameron and Trivedi (1998) for 

additional information about this parameterization. 

Poisson-Gamma Model with Spatial Interaction 

The Poisson-Gamma (or negative binomial model) can also incorporate data that are 

collected spatially. To capture this kind of data, a spatial autocorrelation term needs to be added 

to the model. Using the notation described in Equation C-15, the NB2 model with spatial 

interaction can be defined as: 

(C.20) 

(C.21)

C.5

with the mean of Poisson-Gamma organized as:
 

 

 

The assumption on the uncorrelated error term  i is the same as in the Poisson-Gamma 

model described above; same as before, namely i  exp(x ' 

i

β) . The third term in the expression, 

i , is a spatial random effect, one for each observation. Together, the spatial effects are 

distributed as a complex multivariate normal (or Gausian) density function. In other words, the 

second model is a spatial regression model within a negative binomial model. 

There are two common ways to express the spatial component, either as a Conditional 

Autoregressive (CAR) or as a Simultaneous Autoregressive (SAR) function (De Smith et al., 

2007). The CAR model is expressed as: 

(C.22) 

where μi is the expected value for observation i, wij is a spatial weight between the observation, i, 

and all other observations, j (and for which all weights sum to 1.0), and ρ is a spatial 

autocorrelation parameter that determines the size and nature of the spatial neighborhood effect. 

Note that there are different weight factors that have been proposed, such as the inverse distance 

weight function, exponential distance decay weight function and the Gaussian weighting 

function among others. The summation of the spatial weights times the difference between the 

observed and predicted values is over all other observations (i≠j). The reader is referred to 

Haining (1990) and LeSage (2001) for further details. 

The SAR model has a simpler form and can be expressed as: 

(C.23)

where the terms are as defined above. Note that in the CAR model the spatial weights are 

applied to the difference between the observed and expected values at all other locations whereas 

in the SAR model, the weights are applied directly to the observed value. In practice, the CAR 

and SAR models produce very similar results. Additional information about the PoissonGamma-CAR is described below. 

Estimation Methods 

This section describes two methods that can be used for estimating the coefficients of the 

regression NB models. The two methods are the maximum likelihood estimates (MLE) and the Monte 

Carlo Markov Chain (MCMC). 

C.6
 

 

 

 

 

Maximum Likelihood Estimation 

The characteristics of the MLE method were described in Appendix B for the normal and Poisson 

regression. The same characteristics apply here. The coefficients of the NB regression model are 

estimated by taking the first-order conditions and making them equal to zero. There are two first-order 

equations, one for the model’s coefficients and one for the dispersion parameter (Lawson, 1987). The two 

for the NB2 are as follows: 

 (C.24a) 

(C.24b)

where xi is a vector of covariates. 

Similar to the Poisson model, the series of equations can be solved using the Newton-Raphson 

procedure or the scoring algorithm. The confidence intervals on the s and1 can be calculated 

using the covariance matrix that is assumed to be normally distributed: 

(C.25)

where, 

(C.26a)

(C.26b)

It should be pointed out that the NB2 with spatial interaction model (Poisson-Gamma-CAR) 

cannot be estimated using the MLE method. It needs to be estimated using the MCMC technique, which 

is described next. 

C.7
 

 

 

 

 

parameter  (or the dispersion parameter γ=1/ ). Ideally, samples need to be drawn of each parameter 

from the joint posterior distribution. However, the form in Equation C-28b i s very complex and it is 

difficult to draw samples from such a distribution. Consequently, samples are drawn from the full 

Monte Carlo Markov Chain Estimation 

This section discusses how to draw samples from the posterior distribution of the PoissonGamma model and Poisson-Gamma-Conditional Autoregressive (CAR) model using the MCMC 

technique. 

MCMC Poisson-Gamma Model 

The Poisson-Gamma model can be formulated from a two-stage hierarchical Poisson model: 

(C.27a)

(C.27b) 

where  ( ) is the prior distribution imposed on the Poisson mean, i with a prior parameter  , and ,

  () is the hyper-prior on  with known hyper-parameters (a, b, for example). 

In Equations C-27a and C-27b, if we specify i  i i (where  i ( e i ) ~ Gamma ( , ) in the 

first stage and  ~ Gamma (a, b) in the second stage), these result in exactly the NB2 regression model 

described in the previous section. With this specification, it is also easy to show that i in the first stage 

follows Gamma (, / i ) as shown in Equation C-5. Note that i  exp(x ' 

i

β) as described above. 

For simplicity, if a flat uniform prior is assumed for each  j ( j  0, 1, , J ) and the 

parameters  s and  are mutually independent, the joint posterior distribution for the Poisson-Gamma 

model is defined as: 

(C.28a) 

(C.28b)

The parameters of interest are and the inverse dispersion 

C.8
 

 

 

 

 

 

 

 

 

conditional distribution sequentially (that is, one at a time). This iterative process is called the Gibbs 

sampling method. 

Therefore, once the full conditionals are known for each parameter, Gibbs sampling can be 

implemented by drawing samples of each parameter sequentially. The full conditional distributions for 

each parameter for the Poisson-Gamma model can be easily derived from Equation C-28b and are given 

as (Park, 2010): 

However, unlike Equation C-29a, the full conditional distributions for the  s and  (Equations 

C-29b and C-29c) do not belong to any standard distribution family so it is not easy to draw samples

directly from their full conditional distributions. While there are several approaches to sampling from

such a complex distribution, the particular sampling algorithm used in CrimeStat is a Metropolis-Hastings

(or MH) algorithm with slice sampling of individual parameters.

The MCMC sampling procedure using the slice sampling algorithm within Gibbs sampling, 

therefore, can be summarized as follows: 

λ(0) β(0) (0) 1. Start with initial values , and . Repeat the following steps for

t  1, , T0 , , T0 T .

β(t-1) (t-1) λ(t ) 2. Step 1: Conditional on knowing and , draw from Equation C-29a

independently for i 1, 2, , n .

λ(t ) (t-1) β(t ) 3. Step 2: Conditional on knowing and , draw from Equation C-29b

independently for j  0, 1, , J using the slice sampling method.

λ(t ) β(t ) (t ) 4. Step 3: Conditional on knowing and , draw  from Equation C-29c using the

slice sampling method.

C.9

(C.29a) 

(C.29b)

(C.29c)
 

 

 

5. Step 4: Store the values of all parameters Increase t by one an d

return to Step 1.

6. Step 5: Discard the first k draws as a burn-in period, where k is defined by the user.

After equilibrium is reached at the k

th iteration, sampled values are averaged to provide the 

consistent estimates of the parameters: 

(C.30)

where denotes any parameter of interest in the model. 

MCMC Poisson-Gamma-CAR Model 

way

(C.31)

where i is a spatial random effect, one for each observation. As in the Poisson-Gamma model, we 

specify e

 i ~ Gamma (, ) to model the independent error term. To model the spatial effect, i , we 

assume the following: 

(C.32)

where is the probability of a spatial effect given a lagged spatial effect, which

sums all over all records, j (i.e., all other zones) except for the record of interest, i . This formulation 

ij gives a conditional normal density with mean =   j and variance = . The parameter 

determines the direction and overall magnitude of the spatial effects. The term wij is a spatial weight 

function between zones i and j . In the algorithm, the term for the variance is 2 1/ and the same  

variance is used for all observations. 

C.10

For the Poisson-Gamma-CAR model, the only difference from the Poisson-Gamma model is the 

is structured. The mean of Poisson-Gamma-CAR is organized as:
 

 

 

We define the spatial weight matrix W with the entries wij and the diagonal entries The 

matrix D is defined as a diagonal matrix with the diagonal entries, Sun, Tsutakawa, and Speckman 

(1999) show that if are the smallest and largest eigenvalues of 

respectively, then has Φ a multivariate normal distribution with mean 0 and nonsingular

covariance matrix

(C.33)

Prior Distributions for MCMC Poisson-Gamma-CAR 

For the prior distributions, we assume the following distributions for each parameter: 

Parameter Prior distribution

 j ( j  0, 1, ,  J ) Uniform (, ) 

 Gamma (a , b )  

2  (  )   Gamma (a , b )  

 1 1 Uniform ( ,  ) min max 

The parameters in the Poisson-Gamma-CAR model are

Then, the random samples can be drawn from the full conditional

distributions of each parameter. It can be shown that the full conditional distributions for each parameter 

are given as follows: 

(C.34a) 

(C.34b)

(C.34c) 

C.11
 

 

where are the eigenvalues of 

Since the full conditional distributions were specified, the Gibbs sampling method can be applied 

sequentially. It is easy to generate random samples from Equations C-34a and C-34e. The other full 

conditional distributions are not of closed form, so the slice sampling method should be applied. 

Likelihood Statistics 

There are many measures that can be used for estimating how well the model fits the data. Some 

of them have already been discussed in Appendix B but are also included here for the sake of 

completeness. They fall into three groups. First, there are statistics for indicating the likelihood level of a 

model, that is, how well the model maximizes the likelihood function. Among these statistics are: 

Akaike Information Criterion (AIC) 

The AIC is another measure of fit that can be used to assess models. This measure also uses the 

log-likelihood, but add a penalizing term associated with the number of variables. It is well known that by 

adding variables, one can improve the fit of models. Thus, the AIC tries to balance the goodness-of-fit 

versus the inclusion of variables in the model. The AIC is computed as: 

AIC  2ln L  2 p (C.37)

where p is the number of unknown parameters included in the model (this also includes the inverse 

dispersion parameter  and random spatial effect f i ) and ln L is the log-likelihood described in 

Equations C-13 or C-14. Smaller values are better. 

Bayes Information Criterion (BIC) 

Similar to the AIC, the BIC also employs a penalty term associated with the number of 

parameters ( p ) and the sample size ( n ). This measure is also known as the Schwarz Information 

Criterion. It is computed the following way: 

C.12

(C. )34e

(C.34f)

(C.34d)
 (C.38)

 

 

 

 

 

 Again, smaller values are better. 

Deviance Information Criterion (DIC) 

When the Bayesian estimation method is used, the DIC is often used as a goodness-of-fit 

(GOF) measure instead of the AIC or BIC. The latter ones are generally used for the maximum 

likelihood method. The DIC is defined as follows: 

(C.39)

where D is the average of the deviance (  2 ln L ) over the posterior distribution, and Dˆ is the 

deviance calculated at the posterior mean parameters. As with the AIC and BIC, the DIC uses 

pD  D  Dˆ (effective number of parameters) as a penalty term on the goodness of fit. Differences in

DIC from 5-10 indicate that one model is clearly better (Spiegelhalter et al., 2002). 

Deviance 

The deviance is a measure of goodness of fit that can be used to assess models. It is defined as 

twice the difference between the maximum likelihood achievable ( yi  ˆi ) and the likelihood of the

fitted model (the ^ refers to the estimate of the variable that is based on the data): 

(C.35)

For the NB2 model, the deviance can be computed as: 

(C.36)

Smaller values mean that the model fits the data better. 

Pearson Chi-Square 

Another useful likelihood statistic is the Pearson Chi-square and is defined as 

(C.37) 

C.13
 

 

 

 

If the mean and the variance are properly specified, then

(Cameron and Trivedi, 1998). Values closer to n (the sample size) show a better fit. Recall that the 

variance for the NB2 model is 

Model Error Estimates 

Second, there are statistics for estimating how well the model fit the data and the converse, how 

much error was in the model. Two error statistics are particularly useful. 

Mean Absolute Deviation (MAD) 

This criterion has been proposed by Oh et al. (2003) to evaluate the fit of models. The Mean 

Absolute Deviance (MAD) calculates the absolute difference between the estimated and observed values 

(C.38)

Mean Squared Prediction Error (MSPE) 

The Mean Squared Prediction Error (MSPE) is a traditional indicator of error and calculates the 

difference between the estimated and observed values squared. 

(C.39)

A value closer to 1 means the model fits the data better. 

Over-dispersion Tests 

Third, there are statistics for indicating the degree of over-dispersion in the model, including: 

Adjusted Deviance 

The adjusted deviance is defined as the deviance divided by the degrees of freedom (N-K-1). A 

value closer to 1 indicates a satisfactory GOF. Usually, values greater than 1 indicate signs of overdispersion, while values below 1 show signs of under-dispersion. 

Adjusted Pearson Chi-Square 

The adjusted Pearson Chi-square is defined as the Pearson Chi-square divided by the degrees of 

freedom. A value closer to 1 indicates a satisfactory goodness-of-fit. 

C.14
 

 

Dispersion Multiplier 

The dispersion multiplier measures the extent to which the conditional variance exceeds the 

conditional mean (conditional on the independent variables and the intercept term) and is defined by 

Inverse Dispersion Multiplier 

The inverse dispersion multiplier ( ) is simply the reciprocal of the dispersion multiplier 

( 1/ ) ; some users are more familiar with it in this form. 

It should be pointed out that many GOF measures are not useful when a single model is 

evaluated. The measures are therefore relevant when several models are compared with each other (i.e., 

different functional forms or when different variables are included in the models). 

There are other measures that can be used for estimating the goodness-of-fit and the amount of 

error in models, but are not presented here. Readers can find additional measures in Mitra and 

Washington (2007) and Lord and Park (2008). 

C.15
References 

Cameron, A.C. & Trivedi, P.K. (1998). Regression Analysis of Count Data. Cambridge 

University Press. Cambridge, UK. 

Cook, J.D. (2009). Notes on the Negative Binomial Distribution. 

www.johndcook.com/negative_binomial.pdf (accessed on Aug. 15th, 2010). 

De Smith, M., Goodchild, M.F., & Longley, P.A. (2007). Geospatial Analysis (second edition). 

Matador: Leicester, U.K. 

Greenwood, M., & Yule, G.U. (1920). An Inquiry into the Nature of Frequency Distributions of 

Multiple Happenings, with Particular Reference to the Occurrence of Multiple Attacks of 

Disease or Repeated Accidents. Journal of the Royal Statistical Society A, Vol. 83, pp. 255-279. 

Haining, R. (1990). Spatial Data Analysis in the Social and Environmental Sciences. Cambridge 

University Press, Cambridge, U.K. 

Hilbe, J.M. (2007) Negative Binomial Regression. Cambridge University Press, Cambridge, U.K. 

Lawless, J.F. (1987). Negative binomial and mixed Poisson regression. The Canadian Journal of 

Statistics, Vol. 15, No. 3, pp. 209-225. 

LeSage, J.P. (2001). A Family of Geographically Weighted Regression Models. Working Paper. 

Department of Economics, University of Toledo, Toledo, OH. 

Lord, D., & Park, P.Y.J. (2008). Investigating the effects of the fixed and varying dispersion 

parameters of Poisson-Gamma models on empirical Bayes estimates. Accident Analysis & 

Prevention, Vol. 40, No. 4, pp. 1441-1457. 

C.16
References (continued)

Lord, D., Washington, S.P., & Ivan, J.N. (2005). Poisson, Poisson-Gamma and Zero Inflated 

Regression Models of Motor Vehicle Crashes: Balancing Statistical Fit and Theory. Accident 

Analysis & Prevention. Vol. 37, No. 1, pp. 35-46. 

Maher M.J., & Summersgill, I. (1996) A Comprehensive Methodology for the Fitting Predictive 

Accident Models. Accident Analysis & Prevention, Vol. 28, No. 3, pp.281-296. 

Miaou, S.-P., & Lord, D. (2003). Modeling TraffiC.Flow Relationships at Signalized 

Intersections: Dispersion Parameter, Functional Form and Bayes vs Empirical Bayes. 

Transportation Research Record 1840, pp. 31-40. 

Mitra, S., & Washington, S. (2007). On the nature of over-dispersion in motor vehicle crash 

prediction models. Accident Analysis & Prevention, Vol. 39, No. 3, pp. 459-468. 

Oh, J., Lyon, C., Washington, S., Persaud, B., & Bared, J. (2003). Validation of FHWA crash 

models for rural intersections: lessons learned. Transportation Research Record 1840, pp. 41-49. 

Park, B.-J. (2010). Application of the Finite Mixture Models for Vehicle Crash Data Analysis. 

PhD Dissertation, Texas A&M University, College Station, TX. 

Spiegelhalter, D.J., Best, N.G., Carlin, B.P., van der Linde, A. (2002) Bayesian measures of 

model complexity and fit. Journal of the Royal Statistical Society, Series B, Vol. 64, No. 4, pp. 

583-639.

Sun, D., Tsutakawa, R. K., & Speckman, P. L. (1999). Posterior distribution of hierarchical 

models using CAR(1) distributions. Biometrika, Vol. 86, No. 2, pp. 341-350. 

C.17