Home > Statistics > Exploring results of nonparametric regression models

## 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.

Categories: Statistics Tags: