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.
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
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.
As to your other questions...
Yes, it is more than appropriate, it is expected. One explanation is found here.
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.
Negative Binomial regression does not predict probabilities. Are you sure you're not talking about logistic regression? In any case, both age and age squared are used in making predictions.
