# Interpret negative binomial coefficients and p-value

I fit a negative binomial model for my dataset because of over-dispersion in the poisson model.

DOWFACAT and SSNVOLCLS are categorical. All other variables are numeric.

MDHabitat, MigCorr,WintRg, SumRg, Fence05, Fence16 and Cross05 are numeric with either a 0 or 1 to indicate absence or presence of that variable.

PostFFall is the count

AADTPostFal is the average annual daily traffic/1000

SpeedLmt is numeric and ranges from 0-75

Using the summary function I get these results:

glm.nb(formula = PostFFall ~ AADTPostFal + SpeedLmt + SumRg +
WintRg + MigCorr + MDHabitat + DOWFACAT + SSNVOLCLS + Fence05 +
Cross05 + Fence16, data = roads1, init.theta = 0.4313904984,

Deviance Residuals:
Min        1Q    Median        3Q       Max
-2.38236  -0.05334   0.00000   0.00000   3.04227

Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept)          -3.523e+01  5.275e+06   0.000   1.0000
SpeedLmt             -2.193e-02  1.226e-02  -1.789   0.0737 .
SumRg                 4.842e-01  6.739e-01   0.719   0.4724
WintRg                4.410e-01  6.148e-01   0.717   0.4732
MigCorr               1.753e+00  6.825e-01   2.568   0.0102 *
MDHabitat             3.225e+01  1.677e+06   0.000   1.0000
DOWFACATSuburban     -3.389e-01  9.588e-01  -0.353   0.7238
DOWFACATTransition    1.120e+00  1.192e+00   0.939   0.3475
DOWFACATUrban         8.394e-01  1.326e+07   0.000   1.0000
SSNVOLCLS1-LowVolume -2.809e+00  5.535e+06   0.000   1.0000
SSNVOLCLS2-MidVolume  2.914e+00  5.535e+06   0.000   1.0000
Fence05              -5.679e-01  1.190e+00  -0.477   0.6331
Cross05              -3.945e+01  3.460e+07   0.000   1.0000
Fence16               4.445e+00  1.066e+00   4.169 3.05e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

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

Null deviance: 382.555  on 1641  degrees of freedom
Residual deviance:  82.308  on 1627  degrees of freedom
AIC: 243.56

Number of Fisher Scoring iterations: 1

Theta:  0.431
Std. Err.:  0.143

2 x log-likelihood:  -211.560


And then when I run an ANOVA test, most of the p-values are significant:

Analysis of Deviance Table

Response: PostFFall

Terms added sequentially (first to last)

Df Deviance Resid. Df Resid. Dev  Pr(>Chi)
NULL                         1641     382.55
AADTPostFal  1   25.899      1640     356.66 3.598e-07 ***
SpeedLmt     1    0.793      1639     355.86  0.373231
SumRg        1  109.769      1638     246.09 < 2.2e-16 ***
WintRg       1    8.415      1637     237.68  0.003721 **
MigCorr      1   29.217      1636     208.46 6.470e-08 ***
MDHabitat    1   79.120      1635     129.34 < 2.2e-16 ***
DOWFACAT     3    1.575      1632     127.77  0.664987
SSNVOLCLS    2   27.064      1630     100.70 1.328e-06 ***
Fence05      1    0.740      1629      99.96  0.389572
Cross05      1    1.947      1628      98.02  0.162953
Fence16      1   15.708      1627      82.31 7.392e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1


Do I use the anova test to decide whether the model variables are significant? Also what does a negative coefficient intercept indicate? Also, my Dispersion parameter for Negative Binomial is (0.4314), does this mean my model is now under-dispersed?

It seems you are asking about the discrepancy between the table of coefficients from the summary() function and the analysis of variance table from anova(). In particular, you are wondering why the anova table gives small p-values to variables that have p-value close to 1 in the table of coefficients (e.g. variable sumRg).

The explanation may be that the table of coefficients tells you whether a variable or level of a factor is different from the reference baseline (i.e. from the intercept as explained by Isabella Ghement's answer) while the anova table tells you whether that variable contributes significantly to the quality of the fit.

I suspect that you are more interested in the output of the anova table since you want to select relevant variables. Note, however, that the anova function estimates the contribution of each variable by adding them sequentially in the same order as you listed them in the model. Usually, one wants to estimate the contribution of a variable independently from the order they enter into the model. For this you could use the Anova function in the car package.

As a comment, it looks some variables like MDHabitat have huge standard error relative to the estimated effect (e.g. est= 3.225e+01, se= 1.677e+06). I wonder whether it indicates some weird characteristic of the data that could/should be addressed.

• I guess I do want to use the anova then based on your comment. The MDhabitat variable is just a 0 or 1 for presence and absence. I'm not sure what I can do to address this. Do you have any pointers for addressing standard errors that are larger relative to the estimate effect? Also, my Dispersion parameter for Negative Binomial is (0.4314), does this mean my model is now under-dispersed? – JNN Dec 4 '20 at 17:34

A negative binomial regression model presuposes the dependent variable is a count variable (usually collected over the same units of time or space for each unit in a study; if that is not the case, the model would need to include an offset term). For the purposes of this answer, I will assume the model does not include an offset.

The model relates the log expected value of the count dependent variable to one or more predictor variables. (Another word for expected is average.) Thus, the value of the intercept in the model stands for the log expected value of the count dependent variable when the values of all predictor variables in your dataset were set to 0.

R uses dummy coding for categorical predictor variables included in your model (i.e., those converted to a factor using the factor() function). Setting the values of these categorical predictor variables to 0 when interpreting the intercept therefore means setting their value to the reference category (i.e., the category against which all other categories of the same variable will be compared). You can easily find out the reference category for a categorical predictor variable in R by applying the function levels() to that variable - the level (aka category) listed first by R is treated as the reference category. (Some people might call this the baseline category.) Note that, by default, R sets the reference category of a categorical predictor to the the category that appears first in alphabetical order among all categories of that predictor - this may or may not sense in your case. If it doesn't make sense, you can use the relevel() function of R to chance the reference category to something more meaningful.

For example, if you had a negative binomial regression model which related number of birds (dependent variable) to habitat type (top of mountain, middle of mountain or bottom of mountain), then the intercept in this model would be interpreted as the log expected bird count for bottom of mountain, provided bottom of mountain was declared as the reference category for habitat type.

Setting numeric predictors in your model to 0 may or may not lead to a meaningful interpretation for the intercept. If you have a Year predictor in your model, with values 0, 1, 2, etc., setting year to 0 will enable you to talk about the log expected value of the dependent variable in the first year of study (i.e., year 0). If you have a predictor like "bird size", then setting that predictor to 0 when interpreting the intercept will lead to a non-sensical interpretation - you can't have birds with a size of 0! In this case, you may want to center your variable around its mean or median to ensure a 0 value for the centered variable is meaningful.

The binary variables of the form "presence/absence" do not need any centering. If you convert them to factors in your model, just set their reference level to whatever is meaningful to you (presumably "absence").

If your count predictor variables take the value 0 in your data, you don't need to center them - you can interpret the value of the intercept as the log expected value of the count response variable when those count predictor variables are equal to 0 and all other predictors are equal to 0.

If you have a count predictor variable which does NOT take the value 0 in your data (e.g., its values are greater than 10), you will need to center it. How you center it depends on what the distribution of that variable looks like when you plot it. If it looks approximately normal, you will center the count predictor variable around its mean value observed in the data. If it looks unimodal and skewed, you will center it around its median value observed in the data.

In a negative binomial modelling context, centering a predictor variable $$X$$ around its mean involves replacing $$X$$ with $$Xcen = X - mean(X)$$ and then re-fitting your model with $$Xcen$$ instead of $$X$$. In a model which uses $$Xcen$$ instead of $$X$$, you would interpret the value of the intercept as the log expected value of the count response variable when $$Xcen = 0$$ and all other predictor variables are 0. Since $$Xcen = 0$$ when $$X = mean(X)$$ (presuming you centered $$X$$ around its mean), this is the same as saying that in the model which used $$Xcen$$ as a predictor, the intercept is the log expected value of the count response variable when $$X = mean(X)$$ (i.e., when $$X$$ is set to its average value observed in the data) and all other predictor variables are set to 0.

A similar principe applies for continuous predictor varibles when it comes to centering - just plot their distribution to see if it makes sense to center these variables around their mean or median values observed in the data. Alternatively, if there is a value of one of these variables that is substantively of interest, you can center around that value. As an example, let's say that $$X$$ stands for Age and the age 18 is important in your research as it denotes the Age when your study subjects can begin to vote - then your "centered" Age would be computed as Age - 18 and the intercept in a model with centered Age would be interpreted accordingly.

Scientific notation used by R can make it difficult to see if you have really large or really small standard errors in your model summary. You can turn it off using:

options(scipen=999)


and then print your model summary to see what happens. If you have really large standard errors for some of the dummy variables used to encode the effect of the factors (i.e., categorical variables) in your model, that could be a signal that the categories represented by those dummy variables are sparse (i.e., there ar very few of those categories represented in the data). If that is the case, you will need to merge the sparse categories with other categories that are better represented in the data or perhaps even drop the sparse categories from your data - whichever makes most sense in your context.

• That was really helpful, thank you. I edited my post with a description of what my variables look like. Most are presence (1) and absence (0) with the exception of AADTPostFal, SpeedLmt, and PostFal (the total count). Could you explain a little more on how I would center my variables with this in mind? – JNN Dec 4 '20 at 20:03
• You’re welcome, @James. See my Addendum 1 and Addendum 2 for some further clarifications. – Isabella Ghement Dec 6 '20 at 3:30