## Effects of nonlinear models with interactions of discrete and continuous variables: Estimating, graphing, and interpreting

I want to estimate, graph, and interpret the effects of nonlinear models with interactions of continuous and discrete variables. The results I am after are not trivial, but obtaining what I want using **margins**, **marginsplot**, and factor-variable notation is straightforward.

**Do not create dummy variables, interaction terms, or polynomials**

Suppose I want to use **probit** to estimate the parameters of the relationship

\begin{equation*}

P(y|x, d) = \Phi \left(\beta_0 + \beta_1x + \beta_3d + \beta_4xd + \beta_2x^2 \right)

\end{equation*}

where \(y\) is a binary outcome, \(d\) is a discrete variable that takes on four values, \(x\) is a continuous variable, and \(P(y|x,d)\) is the probability of my outcome conditional on covariates. To fit this model in Stata, I would type

probit y c.x##i.d c.x#c.x

I do not need to create variables for the polynomial or for the interactions between the continuous variable \(x\) and the different levels of \(d\). Stata understands that **c.x#c.x** is the square of \(x\) and that **c.x##i.d** corresponds to the variables \(x\) and \(d\) and their interaction. The result of what I typed would look like this:

. probit y c.x##i.d c.x#c.x Iteration 0: log likelihood = -686.40401 Iteration 1: log likelihood = -643.73527 Iteration 2: log likelihood = -643.44789 Iteration 3: log likelihood = -643.44775 Iteration 4: log likelihood = -643.44775 Probit regression Number of obs = 1,000 LR chi2(8) = 85.91 Prob > chi2 = 0.0000 Log likelihood = -643.44775 Pseudo R2 = 0.0626 ---------------------------------------------------------------------------- y | Coef. Std. Err. z P>|z| [95% Conf. Interval] -------------+-------------------------------------------------------------- x | .3265507 .0960978 3.40 0.001 .1382025 .5148988 | d | uno | .3806342 .1041325 3.66 0.000 .1765382 .5847302 dos | .5815309 .1177613 4.94 0.000 .350723 .8123388 tres | .6864413 .1981289 3.46 0.001 .2981157 1.074767 | d#c.x | uno | -.3163965 .1150582 -2.75 0.006 -.5419065 -.0908865 dos | -.5075077 .1300397 -3.90 0.000 -.7623809 -.2526345 tres | -.7081746 .2082792 -3.40 0.001 -1.116394 -.2999549 | c.x#c.x | -.1870304 .0339825 -5.50 0.000 -.2536349 -.1204259 | _cons | -.0357237 .0888168 -0.40 0.688 -.2098015 .1383541 ----------------------------------------------------------------------------

I did not need to create dummy variables, interaction terms, or polynomials. As we will see below, convenience is not the only reason to use factor-variable notation. Factor-variable notation allows Stata to identify interactions and to distinguish between discrete and continuous variables to obtain correct marginal effects.

This example used **probit**, but most of Stata’s estimation commands allow the use of factor variables.

**Using margins to obtain the effects I am interested in**

I am interested in modeling for individuals the probability of being married (**married**) as a function of years of schooling (**education**), the percentile of income distribution to which they belong (**percentile**), the number of times they have been divorced (**divorce**), and whether their parents are divorced (**pdivorce**). I estimate the following effects:

- The average of the change in the probability of being married when each covariate changes. In other words, the average marginal effect of each covariate.
- The average of the change in the probability of being married when the interaction of
**divorce**and**education**changes. In other words, an average marginal effect of an interaction between a continuous and a discrete variable. - The average of the change in the probability of being married when the interaction of
**divorce**and**pdivorce**changes. In other words, an average marginal effect of an interaction between two discrete variables. - The average of the change in the probability of being married when the interaction of
**percentile**and**education**changes. In other words, an average marginal effect of an interaction between two continuous variables.

I fit the model:

probit married c.education##c.percentile c.education#i.divorce /// i.pdivorce##i.divorce

The average of the change in the probability of being married when the levels of the covariates change is given by

. margins, dydx(*) Average marginal effects Number of obs = 5,000 Model VCE : OIM Expression : Pr(married), predict() dy/dx w.r.t. : education percentile 1.divorce 2.divorce 1.pdivorce ---------------------------------------------------------------------------- | Delta-method | dy/dx Std. Err. z P>|z| [95% Conf. Interval] -------------+-------------------------------------------------------------- education | .02249 .0023495 9.57 0.000 .017885 .027095 percentile | .9873168 .062047 15.91 0.000 .8657069 1.108927 | divorce | 1 | -.0434363 .0171552 -2.53 0.011 -.0770598 -.0098128 2 | -.1239932 .054847 -2.26 0.024 -.2314913 -.0164951 | 1.pdivorce | -.0525977 .0131892 -3.99 0.000 -.0784482 -.0267473 ---------------------------------------------------------------------------- Note: dy/dx for factor levels is the discrete change from the base level.

The first part of the **margins** output states the statistic it is going to compute, in this case, the average marginal effect. Next, we see the concept of an **Expression**. This is usually the default prediction (in this case, the conditional probability), but it can be any other prediction available for the estimator or any function of the coefficients, as we will see shortly.

When **margins** computes an effect, it distinguishes between continuous and discrete variables. This is fundamental because a marginal effect of a continuous variable is a derivative, whereas a marginal effect of a discrete variable is the change of the **Expression** evaluated at each value of the discrete covariate relative to the **Expression** evaluated at the base or reference level. This highlights the importance of using factor-variable notation.

I now interpret a couple of the effects. On average, a one-year change in education increases the probability of being married by 0.022. On average, the probability of being married is 0.043 smaller in the case where everyone has been divorced once compared with the case where no one has ever been divorced, an average treatment effect of \(-\)0.043. The average treatment effect of being divorced two times is \(-\)0.124.

Now I estimate the average marginal effect of an interaction between a continuous and a discrete variable. Interactions between continuous and discrete variables are changes in the continuous variable evaluated at the different values of the discrete covariate relative to the base level. To obtain these effects, I type

. margins divorce, dydx(education) pwcompare Pairwise comparisons of average marginal effects Model VCE : OIM Expression : Pr(married), predict() dy/dx w.r.t. : education -------------------------------------------------------------- | Contrast Delta-method Unadjusted | dy/dx Std. Err. [95% Conf. Interval] -------------+------------------------------------------------ education | divorce | 1 vs 0 | .0012968 .0061667 -.0107897 .0133833 2 vs 0 | .0403631 .0195432 .0020591 .0786672 2 vs 1 | .0390664 .0201597 -.0004458 .0785786 --------------------------------------------------------------

The average marginal effect of education is 0.039 higher when everyone is divorced two times instead of everyone being divorced one time. The average marginal effect of education is 0.040 higher when everyone is divorced two times instead of everyone being divorced zero times. The average marginal effect of education is 0 when everyone is divorced one time instead of everyone being divorced zero times. Another way of obtaining this result is by computing a cross or double derivative. As I mentioned before, we use derivatives for continuous variables and differences with respect to the base level for the discrete variables. I will refer to them loosely as derivatives hereafter. In the appendix, I show that taking a double derivative is equivalent to what I did above.

Analyzing the interaction between two discrete variables is similar to analyzing the interaction between a discrete and a continuous variable. We want to see the change from the base level of a discrete variable for a change in the base level of the other variable. We use the **pwcompare** and **dydx()** options again.

. margins pdivorce, dydx(divorce) pwcompare Pairwise comparisons of average marginal effects Model VCE : OIM Expression : Pr(married), predict() dy/dx w.r.t. : 1.divorce 2.divorce -------------------------------------------------------------- | Contrast Delta-method Unadjusted | dy/dx Std. Err. [95% Conf. Interval] -------------+------------------------------------------------ 1.divorce | pdivorce | 1 vs 0 | -.1610815 .0342655 -.2282407 -.0939223 -------------+------------------------------------------------ 2.divorce | pdivorce | 1 vs 0 | -.0201789 .1098784 -.2355367 .1951789 -------------------------------------------------------------- Note: dy/dx for factor levels is the discrete change from the base level.

The average change in the probability of being married when everyone is once divorced and everyone’s parents are divorced, compared with the case where no one’s parents are divorced and no one is divorced, is a decrease of 0.161. The average change in the probability of being married when everyone is twice divorced and everyone’s parents are divorced, compared with the case where no one’s parents are divorced and no one is divorced, is 0. We could have obtained the same result by typing **margins divorce, dydx(pdivorce) pwcompare**, which again emphasizes the concept of a cross or double derivative.

Now I look at the average marginal effect of an interaction between two continuous variables.

. margins, > expression(normalden(xb())*(_b[percentile] + > education*_b[c.education#c.percentile])) > dydx(education) Average marginal effects Number of obs = 5,000 Model VCE : OIM Expression : normalden(xb())*(_b[percentile] + education*_b[c.education#c.percentile]) dy/dx w.r.t. : education ---------------------------------------------------------------------------- | Delta-method | dy/dx Std. Err. z P>|z| [95% Conf. Interval] -------------+-------------------------------------------------------------- education | .0619111 .021641 2.86 0.004 .0194955 .1043267 ----------------------------------------------------------------------------

The **Expression** is the derivative of the conditional probability with respect to **percentile**. **dydx(education)** specifies that I want to estimate the derivative of this **Expression** with respect to education. The average marginal effect in the marginal effect of income percentile attributable to a change in education is 0.062.

Because **margins** can only take first derivatives of expressions, I obtained a cross derivative by making the expression a derivative. In the appendix, I show the equivalence between this strategy and writing a cross derivative. Also, I illustrate how to verify that your expression for the first derivative is correct.

**Graphing**

After **margins**, we can plot the results in the output table simply by typing **marginsplot**. **marginsplot** works with the conventional graphics options and the Graph Editor. For the first example above, for instance, I could type

. quietly margins, dydx(*) . marginsplot, xlabel(, angle(vertical)) Variables that uniquely identify margins: _deriv

I added the option **xlabel(, angle(vertical))** to obtain vertical labels for the horizontal axis. The result is as follows:

**Conclusion**

I illustrated how to compute, interpret, and graph marginal effects for nonlinear models with interactions of discrete and continuous variables. To interpret interaction effects, I used the concepts of a cross or double derivative and an **Expression**. I used simulated data and the probit model for my examples. However, what I wrote extends to other nonlinear models.

**Appendix**

To verify that your expression for the first derivative is correct, you compare it with the statistic computed by **margins** with the option **dydx(***variable***)**. For the example in the text:

. margins, > expression(normalden(xb())*(_b[percentile] + > education*_b[c.education#c.percentile])) Predictive margins Number of obs = 5,000 Model VCE : OIM Expression : normalden(xb())*(_b[percentile] + education*_b[c.education#c.percentile]) ---------------------------------------------------------------------------- | Delta-method | Margin Std. Err. z P>|z| [95% Conf. Interval] -------------+-------------------------------------------------------------- _cons | .9873168 .062047 15.91 0.000 .8657069 1.108927 ---------------------------------------------------------------------------- . margins, dydx(percentile) Average marginal effects Number of obs = 5,000 Model VCE : OIM Expression : Pr(married), predict() dy/dx w.r.t. : percentile ---------------------------------------------------------------------------- | Delta-method | dy/dx Std. Err. z P>|z| [95% Conf. Interval] -------------+-------------------------------------------------------------- percentile | .9873168 .062047 15.91 0.000 .8657069 1.108927 ----------------------------------------------------------------------------

Average marginal effect of an interaction between a continuous and a discrete variable as a double derivative:

. margins, expression(normal(_b[education]*education + > _b[percentile]*percentile + > _b[c.education#c.percentile]*education*percentile + > _b[1.divorce#c.education]*education + > _b[1.pdivorce]*pdivorce + _b[1.divorce] + > _b[1.pdivorce#1.divorce]*1.pdivorce + > _b[_cons]) - > normal(_b[education]*education + > _b[percentile]*percentile + > _b[c.education#c.percentile]*education*percentile > + _b[1.pdivorce]*pdivorce + > _b[_cons])) dydx(education) Warning: expression() does not contain predict() or xb(). Average marginal effects Number of obs = 5,000 Model VCE : OIM Expression : normal(_b[education]*education + _b[percentile]*percentile + _b[c.education#c.percentile]*education*percentile + _b[1.divorce#c.education]*education + _b[1.pdivorce]*pdivorce + _b[1.divorce] + _b[1.pdivorce#1.divorce]*1.pdivorce + _b[_cons]) - normal(_b[education]*education + _b[percentile]*percentile + _b[c.education#c.percentile]*education*percentile + _b[1.pdivorce]*pdivorce + _b[_cons]) dy/dx w.r.t. : education ---------------------------------------------------------------------------- | Delta-method | dy/dx Std. Err. z P>|z| [95% Conf. Interval] -------------+-------------------------------------------------------------- education | .0012968 .0061667 0.21 0.833 -.0107897 .0133833 ----------------------------------------------------------------------------

Average marginal effect of an interaction between two continuous variables as a double derivative:

. margins, expression(normalden(xb())*(-xb())*(_b[education] > + percentile*_b[c.education#c.percentile] + > 1.divorce*_b[c.education#1.divorce] + > 2.divorce*_b[c.education#2.divorce])*(_b[percentile] + > education*_b[c.education#c.percentile]) + > normalden(xb())*(_b[c.education#c.percentile])) Predictive margins Number of obs = 5,000 Model VCE : OIM Expression : normalden(xb())*(-xb())*(_b[education] + percentile*_b[c.education#c.percentile] + 1.divorce*_b[c.education#1.divorce] + 2.divorce*_b[c.education#2.divorce])*(_b[percentile] + education*_b[c.education#c.percentile]) + normalden(xb())*(_b[c.education#c.percentile]) ---------------------------------------------------------------------------- | Delta-method | Margin Std. Err. z P>|z| [95% Conf. Interval] -------------+-------------------------------------------------------------- _cons | .0619111 .021641 2.86 0.004 .0194955 .1043267 ----------------------------------------------------------------------------

The simulated data:

clear set obs 5000 set seed 111 generate education = int(rbeta(4,2)*15) generate percentile = int(rbeta(1,7)*100)/100 generate divorce = int(rbeta(1,4)*3) generate pdivorce = runiform()<.6 generate e = rnormal() generate xbeta = .35*(education*.06 + .5*percentile + /// .8*percentile*education /// + .07*education*divorce - .5*divorce - /// .2*pdivorce - divorce*pdivorce - .1) generate married = xbeta + e > 0