Multiple regression models often contain interaction terms. This FAQ page covers the situation in which there are two moderator variables which jointly influence the regression of the dependent variable on an independent variable. In other words, a regression model that has a significant three-way interaction of continuous variables.
If you are new to the topic of interactions of continuous variables, you may want to start with the FAQ page the discusses the simpler continuous by continuous two-way interaction, How can I explain a continuous by continuous interaction?
Back to the three-way interaction. There are several approaches that one might use to explain a three-way interaction. The approach that we will demonstrate is to compute simple slopes, i.e., the slopes of the dependent variable on the independent variable when the moderator variables are held constant at different combinations of high and low values. After computing the simple slopes, we will then compute and test the differences among all pairs of the slopes. The method used in this FAQ is adapted from a paper by Dawson and Richter (2004). The Dawson and Richter method computes the differences in simple slopes and their standard errors using a number of moderately-complex formulas. We will take an easier path and let margins do all of the heavy lifting.
We will begin by looking at the regression equation which includes a three-way continuous interaction. In the formula, Y is the response variable, X the predictor (independent) variable with Z and W being the two moderator variables.
Y = b0 + b1X + b2Z + b3W + b4XZ + b5XW + b6ZW + b7XZW
We can reorder the terms into two groups, the first grouping (terms that do not contain X) defines the intercept while the second grouping (all the terms that contain an X) defines the simple slope.
Y = (b0 + b2Z + b3W + b6ZW) + (b1 + b4Z + b5W + b7ZW)X
Next, we will define high values of Z and W as being one standard deviation above their respective means and will denote them as Hz and Hw. The low values will be defined as one standard deviation below their means and will be as Lz and Lw. Altogether there are four possible combinations of conditions: 1) HzHw, 2) HzLw, 3) LzHw and 4) LzLw.
Here is the formula for the simple slope for condition 1) when both Z and W are at their high values. The formulas for the other three conditions are exactly the same except that the values for Z and W are different.
simple slope 1 = b1 + b4Hz + b5Hw + b7HzHw
We will illustrate the simple slopes process using the dataset hsb2 that has a statistically significant (barely) three-way continuous interaction. In order to keep the notation consistent we will temporarily change the names of the variables; y for the response variable, x for the independent variable, z for the first moderator, and w for the second moderator. It is clear from the code below that write is the response variable, read is the IV and math and science are the moderator variables.
use https://stats.idre.ucla.edu/stat/stata/notes/hsb2, clear
rename write y
rename read x
rename math w
rename science z
regress y c.x##c.z##c.w
Source | SS df MS Number of obs = 200
-------------+------------------------------ F( 7, 192) = 26.28
Model | 8748.54751 7 1249.7925 Prob > F = 0.0000
Residual | 9130.32749 192 47.553789 R-squared = 0.4893
-------------+------------------------------ Adj R-squared = 0.4707
Total | 17878.875 199 89.843593 Root MSE = 6.8959
------------------------------------------------------------------------------
y | Coef. Std. Err. t P>|t| [95% Conf. Interval]
-------------+----------------------------------------------------------------
x | -3.013849 1.802977 -1.67 0.096 -6.570035 .5423358
z | -2.435316 1.757363 -1.39 0.167 -5.901532 1.0309
|
c.x#c.z | .0553139 .0325207 1.70 0.091 -.0088297 .1194576
|
w | -3.228097 1.926393 -1.68 0.095 -7.027708 .5715142
|
c.x#c.w | .0737523 .035898 2.05 0.041 .0029472 .1445575
|
c.z#c.w | .0611136 .035031 1.74 0.083 -.0079813 .1302086
|
c.x#c.z#c.w | -.0012563 .0006277 -2.00 0.047 -.0024944 -.0000182
|
_cons | 166.5438 93.12788 1.79 0.075 -17.1413 350.2289
------------------------------------------------------------------------------
Please note that the three-way interaction, xzw, is statistically significant with a p-value of .047.
We will create four global macros that contain the values w and z when that are one standard deviation above the mean and one standard deviation below the mean.
/* determine range of x for graphing */
sum x
Variable | Obs Mean Std. Dev. Min Max
-------------+--------------------------------------------------------
x | 200 52.23 10.25294 28 76
sum w
Variable | Obs Mean Std. Dev. Min Max
-------------+--------------------------------------------------------
w | 200 52.645 9.368448 33 75
/* determined high and low values for w */
global Hw=r(mean)+r(sd)
global Lw=r(mean)-r(sd)
sum z
Variable | Obs Mean Std. Dev. Min Max
-------------+--------------------------------------------------------
z | 200 51.85 9.900891 26 74
/* determine high and low values for z */
global Hz=r(mean)+r(sd)
global Lz=r(mean)-r(sd)
We will use the margins comand to compute the simple slopes for the four values of w and z. We use the dydx option to obtain the slopes. Please note, we have manually annotated the output with the various combinations of w and z in parentheses.
margins, dydx(x) at(w=($Hw $Lw) z=($Hz $Lz)) vsquish
Average marginal effects Number of obs = 200
Model VCE : OLS
Expression : Linear prediction, predict()
dy/dx w.r.t. : x
1._at : z = 61.75089
w = 62.01345
2._at : z = 61.75089
w = 43.27655
3._at : z = 41.94911
w = 62.01345
4._at : z = 41.94911
w = 43.27655
------------------------------------------------------------------------------
| Delta-method
| dy/dx Std. Err. z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
x |
_at |
(HwHz) 1 | .164556 .0953841 1.73 0.084 -.0223935 .3515054
(LwHz) 2 | .236248 .1429934 1.65 0.099 -.0440138 .5165098
(HwLz) 3 | .6119673 .1962032 3.12 0.002 .227416 .9965186
(LwLz) 4 | .2175364 .1275474 1.71 0.088 -.0324519 .4675248
------------------------------------------------------------------------------
The simple slopes for HwHz, LwHz and LwLz are all fairly similar. However, the slope for Hwlz seems to be rather different from the other three.
The next step is to graph the slopes along with a scatterplot of the observations. We will use the margins command again to get the information for plotting. Because we are using linear regression, we need only two points to define each simple slope regression line. We will use the x values 20 and 80 which cover more the full range of variable x.
(introduced in Stata 12) command.
margins, at(x=(20 80) w=($Hw $Lw) z=($Hz $Lz)) vsquish
Adjusted predictions Number of obs = 200
Model VCE : OLS
Expression : Linear prediction, predict()
1._at : x = 20
z = 61.75089
w = 62.01345
2._at : x = 20
z = 61.75089
w = 43.27655
3._at : x = 20
z = 41.94911
w = 62.01345
4._at : x = 20
z = 41.94911
w = 43.27655
5._at : x = 80
z = 61.75089
w = 62.01345
6._at : x = 80
z = 61.75089
w = 43.27655
7._at : x = 80
z = 41.94911
w = 62.01345
8._at : x = 80
z = 41.94911
w = 43.27655
------------------------------------------------------------------------------
| Delta-method
| Margin Std. Err. z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
_at |
1 | 53.29426 4.059313 13.13 0.000 45.33815 61.25037
2 | 44.50291 5.033491 8.84 0.000 34.63745 54.36837
3 | 35.41995 6.802609 5.21 0.000 22.08708 48.75282
4 | 39.98077 3.033838 13.18 0.000 34.03455 45.92698
5 | 63.16762 1.926701 32.79 0.000 59.39135 66.94388
6 | 58.67779 4.006649 14.65 0.000 50.8249 66.53068
7 | 72.13799 5.488377 13.14 0.000 61.38097 82.89501
8 | 53.03295 4.786994 11.08 0.000 43.65062 62.41529
------------------------------------------------------------------------------
Next, we use the marginsplot command (introduced in Stata 12) to do the actual plotting. The addplot option allows us to include a scatterplot of the observations along with the simple slope regression lines.
marginsplot, recast(line) noci addplot(scatter y x, msym(oh) jitter(3))
Now we have a good idea of what is going on with these data.
Hey, wait a minute, my model also has covariates. What about them? We will deal with adding covariates to the model in the next section.
Three-way interaction of continuous variables with covariates
Here is the regression model from the previous page with two covariates V1 and V2 added.
Y = b0 + b1X + b2Z + b3W + b4XZ + b5XW + b6ZW + b7XZW + b8V1 + b9V2
Adding covariates to your model does not change the formulas for the simple slopes as long as the covariates are not interacted with the independent variable or the moderator variables.
We will add two covariates to our previous model, socst which will be renamed as v1 and ses which will be renamed v2. Everything is this section is the same as the previous example with the exception of adding the covariates to the end of the regress command.
rename socst v1
rename ses v2
regress y c.x##c.z##c.w v1 v2
Source | SS df MS Number of obs = 200
-------------+------------------------------ F( 9, 190) = 25.36
Model | 9757.27668 9 1084.14185 Prob > F = 0.0000
Residual | 8121.59832 190 42.7452543 R-squared = 0.5457
-------------+------------------------------ Adj R-squared = 0.5242
Total | 17878.875 199 89.843593 Root MSE = 6.538
------------------------------------------------------------------------------
y | Coef. Std. Err. t P>|t| [95% Conf. Interval]
-------------+----------------------------------------------------------------
x | -3.463037 1.716525 -2.02 0.045 -6.848932 -.0771426
z | -2.902201 1.670106 -1.74 0.084 -6.196533 .3921302
|
c.x#c.z | .0624411 .03093 2.02 0.045 .0014308 .1234514
|
w | -3.696588 1.831366 -2.02 0.045 -7.309009 -.0841671
|
c.x#c.w | .0796921 .0341529 2.33 0.021 .0123245 .1470597
|
c.z#c.w | .0696334 .0332851 2.09 0.038 .0039775 .1352893
|
c.x#c.z#c.w | -.0013838 .0005968 -2.32 0.021 -.0025611 -.0002065
|
v1 | .2789791 .057697 4.84 0.000 .1651703 .392788
v2 | -.9056562 .6914623 -1.31 0.192 -2.269585 .4582727
_cons | 185.3134 88.56438 2.09 0.038 10.6177 360.0092
------------------------------------------------------------------------------
Note that with the inclusion the covariates that the three-way interaction is even “more” significant than before.
Next, we will compute the simple slopes for the various combinations of w and z. The margins command is exactly the same as for the example without the covariates because the covariances are not involved in the 3-way interaction.
/* compute simple slopes */
margins, dydx(x) at(w=($Hw $Lw) z=($Hz $Lz)) vsquish
Average marginal effects Number of obs = 200
Model VCE : OLS
Expression : Linear prediction, predict()
dy/dx w.r.t. : x
1._at : z = 61.75089
w = 62.01345
2._at : z = 61.75089
w = 43.27655
3._at : z = 41.94911
w = 62.01345
4._at : z = 41.94911
w = 43.27655
------------------------------------------------------------------------------
| Delta-method
| dy/dx Std. Err. z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
x |
_at |
1 | .0356648 .0951734 0.37 0.708 -.1508717 .2222012
2 | .1435565 .1371306 1.05 0.295 -.1252146 .4123276
3 | .498484 .1884961 2.64 0.008 .1290384 .8679296
4 | .0929561 .1236172 0.75 0.452 -.1493291 .3352413
------------------------------------------------------------------------------
Let’s plot the four simple slopes on top of the observations. We begin the graphing by computing expected values using the margins command for different values of ww and z. Again, we need only two values of x (20 & 80) to define each simple slope regression line.
margins, at(x=(20 80) w=($Hw $Lw) z=($Hz $Lz)) vsquish
Predictive margins Number of obs = 200
Model VCE : OLS
Expression : Linear prediction, predict()
1._at : x = 20
z = 61.75089
w = 62.01345
2._at : x = 20
z = 61.75089
w = 43.27655
3._at : x = 20
z = 41.94911
w = 62.01345
4._at : x = 20
z = 41.94911
w = 43.27655
5._at : x = 80
z = 61.75089
w = 62.01345
6._at : x = 80
z = 61.75089
w = 43.27655
7._at : x = 80
z = 41.94911
w = 62.01345
8._at : x = 80
z = 41.94911
w = 43.27655
------------------------------------------------------------------------------
| Delta-method
| Margin Std. Err. z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
_at |
1 | 56.98694 3.936984 14.47 0.000 49.27059 64.70329
2 | 47.84013 4.825798 9.91 0.000 38.38174 57.29852
3 | 38.20389 6.500797 5.88 0.000 25.46256 50.94522
4 | 44.62434 3.040215 14.68 0.000 38.66563 50.58306
5 | 59.12683 2.060024 28.70 0.000 55.08925 63.1644
6 | 56.45352 3.835255 14.72 0.000 48.93656 63.97048
7 | 68.11293 5.302337 12.85 0.000 57.72054 78.50532
8 | 50.20171 4.580553 10.96 0.000 41.22399 59.17943
------------------------------------------------------------------------------
marginsplot, recast(line) noci addplot(scatter y x, msym(oh) jitter(3))
With the covariates in the model we have a bit more variation in the simple slopes. However, once again it looks like slope 3 is the one that is the most different from the others.
We can compute tests of differences in simple slopes using the pwcompare(effects) option for margins. This option was introduced in Stata 12.
/* compute pairwise differences in simple slopes */
margins, dydx(x) at(w=($Hw $Lw) z=($Hz $Lz)) vsquish pwcompare(effects)
Pairwise comparisons of average marginal effects
Model VCE : OLS
Expression : Linear prediction, predict()
dy/dx w.r.t. : x
1._at : z = 61.75089
w = 62.01345
2._at : z = 61.75089
w = 43.27655
3._at : z = 41.94911
w = 62.01345
4._at : z = 41.94911
w = 43.27655
------------------------------------------------------------------------------
| Contrast Delta-method Unadjusted Unadjusted
| dy/dx Std. Err. z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
x |
_at |
2 vs 1 | .1078917 .1559956 0.69 0.489 -.197854 .4136375
3 vs 1 | .4628192 .1955057 2.37 0.018 .0796352 .8460033
4 vs 1 | .0572913 .1557822 0.37 0.713 -.2480361 .3626188
3 vs 2 | .3549275 .2449055 1.45 0.147 -.1250785 .8349335
4 vs 2 | -.0506004 .1635762 -0.31 0.757 -.3712038 .2700031
4 vs 3 | -.4055279 .2096843 -1.93 0.053 -.8165015 .0054457
------------------------------------------------------------------------------
Inspecting the table above it appears slope 3 versus slope 1 (LzHw vs HzHw) is the only one of the tests that is significant. But before we reach that conclusion we will need to adjust the p-values to take into account that these are post-hoc tests. To do this, we will add the mcompare(bonferroni) option to the margins command.
/* compute pairwise differences in simple slopes with Bonferroni correction */
margins, dydx(x) at(w=($Hw $Lw) z=($Hz $Lz)) vsquish pwcompare(effects) mcompare(bonferroni)
Pairwise comparisons of average marginal effects
Model VCE : OLS
Expression : Linear prediction, predict()
dy/dx w.r.t. : x
1._at : z = 61.75089
w = 62.01345
2._at : z = 61.75089
w = 43.27655
3._at : z = 41.94911
w = 62.01345
4._at : z = 41.94911
w = 43.27655
---------------------------
| Number of
| Comparisons
-------------+-------------
x |
_at | 6
---------------------------
------------------------------------------------------------------------------
| Contrast Delta-method Bonferroni Bonferroni
| dy/dx Std. Err. z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
x |
_at |
2 vs 1 | .1078917 .1559956 0.69 1.000 -.3036648 .5194483
3 vs 1 | .4628192 .1955057 2.37 0.108 -.052975 .9786135
4 vs 1 | .0572913 .1557822 0.37 1.000 -.3537022 .4682849
3 vs 2 | .3549275 .2449055 1.45 0.884 -.2911963 1.001051
4 vs 2 | -.0506004 .1635762 -0.31 1.000 -.4821564 .3809557
4 vs 3 | -.4055279 .2096843 -1.93 0.319 -.9587289 .1476731
------------------------------------------------------------------------------
With the Bonferroni correction the p-value for 3 vs 1 becomes 0.108, which is not statistically significant. The Bonferroni correction may not be the best choice because it is overly conservative but it is certainly easy to implement.
References
Dawson, J.F. & Richter, A.W. (2004) A significance test of slope differences for three-way interactions in moderated multiple regression analysis. RP0422. Aston University, Birmingham, UK