# Interpreting Quadratic Variables in Negative Binomial Regression

I want to include a quadratic term for age in my negative binomial model as past work has suggested it may be curvilinear. I know how to interpret a quadratic coefficient in OLS, but am unsure with negative binomial.

1. Is it appropriate to include both the original variable (age) and the quadratic (agequad) in the model?
2. How are the coefficients interpreted?
3. When plotting predicted probabilities, is the original (age) or quadratic (agequad) used? Note: I use stata's margins command for this.
• 1. What link function are you using in your negative binomial model? 2. Note that if it's a non-identity link you'll already have a curved relationship between the conditional mean response and age. 3. It's gong to be hard for people to discuss appropriateness when you don't even mention your response variable. 4. Are there other variables besides age in the model? Commented Jan 4, 2022 at 6:22

Generally, if you add high-order polynomial terms, you should also include the lower-level ones (i.e., both $$x$$ and $$x^2)$$. The polynomial terms are used to approximate a complex nonlinear relationship with a Taylor series expansion of the unknown nonlinear function.

You should also let Stata know that these variables are related by using factor variable notation rather than a separate variable. Otherwise, Stata does not know that agequad and age are related and most of the calculations will be incorrect.

Here is an example of a negative binomial model for the number of doctor visits as a function of age and age$$^2$$ and a constant:

. use http://cameron.econ.ucdavis.edu/nhh2017/mus17data.dta, clear

. nbreg docvis c.age##c.age, nolog

Negative binomial regression                            Number of obs =  3,677
LR chi2(2)    =  29.52
Dispersion: mean                                        Prob > chi2   = 0.0000
Log likelihood = -10961.297                             Pseudo R2     = 0.0013

------------------------------------------------------------------------------
docvis | Coefficient  Std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
age |   .2989626   .0660418     4.53   0.000     .1695231    .4284022
|
c.age#c.age |  -.0019345   .0004394    -4.40   0.000    -.0027958   -.0010733
|
_cons |  -9.537893   2.467098    -3.87   0.000    -14.37332   -4.702469
-------------+----------------------------------------------------------------
/lnalpha |  -.1833182   .0282197                     -.2386278   -.1280086
-------------+----------------------------------------------------------------
alpha |   .8325032    .023493                       .787708    .8798458
------------------------------------------------------------------------------
LR test of alpha=0: chibar2(01) = 1.2e+04              Prob >= chibar2 = 0.000


The index function coefficients are hard to interpret on their own, especially with interactions since $$E[y \vert x] = \exp \{a+b\cdot age + c \cdot age^2\}$$. It is much easier to use margins:

. margins, dydx(age)

Average marginal effects                                 Number of obs = 3,677
Model VCE: OIM

Expression: Predicted number of events, predict()
dy/dx wrt:  age

------------------------------------------------------------------------------
|            Delta-method
|      dy/dx   std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
age |   .0714711   .0153976     4.64   0.000     .0412924    .1016499
------------------------------------------------------------------------------


This means that if everyone in your sample was 1 year older, that would be associated with a .07 additional doctor visits. This is called an average marginal effect. It's the average derivative of the expected value with respect to age:

$$AME_{age} = \frac{1}{N} \sum_i^N \exp \{ a+b\cdot age_i + c \cdot age_i^2\} \cdot (b+ 2 \cdot age_i)$$

You can also calculate predictions of various sorts (code below):

The first graph shows that there is an inverted-U relationship between age and expected visits to the doc. Not surprising, the slope in the graph below is at first positive and then turns negative somewhere between 75 and 80. Since most people in this dataset are relatively young, it is not surprising that the overall effect we got above is positive.

You can also use margins to calculate various predictions about the number of visits (second column).

All the calculations include both the effect of age and age$$^2$$.

Stata Code:

use http://cameron.econ.ucdavis.edu/nhh2017/mus17data.dta, clear
nbreg docvis c.age##c.age, nolog
margins, dydx(age)
margins, at(age = (65(5)90))
quietly marginsplot, name(yhat, replace) title("Predicted Visits")
margins, dydx(age) at(age = (65(5)90))
quietly marginsplot, name(slope, replace) title("Marginal Effect on Visits")
margins, predict(pr(7)) at(age = (65(5)90))
quietly marginsplot, name(phat7, replace) title("Probability of Exactly 7 Visits")
margins, predict(pr(7,.)) at(age = (65(5)90))
quietly marginsplot, name(phat7plus, replace) title("Probability of 7+ Visits")
graph combine yhat phat7 slope phat7plus, xcommon

• “The polynomial terms are used to approximate a complex nonlinear relationship with a Taylor series expansion of the unknown nonlinear function.”. This may not be as well-known as it should be. Adding lower order terms if higher order terms are present is often mentioned, but a rationale is too often unstated. Commented Feb 15, 2022 at 9:38

I want to include a quadratic term for age in my negative binomial model as past work has suggested it may be curvilinear.

While I do not doubt the conditional mean is non-linear in age, do you have any specific reason to prefer a quadratic term as opposed to something more flexible such as a natural spline? Polynomials are great tools, but are rather biased because you can only estimate things which look like quadratic functions (or whatever polynomial you decide to use). If you require additional flexibility, a spline is a more flexible but still interpretable method.

2. The coefficients are not directly interpretable, you have to look at the conditional mean in order to interpret the effects of age. In a negative binomial regression, your model will look like $$\log(E(y)) = \beta_0 + \beta_1 x + \beta_2 x^2$$ The typical "a one unit change in $$x$$ leads to a $$\beta$$ unit change in the expectation of the outcome" no longer applies because of the included quadratic term. You could take a partial derivative with respect to $$x$$ in order to determine how varying $$x$$ changes the expectation on the log scale, but my preference is just to plot the estimated function. If all other predictors are linear, then they will simply move the function up or down.