## Exploring results of nonparametric regression models

In his blog post, Enrique Pinzon discussed how to perform regression when we don’t want to make any assumptions about functional form—use the **npregress** command. He concluded by asking and answering a few questions about the results using the **margins** and **marginsplot** commands.

Recently, I have been thinking about all the different types of questions that we could answer using **margins** after nonparametric regression, or really after any type of regression. **margins** and **marginsplot** are powerful tools for exploring the results of a model and drawing many kinds of inferences. In this post, I will show you how to ask and answer very specific questions and how to explore the entire response surface based on the results of your nonparametric regression.

The dataset that we will use includes three covariates—continuous variables **x1** and **x2** and categorical variable **a** with three levels. If you would like to follow along, you can use these data by typing **use http://www.stata.com/users/kmacdonald/blog/npblog**.

Let’s first fit our model.

. npregress kernel y x1 x2 i.a, vce(boot, rep(10) seed(111)) (running npregress on estimation sample) Bootstrap replications (10) ----+--- 1 ---+--- 2 ---+--- 3 ---+--- 4 ---+--- 5 .......... Bandwidth ------------------------------------ | Mean Effect -------------+---------------------- x1 | .4795422 1.150255 x2 | .7305017 1.75222 a | .3872428 .3872428 ------------------------------------ Local-linear regression Number of obs = 1,497 Continuous kernel : epanechnikov E(Kernel obs) = 524 Discrete kernel : liracine R-squared = 0.8762 Bandwidth : cross validation ------------------------------------------------------------------------------ | Observed Bootstrap Percentile y | Estimate Std. Err. z P>|z| [95% Conf. Interval] -------------+---------------------------------------------------------------- Mean | y | 15.64238 .3282016 47.66 0.000 15.17505 16.38781 -------------+---------------------------------------------------------------- Effect | x1 | 3.472485 .2614388 13.28 0.000 2.895202 3.846091 x2 | 4.143544 .2003619 20.68 0.000 3.830656 4.441306 | a | (2 vs 1) | 1.512241 .2082431 7.26 0.000 .9880627 1.717007 (3 vs 1) | -10.23723 .2679278 -38.21 0.000 -10.6352 -9.715858 ------------------------------------------------------------------------------ Note: Effect estimates are averages of derivatives for continuous covariates and averages of contrasts for factor covariates.

Enrique discussed interpretation of these results in his blog, so I will not focus on it here. I should, however, point out that because my goal is to show you how to use **margins**, I used only 10 bootstrap replications (a ridiculously small number) when estimating standard errors. In real research, you would certainly want to use more replications both with the **npregress** command and the **margins** commands that follow.

The **npregress** output includes estimates of the effects of **x1**, **x2**, and the levels of **a** on our outcome, but these estimates are probably not enough to answer some of the important questions that we want to address in our research.

Below, I will first show you how you can explore the nonlinear response surface—the expected value of **y** at different combinations of **x1**, **x2**, and **a**. For instance, suppose your outcome variable is the response to a drug, and you want to know the expected value for a female whose weight is 150 pounds and whose cholesterol level is 220 milligrams per deciliter. What about for a male with the same characteristics? How do these expectations change across a range of weights and cholesterol levels?

I will also demonstrate how to answer questions about population averages, counterfactuals, treatment effects, and more. These are exactly the types of questions that policy makers ask. How, on average, does a variable affect the population they are interested in? For instance, suppose your outcome variable is income of individuals in their 20s. What is the expected value of income for this group, the population average? What is the expected value if, instead of having their observed education level, they were all high school graduates? What if they were all college graduates? What is the difference in these values—the effect of college education?

These are just a few examples of the types of questions that you could answer. I will continue with variable names **x1**, **x2**, and **a** below, but you can imagine relevant questions for your research.

**Exploring the response surface**

Let’s start at the beginning. We might want to know the expected value of the outcome at one specific point. To get the expected value of **y** when **a**=1, **x1**=2, and **x2**=5, we can type

. margins 1.a, at(x1=2 x2=5) reps(10) (running margins on estimation sample) Bootstrap replications (10) ----+--- 1 ---+--- 2 ---+--- 3 ---+--- 4 ---+--- 5 .......... Adjusted predictions Number of obs = 1,497 Replications = 10 Expression : mean function, predict() at : x1 = 2 x2 = 5 ------------------------------------------------------------------------------ | Observed Bootstrap Percentile | Margin Std. Err. z P>|z| [95% Conf. Interval] -------------+---------------------------------------------------------------- 1.a | 12.66575 .4428049 28.60 0.000 11.57809 13.0897 ------------------------------------------------------------------------------

We predict that **y**=12.7 at this point.

We could evaluate it at another point, say, **a**=2, **x1**=2, and **x2**=5.

. margins 2.a, at(x1=2 x2=5) reps(10) (running margins on estimation sample) Bootstrap replications (10) ----+--- 1 ---+--- 2 ---+--- 3 ---+--- 4 ---+--- 5 .......... Adjusted predictions Number of obs = 1,497 Replications = 10 Expression : mean function, predict() at : x1 = 2 x2 = 5 ------------------------------------------------------------------------------ | Observed Bootstrap Percentile | Margin Std. Err. z P>|z| [95% Conf. Interval] -------------+---------------------------------------------------------------- 2.a | 14.80935 .6774548 21.86 0.000 13.48154 15.98991 ------------------------------------------------------------------------------

With **a**=2, the expected value of **y** is now 14.8.

What if our interest is in the effect of going from **a**=1 to **a**=2 when **x1**=2 and **x2**=5? That is just a contrast—the difference in our two previous results. Using the **r.** contrast operator with **margins**, we can perform a hypothesis test of whether these two values are the same.

. margins r(1 2).a, at(x1=2 x2=5) reps(10) (running margins on estimation sample) Bootstrap replications (10) ----+--- 1 ---+--- 2 ---+--- 3 ---+--- 4 ---+--- 5 .......... Contrasts of predictive margins Number of obs = 1,497 Replications = 10 Expression : mean function, predict() at : x1 = 2 x2 = 5 ------------------------------------------------ | df chi2 P>chi2 -------------+---------------------------------- a | 1 27.58 0.0000 ------------------------------------------------ -------------------------------------------------------------- | Observed Bootstrap Percentile | Contrast Std. Err. [95% Conf. Interval] -------------+------------------------------------------------ a | (2 vs 1) | 2.143598 .4081472 1.335078 2.523469 --------------------------------------------------------------

The confidence interval for the difference does not include 0. Using a 5% significance level, we find that the expected value is significantly different for these two points of interest.

But we might be interested in more than these two points. Let’s continue holding **x2** to 5 and look at a range of values for **x1**. And let’s estimate the expected values across all three levels of **a**. In other words, let’s look at a slice of the 3-dimensional response surface (at **x2**=5) and examine the relationship between the other two variables.

. margins a, at(x1=(1(1)4) x2=5) reps(10) (running margins on estimation sample) Bootstrap replications (10) ----+--- 1 ---+--- 2 ---+--- 3 ---+--- 4 ---+--- 5 .......... Adjusted predictions Number of obs = 1,497 Replications = 10 Expression : mean function, predict() 1._at : x1 = 1 x2 = 5 2._at : x1 = 2 x2 = 5 3._at : x1 = 3 x2 = 5 4._at : x1 = 4 x2 = 5 ------------------------------------------------------------------------------ | Observed Bootstrap Percentile | Margin Std. Err. z P>|z| [95% Conf. Interval] -------------+---------------------------------------------------------------- _at#a | 1 1 | 9.075231 2.017476 4.50 0.000 5.984868 12.61052 1 2 | 14.35551 2.488854 5.77 0.000 9.944582 17.73349 1 3 | -.2925357 1.499151 -0.20 0.845 -1.991074 2.248606 2 1 | 12.66575 .3285203 38.55 0.000 11.93955 13.03459 2 2 | 14.80935 .8713405 17.00 0.000 12.92743 15.75046 2 3 | 1.82419 .5384259 3.39 0.001 .6331934 2.422084 3 1 | 17.1207 .2446792 69.97 0.000 16.89714 17.61579 3 2 | 17.28286 .3536435 48.87 0.000 16.78102 17.89018 3 3 | 6.2364 .2689665 23.19 0.000 5.695041 6.5686 4 1 | 21.29851 .1441116 147.79 0.000 21.05866 21.54149 4 2 | 18.96943 .2000661 94.82 0.000 18.50348 19.24232 4 3 | 11.45945 .2037134 56.25 0.000 11.05172 11.69403 ------------------------------------------------------------------------------

Better yet, let’s plot these values.

. marginsplot

We find that when **x2**=5, the expected value of **y** increases as **x1** increases, and the expected value is lower for **a**=3 than for **a**=1 and **a**=2 across all levels of **x1**.

But is this pattern the same for other values of **x2**?

We have only three covariates. So we can explore the entire response surface easily. Let’s look at additional slices at other values of **x2**. Here’s the command:

. margins a, at(x1=(1(1)4) x2=(2 5 8)) reps(10)

This produces a lot of output, so I won’t show it. But here is the graph

. marginsplot, bydimension(x2) byopts(cols(3))

We can now see that the response surface changes as **x2** changes. When **x2**=2, the expected value of **y** increases slightly as **x1** increases, but there is almost no difference among the levels of **a**. For **x2**=8, the differences among levels of **a** are more pronounced and appear to have a different pattern, increasing with **x1** and then beginning to level off.

Earlier, we typed **r(1 2).a** to test for a difference in expected values when **a**=1 and **a**=2. Similarly, we could type **r(1 3).a** to compare **a**=1 with **a**=3. We could make both of those comparisons by simply typing **r.a**. And we can do this across a range of values of **x1** and **x2**. We just change **a** to **r.a** in our previous **margins** command.

. margins r.a, at(x1=(1(1)4) x2=(2 5 8)) reps(10) (running margins on estimation sample) Bootstrap replications (10) ----+--- 1 ---+--- 2 ---+--- 3 ---+--- 4 ---+--- 5 .......... Contrasts of predictive margins Number of obs = 1,497 Replications = 10 Expression : mean function, predict() 1._at : x1 = 1 x2 = 2 2._at : x1 = 1 x2 = 5 (output omitted) -------------------------------------------------------------- | Observed Bootstrap Percentile | Contrast Std. Err. [95% Conf. Interval] -------------+------------------------------------------------ a@_at | (2 vs 1) 1 | -.4316907 .2926308 -.6839362 .2414861 (2 vs 1) 2 | 5.280281 .7986439 3.821515 6.316996 (2 vs 1) 3 | 8.716763 9.371615 -7.472455 19.98571 (2 vs 1) 4 | -1.388968 .1918396 -1.756544 -1.110874 (2 vs 1) 5 | 2.143598 .4487908 1.075932 2.563385 (2 vs 1) 6 | 12.88835 4.544753 6.060425 20.50299 (2 vs 1) 7 | -2.384115 .2101401 -2.736413 -2.082235 (2 vs 1) 8 | .1621524 .3430663 -.3141147 .6630604 (2 vs 1) 9 | 9.285226 1.292208 7.300806 10.96053 (2 vs 1) 10 | -2.355113 .5453254 -2.893817 -1.346538 (2 vs 1) 11 | -2.329082 .2355506 -2.717465 -1.844513 (2 vs 1) 12 | 7.264358 1.152335 5.603027 9.361632 (3 vs 1) 1 | -3.622857 .6110854 -4.141618 -2.277811 (3 vs 1) 2 | -9.367766 .9618594 -10.86855 -8.056798 (3 vs 1) 3 | -4.751952 8.78832 -14.17876 12.3299 (3 vs 1) 4 | -2.500585 .2942087 -2.999642 -2.113881 (3 vs 1) 5 | -10.84156 .3575732 -11.46529 -10.4761 (3 vs 1) 6 | -18.84966 2.783479 -21.87325 -11.56521 (3 vs 1) 7 | -.3080399 .1685899 -.7749465 -.1793496 (3 vs 1) 8 | -10.8843 .271959 -11.28382 -10.48519 (3 vs 1) 9 | -22.70852 1.147957 -23.26496 -19.51463 (3 vs 1) 10 | 3.922399 .6064593 2.976073 4.793276 (3 vs 1) 11 | -9.839058 .2572868 -10.13578 -9.339336 (3 vs 1) 12 | -20.31851 1.294548 -22.1466 -17.29554 --------------------------------------------------------------

The legend at the top of the output tells us that **1._at** corresponds to **x1**=1 and **x2**=2. The values parentheses, such as **(2 vs 1)**, at the first of each line in the table tell us which values of **a** are compared in that line. Thus, the first line in the table provides a test comparing expected values of **y** for **a**=2 versus **a**=1 when **x1**=1 and **x2**=2. Admittedly, this is a lot to look at, and it is probably easier to interpret with a graph. We use **marginsplot** to graph these differences with their confidence intervals. This time, let’s use the **yline(0)** option to add a reference line at 0. This allows us to perform the test visually by checking whether the confidence interval for the difference includes 0.

. marginsplot, bydimension(x2) byopts(cols(3)) yline(0)

In this case, some of the confidence intervals are so narrow that they are hard to see. If we look closely at the blue point on the far left, we see that the confidence interval for the difference comparing **a**=2 versus **a**=1 when **x1**=1 and **x2**=2 (corresponding to the first line in the output above) does include 0. This indicates that there is not a significant difference in these expected values. We can examine each of the other points and confidence intervals in the same way. For example, looking at the red line and points in the third panel, we see that the effect of moving from **a**=1 to **a**=3 is negative and significantly different from 0 for **x1** values of 2, 3, and 4. When **x1** is 1, the point estimate of the effect is still negative, but that effect is not significantly different from 0 at the 95% level. But recall that we should dramatically increase the number of bootstrap replicates to make any real claims about confidence intervals.

So far, we have compared both **a**=2 with **a**=1 and **a**=3 with **a**=1. But we are not limited to making comparisons with **a**=1. We could compare 1 with 2 and 2 with 3, which often makes more sense if levels of **a** have a natural ordering. To do this, we just replace **r.** with **ar.** in our **margins** command. I won’t show that output, but you have the data and can try it if you like.

**Population-averaged results**

So far, we have talked about evaluating individual points on your response surface and how to perform tests comparing expected values at those points. Now, let’s switch gears and talk about population-averaged results.

We are going to need the dataset to be representative of the population. If that is not true of your data, you will want to stop with the analyses we did above. We will assume our data are representative so that we can answer a variety of questions based on average predictions.

First, what is the overall expected population mean from this response surface?

. margins, reps(10) (running margins on estimation sample) Bootstrap replications (10) ----+--- 1 ---+--- 2 ---+--- 3 ---+--- 4 ---+--- 5 .......... Predictive margins Number of obs = 1,497 Replications = 10 Expression : mean function, predict() ------------------------------------------------------------------------------ | Observed Bootstrap Percentile | Margin Std. Err. z P>|z| [95% Conf. Interval] -------------+---------------------------------------------------------------- _cons | 15.64238 .318746 49.07 0.000 15.34973 16.24543 ------------------------------------------------------------------------------

Whatever the process was that generated this, we believe 15.6 is the expected value in the population and [15.3, 16.2] is the confidence intervals for it.

Do the population averages differ when we first set everyone to have **a**=1, then set everyone to have **a**=2, and finally set everyone to have **a**=3? Let’s look at the expected means for all three.

. margins a, reps(10) (running margins on estimation sample) Bootstrap replications (10) ----+--- 1 ---+--- 2 ---+--- 3 ---+--- 4 ---+--- 5 .......... Predictive margins Number of obs = 1,497 Replications = 10 Expression : mean function, predict() ------------------------------------------------------------------------------ | Observed Bootstrap Percentile | Margin Std. Err. z P>|z| [95% Conf. Interval] -------------+---------------------------------------------------------------- a | 1 | 18.39275 .4847808 37.94 0.000 17.63322 19.31199 2 | 19.90499 .6840056 29.10 0.000 18.17305 20.57171 3 | 8.155515 .3239624 25.17 0.000 7.59187 8.70898 ------------------------------------------------------------------------------

We get 18.4, 19.9, and 8.2. They sure don’t look the same. Let’s test them.

. margins r.a, reps(10) (running margins on estimation sample) Bootstrap replications (10) ----+--- 1 ---+--- 2 ---+--- 3 ---+--- 4 ---+--- 5 .......... Contrasts of predictive margins Number of obs = 1,497 Replications = 10 Expression : mean function, predict() ------------------------------------------------ | df chi2 P>chi2 -------------+---------------------------------- a | (2 vs 1) | 1 17.10 0.0000 (3 vs 1) | 1 618.86 0.0000 Joint | 2 1337.66 0.0000 ------------------------------------------------ -------------------------------------------------------------- | Observed Bootstrap Percentile | Contrast Std. Err. [95% Conf. Interval] -------------+------------------------------------------------ a | (2 vs 1) | 1.512241 .3656855 1.051265 2.016566 (3 vs 1) | -10.23723 .4115157 -10.82644 -9.4265 --------------------------------------------------------------

In the causal-inference or treatment-effects literature, the means would be considered potential-outcome means, and these differences would be the average treatment effects of a multivalued treatment. Here the average treatment effect of **a**=2 (compared with **a**=1) is 1.5.

We saw in the previous section that differences in expected values for levels of **a** varied across values of **x2**. Let’s estimate the potential-outcome means and treatment effects of **a** at different values of **x2**. Notice that these are still population averages because, unlike in the previous section, we are not setting **x1** to any specific value. Instead, predictions use the observed values of **x1** in the data.

. margins a, at(x2=(2 5 8)) reps(10) (running margins on estimation sample) Bootstrap replications (10) ----+--- 1 ---+--- 2 ---+--- 3 ---+--- 4 ---+--- 5 .......... Predictive margins Number of obs = 1,497 Replications = 10 Expression : mean function, predict() 1._at : x2 = 2 2._at : x2 = 5 3._at : x2 = 8 ------------------------------------------------------------------------------ | Observed Bootstrap Percentile | Margin Std. Err. z P>|z| [95% Conf. Interval] -------------+---------------------------------------------------------------- _at#a | 1 1 | 3.707771 .6552861 5.66 0.000 4.619393 6.70335 1 2 | 1.340815 .9491373 1.41 0.158 2.077606 5.010682 1 3 | 3.677883 .935561 3.93 0.000 4.33317 7.380818 2 1 | 17.16092 .27028 63.49 0.000 16.8047 17.75164 2 2 | 17.22508 .1993421 86.41 0.000 16.91111 17.57034 2 3 | 6.933778 .1929387 35.94 0.000 6.751693 7.332845 3 1 | 30.61528 .9961323 30.73 0.000 28.66033 31.88367 3 2 | 40.17685 1.687315 23.81 0.000 36.18359 41.73811 3 3 | 10.97667 .9771489 11.23 0.000 9.510952 12.59133 ------------------------------------------------------------------------------

Instead of looking at the output, let’s graph these potential-outcome means.

. marginsplot

The effect of **a** increases as **x2** increases. The effect is largest when **x2**=8.

Now, we can test for differences at each level of **x2**

. margins r.a, at(x2=(2 5 8)) reps(10) (running margins on estimation sample) Bootstrap replications (10) ----+--- 1 ---+--- 2 ---+--- 3 ---+--- 4 ---+--- 5 .......... Contrasts of predictive margins Number of obs = 1,497 Replications = 10 Expression : mean function, predict() 1._at : x2 = 2 2._at : x2 = 5 3._at : x2 = 8 ------------------------------------------------ | df chi2 P>chi2 -------------+---------------------------------- a@_at | (2 vs 1) 1 | 1 51.80 0.0000 (2 vs 1) 2 | 1 0.13 0.7184 (2 vs 1) 3 | 1 43.16 0.0000 (3 vs 1) 1 | 1 0.01 0.9322 (3 vs 1) 2 | 1 834.54 0.0000 (3 vs 1) 3 | 1 245.45 0.0000 Joint | 6 2732.16 0.0000 ------------------------------------------------ -------------------------------------------------------------- | Observed Bootstrap Percentile | Contrast Std. Err. [95% Conf. Interval] -------------+------------------------------------------------ a@_at | (2 vs 1) 1 | -2.366956 .3288861 -2.737075 -1.466043 (2 vs 1) 2 | .0641596 .1779406 -.224935 .3180577 (2 vs 1) 3 | 9.561575 1.455483 6.847319 11.96769 (3 vs 1) 1 | -.0298883 .3513035 -.6345751 .4584875 (3 vs 1) 2 | -10.22714 .3540221 -10.89414 -9.808236 (3 vs 1) 3 | -19.63861 1.253521 -21.45499 -18.0223 --------------------------------------------------------------

Again, let’s look at the graph.

. marginsplot

The difference in means when **a**=3 and **a**=1, the treatment effect, is not significant when **x2**=2. Neither is the effect of **a**=2 versus **a**=1 when **x2**=5. All the other effects are significantly different from 0.

**Conclusion**

In this blog, we have explored the response surface of a nonlinear function, we have estimated a variety of population averages based on our nonparametric regression model, and we have performed several tests comparing expected values at specific points of the response surface and tests comparing population averages. Yet, we have only scratched the surface of the types of estimates and tests you can obtain using **margins** after **npregress**. There are additional contrasts operators that will allow you to test differences from a grand mean, differences from means of previous or subsequent levels, and more. See [R] contrast for details of the available contrast operators. You can also use **marginsplot** to see the results of **margins** commands from different angles. For instance, if we type **marginsplot, bydimension(x1)** instead of **marginsplot, bydimension(x2)**, we see our nonlinear response surface from a different perspective. See [R] marginsplot for details and examples of this command.

Whether you use nonparametric regression or another model, **margins** and **marginsplot** are the solution for exploring the results, making inferences, and understanding relationships among the variables you are studying.