Negative Binomial Regression

Michael L. Zwilling

Negative binomial regression is implemented using maximum likelihood estimation. The traditional model and the rate model with offset are demonstrated, along with regression diagnostics.

Traditional Model

Negative binomial regression is a type of generalized linear model in which the dependent variable Y is a count of the number of times an event occurs. A convenient parametrization of the negative binomial distribution is given by Hilbe []:

p(y)=P(Y=y)=((Γ(y+1/α))/(Γ(y+1)Γ(1/α)))((((1/(1+α μ))))^(1/α))(((((α μ)/(1+α μ))))^y) ,

where μ>0 is the mean of Y and α>0 is the heterogeneity parameter. Hilbe [] derives this parametrization as a Poisson-gamma mixture, or alternatively as the number of failures before the (((1/α))^th) success, though we will not require 1/α to be an integer.

The traditional negative binomial regression model, designated the NB2 model in [], is

ln μ=(Subscript[β, 0])+(Subscript[β, 1])(Subscript[x, 1])+(Subscript[β, 2])(Subscript[x, 2])+⋯+(Subscript[β, p])(Subscript[x, p]),

where the predictor variables (Subscript[x, 1]),(Subscript[x, 2]),…,(Subscript[x, p]) are given, and the population regression coefficients (Subscript[β, 0]),(Subscript[β, 1]),(Subscript[β, 2]),…,(Subscript[β, p]) are to be estimated.

Given a random sample of n subjects, we observe for subject i the dependent variable (Subscript[y, i]) and the predictor variables (Subscript[x, 1,i]), (Subscript[x, 2,i]),…,(Subscript[x, p,i]). Utilizing vector and matrix notation, we let β=(({{(Subscript[β, 0]), (Subscript[β, 1]), (Subscript[β, 2]), ⋯, (Subscript[β, p])}}))ᵀ, and we gather the predictor data into the design matrix X as follows:

X=(({{1, (Subscript[x, 1,1,]), (Subscript[x, 1,2]), ⋯, (Subscript[x, 1,p])}, {1, (Subscript[x, 2,1]), (Subscript[x, 2,2]), ⋯, (Subscript[x, ,2,p])}, {⋮, ⋮, ⋮, ,, ⋮}, {1, (Subscript[x, n,1]), (Subscript[x, n,2]), ⋯, (Subscript[x, n,p])}})).

Designating the (i^th) row of X to be (Subscript[x, i]), and exponentiating (), we can then write the distribution () as

p((Subscript[y, i]))=((Γ((Subscript[y, i])+1/α))/(Γ((Subscript[y, i])+1)Γ(1/α)))((((1/(1+α (e^((Subscript[x, i])·β))))))^(1/α))(((((α (e^((Subscript[x, i]) ·β)))/(1+α (e^((Subscript[x, i]) ·β))))))^(Subscript[y, i])), i=1, 2, …, n.

We estimate α and β using maximum likelihood estimation. The likelihood function is

L(α,β)=(∏_(i=1)^n)p((Subscript[y, i]))=(∏_(i=1)^n)((Γ((Subscript[y, i])+1/α))/(Γ((Subscript[y, i])+1)Γ(1/α)))((((1/(1+α (e^((Subscript[x, i])·β))))))^(1/α))(((((α (e^((Subscript[x, i])·β)))/(1+α (e^((Subscript[x, i])·β))))))^(Subscript[y, i])),

and the log-likelihood function is

ln L(α,β)=(∑_(i=1)^n)((Subscript[y, i])ln α+(Subscript[y, i]) ((Subscript[x, i])·β)-((Subscript[y, i])+(1/α))ln(1+α (e^((Subscript[x, i]) ·β)))+ln Γ((Subscript[y, i])+(1/α))-ln Γ((Subscript[y, i])+1)-ln Γ((1/α))).

The values of α and β that maximize ln L(α,β) will be the maximum likelihood estimates we seek, and the estimated variance-covariance matrix of the estimators is Σ=-(H^(-1)), where H is the Hessian matrix of second derivatives of the log-likelihood function. Then the variance-covariance matrix can be used to find the usual Wald confidence intervals and p-values of the coefficient estimates.

Example 1: Traditional Model with Simulated Data

We will use Mathematica to replicate some examples given by Hilbe [], who uses R and Stata. We start with simulated data generated with known regression coefficients, then recover the coefficients using maximum likelihood estimation. We will generate a sample of n=5000 observations of a dependent random variable Y that has a negative binomial distribution with mean given by (), using p=2, α=0.50, and β=(2.00 ).75 -1.25)ᵀ. The design matrix (Subscript[X, 1]) will contain independent standard normal variates.

SeedRandom[12345];
n1=5000;
X1=Prepend[#,1]&/@RandomVariate[NormalDistribution[],{n1,2}];
y1=Table[RandomVariate[NegativeBinomialDistribution[(1/α),(1/(1+α (E^(X1[[i]].β))))]/.{α->0.50,
β->{2.00,0.75,-1.25}}],{i,n1}];

Now we define and maximize the log-likelihood function (), obtaining the estimates of α and β. Some experimentation with starting values for the search may be required, and the accuracy goal may need to be lowered; we could obtain good starting values for β using Poisson regression via GeneralizedLinearModelFit, while α is usually between 0.0 and 4.0 [].

GeneralizedLinearModelFit[{X1,y1},ExponentialFamily->"Poisson"]

But we arbitrarily set all starting values to 1.0 and successfully find the correct estimates.

lnL1[X_,y_,α_,β_]:=Module[{n=Length[y]},(∑_(i=1)^n)(y[[i]] Log[α]+y[[i]] (X[[i]].β)-(y[[i]]+(1/α))Log[1+α (E^(X[[i]].β))]+LogGamma[y[[i]]+(1/α)]-LogGamma[y[[i]]+1]-LogGamma[(1/α)])]

results1=FindMaximum[lnL1[X1,y1,α,β],{{α,1.0,0.01,4.0},{β,{1.0,1.0,1.0}}},AccuracyGoal->6]

Define two helper functions.

βs[X_]:=Array[β,Length[X[[2]]],0]

flattenrules[{_,{xx_->xxval_,yy_->yyvals_}}]:={xx->xxval,Thread[Array[yy,Length[yyvals],0]->yyvals]}//Flatten

flattenrules[results1]

Next, we find the standard errors of the estimates. The standard errors are the square roots of the diagonal elements of the variance-covariance matrix Σ, which as mentioned above is given by Σ=-(H^(-1)), where H is the Hessian matrix of second derivatives of the log-likelihood function. First, define the Hessian for any function.

HessianH[f_,x_]:=D[f,{Flatten[x],2}]

Then we find the Hessian and Σ at the values of our parameter estimates.

VarianceCovarianceMatrix[X_,y_,results_]:=-Inverse[HessianH[lnL1[X,y,α,β]/.β->βs[X],{α,βs[X]}]/.flattenrules[results]]

(v1=VarianceCovarianceMatrix[X1,y1,results1])//MatrixForm

Finally, these are our standard errors.

StdErr1=(Sqrt[Diagonal[v1]])

We can now print a table of the results: the estimates of the coefficients, their standard errors, and the Wald z statistics, p-values, and confidence intervals.

coefficients1=Last/@flattenrules[results1]

CoefficientsTable[coefficients_,StdErr_,labels_]:=Module[{zStatistics,pValues,displaypValues,lower,upper},
zStatistics=coefficients/StdErr;
pValues=2CDF[NormalDistribution[],-Abs[zStatistics]];
displaypValues=Table[{NumberForm[Round[pValues[[j]],0.0001],{5,4}]},{j,Length[labels]}]//Flatten;
lower=coefficients-1.96 StdErr;
upper=coefficients+1.96 StdErr;
Text@Column[{TableForm[{coefficients,StdErr,zStatistics,displaypValues}//Transpose,TableHeadings->{labels,{"Estimate","Std. Err.",Row[{Style["z",Italic],"\[Hyphen]Statistic"}],TraditionalForm[P>Abs[z]]}}],
"",
TableForm[{coefficients,lower,upper}//Transpose,TableHeadings->{labels,{"Estimate","95% CI lower","95% CI upper"}}]
}]
]

CoefficientsTable[coefficients1,StdErr1,Prepend[βs[X1],α]]

We see that in each case the confidence interval has captured the population parameter.

Traditional Model for Rates, Using Offset

If the dependent variable Y counts the number of events during a specified time interval t, then the observed rate Y/t can be modeled by using the traditional negative binomial model above, with a slight adjustment. We note that t can also be thought of as area or subpopulation size, among other interpretations that lead to considering Y/t a rate.

Since E(Y/t)=μ/t, we make the following adjustment to model () above:

ln(μ/t)=(Subscript[β, 0])+(Subscript[β, 1])(Subscript[x, 1])+(Subscript[β, 2])(Subscript[x, 2])+⋯+(Subscript[β, p])(Subscript[x, p]),

which can also be written as:

ln μ=(Subscript[β, 0])+(Subscript[β, 1])(Subscript[x, 1])+(Subscript[β, 2])(Subscript[x, 2])+⋯+(Subscript[β, p])(Subscript[x, p])+ln t.

This last term, ln t, is called the offset. So in our log-likelihood function, instead of replacing μ with (e^((Subscript[x, i])·β)), we replace μ with (e^((Subscript[x, i])·β+ln t)), resulting in the following:

ln L(α,β)=(∑_(i=1)^n)((Subscript[y, i])ln α+(Subscript[y, i]) ((Subscript[x, i])·β+ln (Subscript[t, i]))-((Subscript[y, i])+(1/α))ln(1+α (e^((Subscript[x, i])·β+ln (Subscript[t, i]))))+ln Γ((Subscript[y, i])+(1/α))-ln Γ((Subscript[y, i])+1)-ln Γ((1/α))).

Then we proceed as before, maximizing the new log-likelihood function in order to estimate the parameters.

Example 2: Traditional Model with Offset for the Titanic Data

The Titanic survival data, available from [] and analyzed in [] using R and Stata, is summarized in Table , with crew members deleted.

Why did fewer first-class children survive than second class or third class? Was it because first-class children were at extra risk? No, it was because there were fewer first-class children on board the Titanic in the first place. So we do not want to model the raw number (Y) of survivors; instead, we want to model the proportion (Y/cases) of survivors, which is the survival rate. So in () we need t to be the number of cases.

We set up the design matrix, with indicators 1 for adults and males, and using indicator variables for second class and third class, which means first class will be a reference.

ones=ConstantArray[1,12];
age={0,0,0,0,0,0,1,1,1,1,1,1};
sex={0,0,0,1,1,1,0,0,0,1,1,1};
class2={0,1,0,0,1,0,0,1,0,0,1,0};
class3={1,0,0,1,0,0,1,0,0,1,0,0};
X2=Join[Thread[List[ones,age,sex,class2,class3]]];

Then we set up the dependent variable and the offset.

y2={14,13,1,13,11,5,76,80,140,75,14,57};
t2={31,13,1,48,11,5,165,93,144,462,168,175};

We define the log-likelihood ().

lnL2[X_,y_,t_,α_,β_]:=
Module[{n=Length[y]},(∑_(i=1)^n)(y[[i]] Log[α]+y[[i]] (X[[i]].β+Log[t[[i]]])-(y[[i]]+(1/α))Log[1+α (E^(X[[i]].β+Log[t[[i]]]))]+LogGamma[y[[i]]+(1/α)]-LogGamma[y[[i]]+1]-LogGamma[(1/α)])
]

Now we maximize it to find the coefficients.

results2=FindMaximum[lnL2[X2,y2,t2,α,β],{{α,1.0,0.01,4.0},{β,{1.0,1.0,1.0,1.0,1.0}}},AccuracyGoal->6]

Then we find the standard errors of the coefficients.

VarianceCovarianceMatrix2[X_,y_,results_]:=-Inverse[HessianH[lnL2[X,y,t2,α,β]/.β->βs[X],{α,βs[X]}]/.flattenrules[results]]

StdErr2=(Sqrt[Diagonal[VarianceCovarianceMatrix2[X2,y2,results2]]])

And again we can print a table of the results.

CoefficientsTable[Last/@flattenrules[results2],StdErr2,Prepend[βs[X2],α]]

But perhaps more useful for interpretation of the coefficients would be the Incidence Rate Ratio (IRR) for each variable, which is obtained by exponentiating each coefficient. For example, out of a sample of t adults, we expect that the survival rate, from our model (), will be (Subscript[μ, adults])/t=exp((Subscript[β, 0])+(Subscript[β, 1])(1)+(Subscript[β, 2])(sex)+(Subscript[β, 2])(class2)+(Subscript[β, 3])(class3)), while for an identical number t of children we expect their survival rate to be (Subscript[μ, children])/t=exp((Subscript[β, 0])+(Subscript[β, 1])(0)+(Subscript[β, 2])(sex)+(Subscript[β, 2])(class2)+(Subscript[β, 3])(class3)). So by dividing the two rates, we obtain the ratio of rates (IRR) to be

IRR=(((Subscript[μ, adults])/t)/((Subscript[μ, children])/t))=((exp((Subscript[β, 0])+(Subscript[β, 1])(1)+(Subscript[β, 2])(sex)+(Subscript[β, 2])(class2)+(Subscript[β, 3])(class3)))/(exp((Subscript[β, 0])+(Subscript[β, 1])(0)+(Subscript[β, 2])(sex)+(Subscript[β, 2])(class2)+(Subscript[β, 3])(class3))))=(e^(Subscript[β, 1])),

which we estimate to be (e^(-0.670034))=0.51. Thus, our interpretation is that adults survived at roughly half the rate at which children survived, among those of the same sex and class. The standard error of IRR is found by multiplying the estimated IRR by the standard error of the coefficient (see []), while a confidence interval for IRR is found by exponentiating the confidence interval for the coefficient. Thus we obtain the following.

coefficients2=Last/@flattenrules[results2];
lower2=coefficients2-1.96 StdErr2;
upper2=coefficients2+1.96 StdErr2;

IRR=(E^coefficients2);
IRRStdErr=IRR StdErr2;
IRRlower=(E^lower2);
IRRupper=(E^upper2);

We do not need IRR for α or (Subscript[β, 0]), so we drop them and then print the resulting table.

IRR=Drop[IRR,2];
IRRStdErr=Drop[IRRStdErr,2];
IRRlower=Drop[IRRlower,2];
IRRupper=Drop[IRRupper,2];
Text@TableForm[{IRR,IRRStdErr,IRRlower,IRRupper}//Transpose,TableHeadings->{{"age","sex","class2","class3"},{"IRR","IRR\nStd. Err.","95% IRR\nCI lower","95% IRR\nCI upper"}}]

The confidence interval for the variable class2 contains 1.0, consistent with the lack of significance of its coefficient, and indicating that the survival rate of second-class passengers was not significantly different than that of first-class passengers. We will address this after computing some model assessment statistics and residuals.

Model Assessment

Various types of model fit statistics and residuals are readily computed. We use definitions given in []; alternate definitions exist and would require only minor changes.

Commonly used model fit statistics include the log-likelihood, deviance, Pearson chi-square dispersion, Akaike Information Criterion (AIC), and Bayesian Information Criterion (BIC).

We already have the log-likelihood ℒ as a byproduct of the maximization process. The deviance D is defined as

D=2(∑_(i=1)^n)(ℒ((Subscript[y, i]);(Subscript[y, i]))-ℒ((Subscript[μ, i]);(Subscript[y, i]))),

where ℒ((Subscript[μ, i]);(Subscript[y, i])) is our log-likelihood function (), and ℒ((Subscript[y, i]);(Subscript[y, i])) is the log-likelihood function with (Subscript[y, i]) replacing (Subscript[μ, i]). For our NB2 model, this simplifies to D=(∑_(i=1)^n)((Subscript[d, i])^2), where

((Subscript[d, i])^2)=({{{, ({{2((Subscript[y, i]) ln(((Subscript[y, i])/(Subscript[μ, i])))-((Subscript[y, i])+(1/α))ln(((1+α (Subscript[y, i]))/(1+α (Subscript[μ, i]))))), if (Subscript[y, i])>0}, {(2/α)ln(1+α (Subscript[μ, i])), if (Subscript[y, i])=0.}})}})

The Pearson chi-square dispersion statistic is given by (∑_(i=1)^n)((((Subscript[y, i])-(Subscript[μ, i])))^2)/((Subscript[μ, i])+α ((Subscript[μ, i])^2)), while AIC and BIC are defined as

AIC=((2(ℒ-(p+1)))/n)

and

BIC=D-(n-p-1)ln n.

We compute these for the Titanic data above and display them.

dSquared[i_,y_,μ_,Α_]:=(Piecewise[{{2(y[[i]] Log[((y[[i]])/(μ[[i]]))]-(y[[i]]+(1/Α))Log[((1+Α y[[i]])/(1+Α μ[[i]]))]), y[[i]]≠0}, {(2/Α)Log[1+Α μ[[i]]], y[[i]]==0}}]);

ModelAssessment[coefficients_,X_,y_,t_]:=Module[{μ,n,results,deviance,Pearson,p,AIC,BIC},
μ=(E^(X.Drop[coefficients,1]+Log[t]));
n=Length@y;
results=FindMaximum[lnL2[X,y,t,α,β],{{α,1.0,0.01,4.0},{β,Table[1.0,{Length@X[[1]]}]} },AccuracyGoal->6];
ℒ=results[[1]];
deviance=(∑_(i=1)^n)dSquared[i,y,μ,coefficients[[1]]];
Pearson=(∑_(i=1)^n)((((y[[i]]-μ[[i]]))^2)/(μ[[i]]+coefficients[[1]] ((μ[[i]])^2)));
p=Length[X[[2]]]-1;
AIC=-((2(ℒ-(p+1)))/n);
BIC=deviance-(n-p-1) Log[n];
{ℒ,deviance,Pearson,AIC,BIC}
]

ModelAssessmentTable[coefficients_,X_,y_,t_]:=Text@TableForm[ModelAssessment[coefficients,X,y,t],TableHeadings->{{"Log\[Hyphen]likelihood","Deviance","Pearson","AIC","BIC"}}]

ModelAssessmentTable[coefficients2,X2,y2,t2]

These model assessment statistics are most useful when compared to those of a competing model, which we pursue in the next section after computing residuals.

Residuals

The raw residuals are of course (Subscript[y, i])-(Subscript[μ, i]), while the Pearson residuals are ((Subscript[y, i])-(Subscript[μ, i]))/(Sqrt[(Subscript[μ, i])+α ((Subscript[μ, i])^2)]), and the deviance residuals are (sgn((Subscript[y, i])-(Subscript[μ, i])))(Subscript[d, i]), as defined in ().

These residuals can be standardized by dividing by (Sqrt[1-(Subscript[h, i])]), where the (Subscript[h, i]) are the leverages obtained from the diagonal of the hat matrix (W^(1/2)) ((X(X' W X))^(-1))X' (W^(1/2)), for W equal to the n×n diagonal matrix, with (Subscript[μ, i])/(1+α (Subscript[μ, i])) as the (i^th) element of the diagonal.

Here are the unstandardized residuals for the Titanic data.

UnstandardizedResiduals[coefficients_,X_,y_,t_]:=Module[{rawResiduals,μ,PearsonResiduals,DevianceResiduals},
μ=(E^(X.Drop[coefficients,1]+Log[t]));
rawResiduals=y-μ;
PearsonResiduals=((y-μ)/(Sqrt[μ+coefficients[[1]] (μ^2)]));
DevianceResiduals=Table[Sign[y2[[i]]-μ[[i]]](Sqrt[dSquared[i,y,μ,coefficients[[1]]]]),{i,Length@y}];
{rawResiduals,PearsonResiduals,DevianceResiduals}
]

UnstandardizedResidualsTable[coefficients_,X_,y_,t_]:=Text@TableForm[UnstandardizedResiduals[coefficients,X,y,t]//Transpose,TableHeadings->{{},{"Raw\nresiduals","Pearson\nresiduals","Deviance \nresiduals"}}]

UnstandardizedResidualsTable[coefficients2,X2,y2,t2]

And here are the leverages and the standardized residuals.

LeveragesAndStandardizedResiduals[coefficients_,X_,y_,t_]:=Module[{μ,h,PearsonResiduals,DevianceResiduals,StandardizedPearsonResiduals,StandardizedDevianceResiduals},
μ=(E^(X.Drop[coefficients,1]+Log[t]));
W=DiagonalMatrix[(μ/(1+coefficients[[1]] μ))];
(*W=DiagonalMatrix[(μ/(1+α μ))]/.α->coefficients[[1]];*)
h=Diagonal[(W^(1/2)).X.Inverse[Transpose[X].W.X].Transpose[X].(W^(1/2))] ;
{PearsonResiduals,DevianceResiduals}=Rest@UnstandardizedResiduals[coefficients,X,y,t];
StandardizedPearsonResiduals=(PearsonResiduals/(Sqrt[1-h]));
StandardizedDevianceResiduals=(DevianceResiduals/(Sqrt[1-h]));
{h,StandardizedPearsonResiduals,StandardizedDevianceResiduals}
]

LeveragesAndStandardizedResidualsTable[coefficients_,X_,y_,t_]:=Text@TableForm[LeveragesAndStandardizedResiduals[coefficients,X,y,t]//Transpose,TableHeadings->{{},{Row[{"Leverages\n(",Style["h",Italic],")"}],"Standardized\nPearson\nresiduals","Standardized\nDeviance\nresiduals"}}]

LeveragesAndStandardizedResidualsTable[coefficients2,X2,y2,t2]

Hilbe recommends plotting the Standardized Pearson residuals versus h, with a poor model fit indicated by residuals that are outside the interval ±2 when the leverage is high.

ListPlot[Most@LeveragesAndStandardizedResiduals[coefficients2,X2,y2,t2]//Transpose,AxesLabel->{Style["h",Italic],"Std. Pearson Residuals"}]

We have two Standardized Pearson residuals that are not within the range ±2, one of which has a high leverage. We also recall that the variable class2 was not significant. Perhaps the model will be improved if we remove class2. All that is required is to remove class2 from the design matrix X, remove the corresponding starting value from the maximizing command, and run the model again. We obtain the following assessment statistics and standardized residuals for the revised model with class2 removed.

We set up design matrix X and find the coefficients.

X3=Join[Thread[List[ones,age,sex,class3]]];results3=FindMaximum[lnL2[X3,y2,t2,α,β],{{α,1.0,0.01,4.0},{β,{1.0,1.0,1.0,1.0}}},AccuracyGoal->6];
coefficients3=Last/@flattenrules[results3]

ModelAssessmentTable[coefficients3,X3,y2,t2]

LeveragesAndStandardizedResidualsTable[coefficients3,X3,y2,t2]

ListPlot[Most@LeveragesAndStandardizedResiduals[coefficients3,X3,y2,t2]//Transpose,AxesLabel->{Style["h",Italic],"Std. Pearson Residuals"}]

Comparing to the full model, we see that the assessment statistics have improved (they are smaller, indicating a better fit), and the Standardized Pearson residuals with high leverages are within the recommended boundaries. It appears that the model has been improved by dropping class2.

Conclusion

The traditional negative binomial regression model (NB2) was implemented by maximum likelihood estimation without much difficulty, thanks to the maximization command and especially to the automatic computation of the standard errors via the Hessian.

Other negative binomial models, such as the zero-truncated, zero-inflated, hurdle, and censored models, could likewise be implemented by merely changing the likelihood function.

Acknowledgments

The author acknowledges suggestions and assistance by the editor and the referee that helped to improve this article.

References

J. Hilbe, Negative Binomial Regression, 2nd ed., New York: Cambridge University Press, 2011.

“JSE Data Archive.” Journal of Statistics Education. (Nov 19, 2012) www.amstat.org/publications/jse/jse_data_archive.htm.