🍓 课程推荐：2024 机器学习与因果推断专题
主讲老师：司继春；张宏亮
课程时间：2024 年 11 月 9-10 日 ；16-17日
课程咨询：王老师 18903405450（微信）
作者：陈希 (中国农业科学院农经所)
邮箱：chenxi_alice@163.com
温馨提示： 文中链接在微信中无法生效。请点击底部「阅读原文」。或直接长按/扫描如下二维码，直达原文：
目录
1. 写在前面
2. 对三项交互的解释
3. Stata 实操
3.1 基础模型
3.2 加入协变量的情况
4. 参考资料
5. 相关推文
1. 写在前面
多元回归模型中经常包含交互项，也会出现三个连续变量交互的情况，但是应该如何理解三项交互的含义呢？
我们摘译了 UCLA FAQs 中的相关介绍，从斜率的角度加以解释：How can I understand a 3-way continuous interacation?。
2. 对三项交互的解释
三项交互的情况，即两个调节变量共同影响因变量在自变量上的回归。
我们可以从多个角度解读三项交互的含义，本文主要介绍基于边际效应 (斜率) 视角的理解，即当两个调节变量在高值和低值的不同组合下保持恒定时，因变量在自变量上的斜率。在计算出简单斜率之后，我们将计算并检验所有斜率之间的差异。本文使用的方法主要源于 Dawson 和 Richter (2004) 的论文。Dawson 和 Richter 使用许多复杂公式计算了简单斜率的差异及其标准误差。我们将采取一种更加简单的方式来进行说明。
我们将首先研究回归方程，其中包括三项连续相互作用。在公式中， 是因变量， 是核心变量， 和 分别是两个调节变量。
我们可以将方程中的各项重新排列组合为两组，第一组（不包含 的项）为截距，而第二组（所有包含 的项）定义为简单斜率。
接下来，我们将 和 的高值定义为高于其各自均值一个标准差，并将它们表示为 和 ；低值定义为低于其均值一个标准差，分别为 和 。由此，共形成四种可能的条件组合：
其中，条件 1 表示为当 和 都处于其高值时的斜率公式，具体为：
3. Stata 实操
3.1 基础模型
我们将使用数据集 hsb2.dta 来说明。为了保持符号的一致性，我们将暂时更改变量的名称：上文中， 表示因变量， 表示自变量， 表示第一个调节变量， 表示第二个调节变量。从下面的代码中可以清楚地看出，write 是因变量，read 是自变量，math 和 science 是调节变量。
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
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.014 1.803 -1.67 0.096 -6.570 0.542
z | -2.435 1.757 -1.39 0.167 -5.902 1.031
c.x#c.z | 0.055 0.033 1.70 0.091 -0.009 0.119
w | -3.228 1.926 -1.68 0.095 -7.028 0.572
c.x#c.w | 0.074 0.036 2.05 0.041 0.003 0.145
c.z#c.w | 0.061 0.035 1.74 0.083 -0.008 0.130
c.x#c.z#c.w | -0.001 0.001 -2.00 0.047 -0.002 -0.000
_cons | 166.544 93.128 1.79 0.075 -17.141 350.229
------------------------------------------------------------------------------
请注意这个 xzw 的 P 值为0.047，在5%的显著性水平下显著。根据上文，我们将创建四个条件，分别是当 w 和 z 比平均值高一个标准差，以及其比平均值低一个标准差的排列组合。
/* determine range of x for graphing */
sum x
sum w
/* determined high and low values for w */
global Hw=r(mean)+r(sd)
global Lw=r(mean)-r(sd)
sum z
/* determine high and low values for z */
global Hz=r(mean)+r(sd)
global Lz=r(mean)-r(sd)
我们将使用 margins 命令来计算 w 和 z 这四个值的简单斜率。我们使用 dydx 选项来获得斜率。值得注意的是，在命令的括号中我们加入了 w 和 z 的各种组合。
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. t P>|t| [95% Conf. Interval]
---------+--------------------------------------------------------
x |
_at |
1 | 0.165 0.095 1.73 0.086 -0.024 0.353
2 | 0.236 0.143 1.65 0.100 -0.046 0.518
3 | 0.612 0.196 3.12 0.002 0.225 0.999
4 | 0.218 0.128 1.71 0.090 -0.034 0.469
------------------------------------------------------------------
根据上述结果，我们发现 HwHz, LwHz 和 LwLz 条件下的简单斜率都非常相似。但是， Hwlz的斜率似乎与其他三种都不太相同。下一步是绘制斜率以及观测值的散点图。我们将再次使用 margins 命令来获取绘图信息。由于我们使用的是线性回归，因此我们只需要两个点来定义每个简单的斜率回归线。我们将使用 x 值 20 和 80，它们涵盖了变量 x 的全部范围。
margins, at(x=(20 80) w=($Hw $Lw) z=($Hz $Lz)) vsquish
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. t P>|t| [95% Conf. Interval]
-----+----------------------------------------------------------
_at |
1 | 53.294 4.059 13.13 0.000 45.288 61.301
2 | 44.503 5.033 8.84 0.000 34.575 54.431
3 | 35.420 6.803 5.21 0.000 22.003 48.837
4 | 39.981 3.034 13.18 0.000 33.997 45.965
5 | 63.168 1.927 32.79 0.000 59.367 66.968
6 | 58.678 4.007 14.65 0.000 50.775 66.580
7 | 72.138 5.488 13.14 0.000 61.313 82.963
8 | 53.033 4.787 11.08 0.000 43.591 62.475
----------------------------------------------------------------
接下来，我们使用 marginsplot 命令进行绘图。addplot 选项允许包括观测值的散点图以及简单的斜率回归线。
marginsplot, recast(line) noci addplot(scatter y x, msym(oh) jitter(3))
marginsplot, recast(line) noci addplot(scatter y x, msym(oh) jitter(3))
3.2 加入协变量的情况
我们将在这一节中处理向模型添加协变量的情况。下面是上一页的回归模型，其中添加了两个协变量 V~1~ 和 **V~2~**。
向模型添加协变量不会更改简单斜率的公式，只要协变量不与自变量或调节变量交互即可。我们在之前的模型中添加两个协变量，socst 重命名为 V~1~，ses 将重命名为 **V~2~**。本节与上一示例完全相同，只是将协变量添加到 regress 命令中。
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.463 1.717 -2.02 0.045 -6.849 -0.077
z | -2.902 1.670 -1.74 0.084 -6.197 0.392
c.x#c.z | 0.062 0.031 2.02 0.045 0.001 0.123
w | -3.697 1.831 -2.02 0.045 -7.309 -0.084
c.x#c.w | 0.080 0.034 2.33 0.021 0.012 0.147
c.z#c.w | 0.070 0.033 2.09 0.038 0.004 0.135
c.x#c.z#c.w | -0.001 0.001 -2.32 0.021 -0.003 -0.000
v1 | 0.279 0.058 4.84 0.000 0.165 0.393
v2 | -0.906 0.691 -1.31 0.192 -2.270 0.458
_cons | 185.313 88.564 2.09 0.038 10.618 360.009
----------------------------------------------------------------------
值得注意的是，包含协变量后，三项相互作用甚至比之前更加显著。接下来，我们将计算 w 和 z 的各种组合的简单斜率。margins 命令与没有协变量的示例完全相同，因为协变量不参与交互作用。
/* 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. t P>|t| [95% Conf. Interval]
---------+------------------------------------------------------
x |
_at |
1 | 0.036 0.095 0.37 0.708 -0.152 0.223
2 | 0.144 0.137 1.05 0.296 -0.127 0.414
3 | 0.498 0.188 2.64 0.009 0.127 0.870
4 | 0.093 0.124 0.75 0.453 -0.151 0.337
----------------------------------------------------------------
让我们在观测值的顶部绘制四个简单的斜率。同样的，我们只需要两个 x=20 和 x= 80 两个点来定义每条斜率回归线。
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. t P>|t| [95% Conf. Interval]
-----+----------------------------------------------------------
_at |
1 | 56.987 3.937 14.47 0.000 49.221 64.753
2 | 47.840 4.826 9.91 0.000 38.321 57.359
3 | 38.204 6.501 5.88 0.000 25.381 51.027
4 | 44.624 3.040 14.68 0.000 38.627 50.621
5 | 59.127 2.060 28.70 0.000 55.063 63.190
6 | 56.454 3.835 14.72 0.000 48.888 64.019
7 | 68.113 5.302 12.85 0.000 57.654 78.572
8 | 50.202 4.581 10.96 0.000 41.166 59.237
----------------------------------------------------------------
marginsplot, recast(line) noci addplot(scatter y x, msym(oh) jitter(3))
marginsplot, recast(line) noci addplot(scatter y x, msym(oh) jitter(3))
如图所示，模型中加入协变量后，我们发现简单斜率有了一些变化。但是，斜率 3 还是与其他斜率有着较为不同的坡度。我们可以通过 margins 命令的 pwcompare(effects) 选项计算斜率之间的差异。
/* 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 Number of obs = 200
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. t P>|t| [95% Conf. Interval]
----------+-----------------------------------------------------------
x |
_at |
2 vs 1 | 0.108 0.156 0.69 0.490 -0.200 0.416
3 vs 1 | 0.463 0.196 2.37 0.019 0.077 0.848
4 vs 1 | 0.057 0.156 0.37 0.713 -0.250 0.365
3 vs 2 | 0.355 0.245 1.45 0.149 -0.128 0.838
4 vs 2 | -0.051 0.164 -0.31 0.757 -0.373 0.272
4 vs 3 | -0.406 0.210 -1.93 0.055 -0.819 0.008
----------------------------------------------------------------------
根据该表显示，似乎只有斜率 3 与斜率 1 (LzHw vs HzHw) 之间有显著差异。但在我们下此结论之前，我们需要调整 p 值，以考虑到这些是事后检验。为了检验，我们需要在 margins 命令中加入 mcompare(bonferroni) 选项。
/* 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 Number of obs = 200
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. t P>|t| [95% Conf. Interval]
-----------+------------------------------------------------------------
x |
_at |
2 vs 1 | 0.108 0.156 0.69 1.000 -0.308 0.524
3 vs 1 | 0.463 0.196 2.37 0.114 -0.058 0.984
4 vs 1 | 0.057 0.156 0.37 1.000 -0.358 0.473
3 vs 2 | 0.355 0.245 1.45 0.894 -0.298 1.008
4 vs 2 | -0.051 0.164 -0.31 1.000 -0.487 0.386
4 vs 3 | -0.406 0.210 -1.93 0.328 -0.965 0.154
------------------------------------------------------------------------
使用 Bonferroni 校正后，斜率 3 vs 斜率 1 的 p 值变为 0.114，在统计意义上并不显著。由于 Bonferroni 校正过于保守，可能不是最佳选择，但它比较容易实施。
4. 参考资料
How can I understand a 3-way continuous interaction?, UCLA FAQs, -Link-
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.
Dawson, J. F., A. W. Richter, 2006, Probing three-way interactions in moderated multiple regression: Development and application of a slope difference test, Journal of Applied Psychology, 91 (4): 917. -Link-, -PDF-, -cited-2000+
5. 相关推文
Note：产生如下推文列表的 Stata 命令为：
lianxh 交乘 交互 边际效应
安装最新版lianxh命令：
ssc install lianxh, replace
专题：回归分析 调节效应是否需要考虑对控制变量交乘？ Stata：边际效应知多少？f_able命令（上） Stata：边际效应知多少？f_able命令（下） Stata: 手动计算和图示边际效应 Stata因子变量：虚拟变量-交乘项批量处理 Stata: 边际效应分析 专题：内生性-因果推断 Stata：内生变量的交乘项如何处理？ 专题：交乘项-调节-中介 我为什么越跑越肥？交乘项系数的直观解释 Logit-Probit中的交乘项及边际效应图示 内生变量的交乘项如何处理？ interactplot：图示交乘项-交互项-调节效应 Stata：交乘项的对称效应与图示 Stata：交乘项该如何使用？-黄河泉老师PPT Stata：虚拟变量交乘项生成和检验的简便方法 Stata：图示连续变量的连续边际效应 Stata：内生变量和它的交乘项 交乘项-交叉项的中心化问题 交乘项专题：主效应项可以忽略吗？ Stata：交乘项该这么分析! 专题：Probit-Logit Logit-Probit：非线性模型中交互项的边际效应解读 专题：分位数回归 Stata分位数回归I：理解边际效应和条件边际效应
🍉 扫码加入连享会微信群，提问交流更方便
