Stata FAQ
Part 1: How can I understand a 3-way continuous interaction? (Stata 10 and earlier)

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 lincom 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 are denoted as Lz and Lw. Altogether there are four possible combinations of conditions: 1) HzHw, 2) HzLw, 3) LzHw and 4) LzLw.

Here are the formulas for the simple slope and intercept 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

intercept 1 =  b0 + b2Hz + b3Hw + b6HzHw
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. After renaming, we will create the interactions manually. We could have used xi3 but chose to do it manually to keep the notation as concise and consistent as possible.

Then, after creating the interaction terms, we will run the regression model. 

use http://www.ats.ucla.edu/stat/stata/notes/hsb2, clear

rename write y
rename read x
rename math w
rename science z

generate xz  = x*z
generate xw  = x*w
generate zw  = z*w
generate xzw = x*z*w

regress y x z w xz xw zw xzw

      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
           w |  -3.228097   1.926393    -1.68   0.095    -7.027708    .5715142
          xz |   .0553139   .0325207     1.70   0.091    -.0088297    .1194576
          xw |   .0737523    .035898     2.05   0.041     .0029472    .1445575
          zw |   .0611136    .035031     1.74   0.083    -.0079813    .1302086
         xzw |  -.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 need to create two global macro variables to hold the range of x that will be used for graphing. We will also 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. 

sum x
global hi=r(max)
global lo=r(min)
sum w
global Hw=r(mean)+r(sd)
global Lw=r(mean)-r(sd)
sum z
global Hz=r(mean)+r(sd)
global Lz=r(mean)-r(sd)
In order to keep the lincom commands as simple as possible we will create additional global macros that contain the code for the four simple slopes. Note that each of the macros have the same terms from the model; x, xz, xw and xzw. They only differ in the four combinations of high and low values for z and w. 
global HzHw "x + ($Hz)*xz + ($Hw)*xw + ($Hz)*($Hw)*xzw" 
global HzLw "x + ($Hz)*xz + ($Lw)*xw + ($Hz)*($Lw)*xzw"
global LzHw "x + ($Lz)*xz + ($Hw)*xw + ($Lz)*($Hw)*xzw"
global LzLw "x + ($Lz)*xz + ($Lw)*xw + ($Lz)*($Lw)*xzw"
Please note that the x, xz, xw and xzw above refer not to the values of the variables themselves but refer to the coefficients b1, b4, b5 and b7 respectively.

Using the global macro variables created in the last step, we will compute each of the four simple slopes. Note that lincom substitutes in the actual numeric values for the high and low values of z and w (61.75, 62.01, 3829.38, etc.). We will save each of the simple slopes in a global macro variable starting with the letter "b." 

/* simple slopes */

lincom $HzHw /* slope 1 */

 ( 1)  x + 61.75089 xz + 62.01345 xw + 3829.386 xzw = 0

------------------------------------------------------------------------------
           y |      Coef.   Std. Err.      t    P>|t|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
         (1) |    .164556   .0953841     1.73   0.086    -.0235794    .3526913
------------------------------------------------------------------------------

global b1 = r(estimate)

lincom $HzLw /* slope 2 */


 ( 1)  x + 61.75089 xz + 43.27655 xw + 2672.366 xzw = 0

------------------------------------------------------------------------------
           y |      Coef.   Std. Err.      t    P>|t|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
         (1) |    .236248   .1429934     1.65   0.100    -.0457916    .5182876
------------------------------------------------------------------------------

global b2 = r(estimate) 

lincom $LzHw /* slope 3 */

 ( 1)  x + 41.94911 xz + 62.01345 xw + 2601.409 xzw = 0

------------------------------------------------------------------------------
           y |      Coef.   Std. Err.      t    P>|t|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
         (1) |   .6119673   .1962032     3.12   0.002     .2249767    .9989579
------------------------------------------------------------------------------

global b3 = r(estimate)

lincom $LzLw /* slope 4 */

 ( 1)  x + 41.94911 xz + 43.27655 xw + 1815.413 xzw = 0

------------------------------------------------------------------------------
           y |      Coef.   Std. Err.      t    P>|t|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
         (1) |   .2175364   .1275474     1.71   0.090    -.0340377    .4691105
------------------------------------------------------------------------------

global b4 = r(estimate)
With some copying and pasting from above, we can put together a table of the four simple slopes. This will make it easier to visually compare the values of the simple slopes with one another. 
               |      Coef.   Std. Err.      t    P>|t|     [95% Conf. Interval]
---------------+----------------------------------------------------------------
(slope 1 HzHw) |   .164556    .0953841     1.73   0.086    -.0235794    .3526913
(slope 2 HzLw) |   .236248    .1429934     1.65   0.100    -.0457916    .5182876
(slope 3 LzHw) |   .6119673   .1962032     3.12   0.002     .2249767    .9989579
(slope 4 LzLw) |   .2175364   .1275474     1.71   0.090    -.0340377    .4691105
If we are going to make a graph of the simple regression lines we will also need to compute the four intercepts (constants) to go along with the slopes. We will save each of the constants in a global macro variable starting with the letter "c." 
/* intercepts */

lincom _cons + $Hz*z + $Hw*w + $Hz*$Hw*zw

 ( 1)  61.75089 z + 62.01345 w + 3829.386 zw + _cons = 0

------------------------------------------------------------------------------
           y |      Coef.   Std. Err.      t    P>|t|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
         (1) |   50.00314   5.940565     8.42   0.000     38.28599    61.72029
------------------------------------------------------------------------------

global c1 = r(estimate)

lincom _cons + $Hz*z + $Lw*w + $Hz*$Lw*zw

 ( 1)  61.75089 z + 43.27655 w + 2672.366 zw + _cons = 0

------------------------------------------------------------------------------
           y |      Coef.   Std. Err.      t    P>|t|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
         (1) |   39.77795   7.819569     5.09   0.000     24.35466    55.20124
------------------------------------------------------------------------------

global c2 = r(estimate)

lincom _cons + $Lz*z + $Hw*w + $Lz*$Hw*zw

 ( 1)  41.94911 z + 62.01345 w + 2601.409 zw + _cons = 0

------------------------------------------------------------------------------
           y |      Coef.   Std. Err.      t    P>|t|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
         (1) |   23.18061   10.64204     2.18   0.031     2.190283    44.17093
------------------------------------------------------------------------------

global c3 = r(estimate)

lincom _cons + $Lz*z + $Lw*w + $Lz*$Lw*zw

 ( 1)  41.94911 z + 43.27655 w + 1815.413 zw + _cons = 0

------------------------------------------------------------------------------
           y |      Coef.   Std. Err.      t    P>|t|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
         (1) |   35.63004   5.537431     6.43   0.000     24.70803    46.55205
------------------------------------------------------------------------------

global c4 = r(estimate)
Now we have the slopes and constants for each of the four regression lines we can plot the four simple regression lines. We do this by holding the values of the moderator variables, z and w, constant at either their high (mean + 1 SD) or low values (mean - 1 SD) and allowing the independent variable, x, to vary.

Please note: You must use the names y and x with the twoway function plot command. Do not use variable names from your data. 

twoway (function y=$c1+$b1*x, range($lo $hi))  ///
       (function y=$c2+$b2*x, range($lo $hi))  ///
       (function y=$c3+$b3*x, range($lo $hi))  ///
       (function y=$c4+$b4*x, range($lo $hi))  ///
       (scatter y x, msym(oh) jitter(3)),     ///
       legend(order(1 "1 HzHw" 2 "2 HzLw" 3 "3 LzHw" 4 "4 LzLw")) /// 
       xtitle(X) ytitle(Y) name(con3way, replace)
Inspecting the graph above, it seems that the biggest differences in slopes is between slopes 1 & 3, 2 & 3, and 4 & 3, i.e., the slope for LzHw (3) seems to be the most different from the others.

It should be noted that the lines on this graph were not plotted from the actual data in the dataset, rather they are theoretical regression lines that have the same slopes and intercepts but are plotted using the entire range of x.

Next, we will compute the tests of differences in simple slopes. There are six possible comparisons among the four slopes. 

/* differences in simple slopes */

lincom ($HzHw)-($HzLw)  /* a) HzHw vs HzLw (1 & 2) */

 ( 1)  18.7369 xw + 1157.02 xzw = 0

------------------------------------------------------------------------------
           y |      Coef.   Std. Err.      t    P>|t|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
         (1) |   -.071692   .1642481    -0.44   0.663    -.3956544    .2522703
------------------------------------------------------------------------------

lincom ($HzHw)-($LzHw)  /* b) HzHw vs LzHw (1 & 3) */

 ( 1)  19.80178 xz + 1227.977 xzw = 0

------------------------------------------------------------------------------
           y |      Coef.   Std. Err.      t    P>|t|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
         (1) |  -.4474113   .2060755    -2.17   0.031    -.8538739   -.0409487
------------------------------------------------------------------------------

lincom ($HzLw)-($LzLw)  /* c) HzLw vs LzLw (2 & 4) */

 ( 1)  19.80178 xz + 856.9528 xzw = 0

------------------------------------------------------------------------------
           y |      Coef.   Std. Err.      t    P>|t|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
         (1) |   .0187116   .1722173     0.11   0.914    -.3209691    .3583924
------------------------------------------------------------------------------

lincom ($LzHw)-($LzLw)  /* d) LzHw vs LzLw (3 & 4) */

 ( 1)  18.7369 xw + 785.9961 xzw = 0

------------------------------------------------------------------------------
           y |      Coef.   Std. Err.      t    P>|t|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
         (1) |   .3944309   .2202323     1.79   0.075    -.0399546    .8288163
------------------------------------------------------------------------------

lincom ($HzHw)-($LzLw)  /* e) HzHw vs LzLw (1 & 4) */

 ( 1)  19.80178 xz + 18.7369 xw + 2013.973 xzw = 0

------------------------------------------------------------------------------
           y |      Coef.   Std. Err.      t    P>|t|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
         (1) |  -.0529804    .163752    -0.32   0.747    -.3759642    .2700034
------------------------------------------------------------------------------

lincom ($HzLw)-($LzHw)  /* f) HzLw vs LzHw (2 & 3) */

 ( 1)  19.80178 xz - 18.7369 xw + 70.95675 xzw = 0

------------------------------------------------------------------------------
           y |      Coef.   Std. Err.      t    P>|t|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
         (1) |  -.3757192   .2579781    -1.46   0.147    -.8845543    .1331159
------------------------------------------------------------------------------
Note how lincom has simplified the equations by combining like terms and canceling whenever possible so that the each lincom test has only two or, at most, three terms.

Again, with some cutting and pasting we can put together a table including each of the tests of differences in simple slopes. 

   Slopes      |      Diff.   Std. Err.      t    P>|t|     [95% Conf. Interval]
---------------+----------------------------------------------------------------
(HzHw vs HzLw) |   -.071692   .1642481    -0.44   0.663    -.3956544    .2522703
(HzHw vs LzHw) |  -.4474113   .2060755    -2.17   0.031    -.8538739   -.0409487
(HzLw vs LzLw) |   .0187116   .1722173     0.11   0.914    -.3209691    .3583924
(LzHw vs LzLw) |   .3944309   .2202323     1.79   0.075    -.0399546    .8288163
(HzHw vs LzLw) |  -.0529804    .163752    -0.32   0.747    -.3759642    .2700034
(HzLw vs LzHw) |  -.3757192   .2579781    -1.46   0.147    -.8845543    .1331159
The p-values in the table above cannot be taken at face value. We need to make some kind of an adjustment to account for the fact that there are six post-hoc tests. The easiest adjustment would be a Bonferroni. This is an extremely conservative adjustment and is just a matter of multiplying each of the p-values by the number of tests (6). In this case, the comparison with the smallest p-value, HzHw versus LzHw, would have a Bonferroni adjusted p-value of .031*6 = .186.

In this example, the p-value for the three-way interaction was barely significant (.047), so it is not surprising that none of the post-hoc tests of differences in simple slopes were statistically significant.

Below are all the commands that we used collected into a do-file to make it easier to use. You might be wondering, wouldn't it be even easier if there was an ado-file for that would do everything with a single command? It wouldn't be very difficult to create such a command but you would not be able to see how all of the steps work to create the tests of differences in simple slopes. In other words, it builds character. 

use http://www.ats.ucla.edu/stat/stata/notes/hsb2, clear

rename write y
rename read x
rename math w
rename science z

generate xz  = x*z
generate xw  = x*w
generate zw  = z*w
generate xzw = x*z*w

regress y x z w xz xw zw xzw

sum x
global hi=r(max)
global lo=r(min)
sum w
global Hw=r(mean)+r(sd)
global Lw=r(mean)-r(sd)
sum z
global Hz=r(mean)+r(sd)
global Lz=r(mean)-r(sd)

global HzHw "x + ($Hz)*xz + ($Hw)*xw + ($Hz)*($Hw)*xzw"
global HzLw "x + ($Hz)*xz + ($Lw)*xw + ($Hz)*($Lw)*xzw"
global LzHw "x + ($Lz)*xz + ($Hw)*xw + ($Lz)*($Hw)*xzw"
global LzLw "x + ($Lz)*xz + ($Lw)*xw + ($Lz)*($Lw)*xzw"

/* simple slopes */
lincom $HzHw /* slope 1 */
global b1 = r(estimate)
lincom $HzLw /* slope 2 */
global b2 = r(estimate)
lincom $LzHw /* slope 3 */
global b3 = r(estimate)
lincom $LzLw /* slope 4 */
global b4 = r(estimate)

/* intercepts */
lincom _cons + $Hz*z + $Hw*w + $Hz*$Hw*zw
global c1 = r(estimate)
lincom _cons + $Hz*z + $Lw*w + $Hz*$Lw*zw
global c2 = r(estimate)
lincom _cons + $Lz*z + $Hw*w + $Lz*$Hw*zw
global c3 = r(estimate)
lincom _cons + $Lz*z + $Lw*w + $Lz*$Lw*zw
global c4 = r(estimate)


twoway (function y=$c1+$b1*x, range($lo $hi))  ///
       (function y=$c2+$b2*x, range($lo $hi))  ///
       (function y=$c3+$b3*x, range($lo $hi))  ///
       (function y=$c4+$b4*x, range($lo $hi)), ///
       legend(order(1 "1 HzHw" 2 "2 HzLw" 3 "3 LzHw" 4 "4 LzLw")) name(con3way, replace)

/* differences in simple slopes */
lincom ($HzHw)-($HzLw)  /* a) HzHw vs HzLw */
lincom ($HzHw)-($LzHw)  /* b) HzHw vs LzHw */
lincom ($HzLw)-($LzLw)  /* c) HzLw vs LzLw */
lincom ($LzHw)-($LzLw)  /* d) LzHw vs LzLw */
lincom ($HzHw)-($LzLw)  /* e) HzHw vs LzLw */
lincom ($HzLw)-($LzHw)  /* f) HzLw vs LzHw */
Hey, wait a minute, my model also has covariates. What about them? We will deal with adding covariates to the model in Part 2: How can I explain a 3-way continuous interaction?

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

How to cite this page

Report an error on this page or leave a comment

The content of this web site should not be construed as an endorsement of any particular web site, book, or software product by the University of California.