Stata: 手动计算和图示边际效应

作者：高娜娜（中南财经政法大学）

特别说明

文中包含的链接在微信中无法生效。请点击本文底部左下角的【阅读原文】，转入本文【简书版】。

连享会-交乘项 (调节效应) 专题系列

• Stata：交乘项该如何使用？

• Stata：边际效应分析

• Stata：图示连续变量的边际效应(交乘项)

• Stata：图示交互效应\调节效应

在建立模型时，我们往往会根据一定的经济理论加入一些交乘项，在 Stata 里，用 margins 命令可以直接得到相关变量的边际效应，进而可以使用 marginsplot 图示边际效应，详情参见 [Stata：边际效应分析]。如果希望做更为深入的分析，并采用更为酷炫的图形展示，则可以使用外部命令 interflex，详情参见 [Stata：交叉项\交乘项该这么分析！]。至于有关交乘项更一般化的介绍，则可以参考 [Stata：交乘项该如何使用？] 和 [加入交乘项后符号变了!?]。

然而，「百闻不如一见，百见不如一干」。若能亲手计算一次边际效应，必能大大加深对齐背后原理的理解，在论文写作过程中也能够自信地应对各种奇葩结果的解释。

为此，本文采用 庖丁解牛 的方式，介绍一下如何手动计算交乘项的边际效应。在下一篇推文中，我们将进一步对包含三个交乘项的情形进行拆解。

1. 边际效应计算原理

1.1 两个变量的交乘项

在线性模型中，假设 X 和 Y 的关系会受到第三者 Z (通常被形象地称为 调节变量) 的影响，模型的基本形式如下：

所谓 X 的边际效应 是指 X 每增加 1 单位时 Y 的变化，即：

由上式可以看出， X 的边际效应不再是常数，它会随着 Z 的变化而变化。因此，在讨论 X 对 Y 的边际影响时，我们必须考虑 Z 的取值。

例如，假设 ，，则当  时，「X 对 Y 的边际影响」为：

多数情况下，我们会更多地关注当  (样本均值) 或特定分位数 (如第 25，50, 75 百分位数) 时， X 对 Y 的边际效应。

在该模型中，X 对 Y 的边际效应的标准误，即  的标准误也与 Z 的取值有关：

利用模型估计得到的 的方差、  的方差、 和 的协方差以及给定的 Z 值，可以得到 X 的标准误。计算标准误有助于在绘制边际效应图时附加上置信区间。

1.2 三个变量的交乘项

对于含有三个变量交乘项的模型（如下模型），也可以用上述原理来手动计算变量的边际效应。

该模型中 X 的边际效应和 ∂y/∂x 的方差的表达式如下所示：

特别地，一种特殊的情况是，在三个变量中有一个是 0-1 变量，比如下面模型中的 D 变量。

在这种情况下，我们可以根据 D 将样本分成两个样本组，即 G1 组 (D=1) 和 G2 组 (D=0) ，使每个组中只包含一个交乘项，然后对两个样本组分别进行回归分析并求得 X 的边际效应。

当我们把该模型拆成两个含有两个交乘项的模型时，变量 X 的边际效应和 的标准误的表达式与含有两个变量交乘项的模型一致，此处不再赘述。

2. Stata 实现

下面我们使用 Stata 自带的 nlsw88.dta 数据为例，手动计算交乘项的边际效应。

2.1 两个变量的交乘项

其中， wage 表示个体的工资水平， grade 表示已完成受教育的年限， ttl 表示个体的工作时间，而 是 grade 与 ttl 的乘积。

• 计算边际效应

在上述模型表示，已完成受教育年限会影响工资水平，而个体的工作经验又充当了调节变量，我们要计算的是已完成受教育年限对工资的边际效应。

首先，对上述模型进行估计。

1. sysuse "nlsw88.dta", clear

2. gen grade_x_ttl = grade*ttl

3. reg wage grade ttl grade_x_ttl

4. Source | SS df MS Number of obs = 2,244

5. -------------+---------------------------------- F(3, 2240) = 133.83

6. Model | 11301.2662 3 3767.08872 Prob > F = 0.0000

7. Residual | 63053.0643 2,240 28.1486894 R-squared = 0.1520

8. -------------+---------------------------------- Adj R-squared = 0.1509

9. Total | 74354.3305 2,243 33.1495009 Root MSE = 5.3055

10. ------------------------------------------------------------------------------

11. wage | Coef. Std. Err. t P>|t| [95% Conf. Interval]

12. -------------+----------------------------------------------------------------

13. grade | 0.252 0.132 1.91 0.056 -0.006 0.509

14. ttl_exp | -0.144 0.128 -1.12 0.264 -0.396 0.108

15. grade_x_ttl | 0.032 0.010 3.21 0.001 0.013 0.052

16. _cons | 0.934 1.658 0.56 0.573 -2.317 4.184

17. ------------------------------------------------------------------------------

模型中变量的系数，以及系数的方差协方差存储分别在矩阵 b 和矩阵 V 中。

1. . ereturn list // 列示内存中的返回值

2. matrices:

3. e(b) : 1 x 4

4. e(V) : 4 x 4

5. . matrix list e(V) // 呈现系数的方差-协方差矩阵

6. symmetric e(V)[4,4]

7. grade ttl_exp grade_x_ttl _cons

8. grade .0173019

9. ttl_exp .01534545 .01651049

10. grade_x_ttl -.00123247 -.00125844 .00009963

11. _cons -.21378964 -.19844023 .01533119 2.7477949

12. . dis %4.3f sqrt(.0173019) // 变量 grade 的系数标准误

13. 0.132

然后，从矩阵中提取我们需要的 β1、β3、var(β1)、var(β3)、cov(β1 β3)，分别赋给标量b1、b3、varb1、varb3、covb1b3。

1. matrix b = e(b)

2. matrix V = e(V)

3. scalar b1 = b[1,1]

4. scalar b3 = b[1,3]

5. scalar varb1 = V[1,1]

6. scalar varb3 = V[3,3]

7. scalar covb1b3= V[1,3]

最后，计算边际效应。如果我们要计算在 ttl 的均值下 grade 的边际效应，我们需要找到 ttl 的均值，然后带入到边际效应的计算公式里即可。当然，我们也可以带入其他的 ttl 值来得到 grade 的边际效应。

1. sum ttl

2. scalar ttl_mean = r(mean)

3. dis "Margins：" b1+b3*ttl_mean

4. dis "Std.Err：" sqrt(varb1+varb3*(ttl_mean^2)+2*covb1b3*ttl_mean)

1. Margins：.65359238

2. Std.Err：.04536131

• 绘制边际效应图

为了绘制 grade 的边际效应图，我们需要知道 grade 的边际效应如何随着 ttl 变化。

首先，我们需要得到一系列连续变化的 ttl 值，在数据中，ttl 的最小值是 0.16，最大值是 28.9，用以下命令得到以 0.01 为间隔从 0.16 变化到 28.9 的序列。

1. set obs 2875

2. generate MVZ=(_n+15)/100

然后，分别计算每一个 MVZ 上 X 的边际效应 conbx、标准误 consx、5%置信区间的上界 upperx和下界 lowerx。

1. gen conbx=b1+b3*MVZ

2. gen consx=sqrt(varb1+varb3*(MVZ^2)+2*covb1b3*MVZ)

3. gen ax=1.96*consx

4. gen upperx=conbx+ax

5. gen lowerx=conbx-ax

为了图形的美观，我们构建了一些辅助值。这些辅助值产生的效果可以参见后文绘制的图形。

1. gen where=-0.045

2. gen pipe = "|"

3. gen yline=0

1. graph twoway hist ttl, width(0.5) percent color(gs14) yaxis(2) ///

2. || scatter where ttl , plotr(m(b 4)) ms(none) mlabcolor(gs5) mlabel(pipe) mlabpos(6) legend(off) ///

3. || line conbx MVZ,clpattern(solid) clwidth(medium) clcolor(black) yaxis(1) ///

4. || line upperx MVZ,clpattern(dash) clwidth(thin) clcolor(black) ///

5. || line lowerx MVZ,clpattern(dash) clwidth(thin) clcolor(black) ///

6. || line yline MVZ,clpattern(solid) clwidth(thin) clcolor(black) ///

7. || , ///

8. xlabel( 0 5 10 15 20 25 30, nogrid labsize(2)) ///

9. ylabel(-.2 0 .2 .4 .6 0.8 1 1.2 1.4 1.6 , axis(1) nogrid labsize(2)) ///

10. ylabel(0 1 2 3 4 5, axis(2) nogrid labsize(2)) ///

11. yscale(noline alt) ///

12. yscale(noline alt axis(2)) ///

13. xscale(noline) ///

14. legend(off) ///

15. xtitle("" , size(2.5)) ///

16. ytitle("" , axis(2) size(2.5)) ///

17. xsca(titlegap(2)) ///

18. ysca(titlegap(2)) ///

19. scheme(s2mono) graphregion(fcolor(white) ilcolor(white) lcolor(white))

连享会-手动绘制边际效应图形

（2）含有三个变量的交乘项

其中，wage、grade、ttl、grade_x_ttl 的含义同模型一。此外，该模型中加入了 union 变量，该变量表示该妇女是否为工会成员，分别用 1 和 0 表示。并加入 了 grade 和 union、ttl 和 union 的交乘项 grade_union 、ttl_union，以及 grade、ttl 和 union 三个变量的交乘项 grade_x_ttl_union。

在该模型中，union 是 0-1 变量，我们根据 union 将模型二改写成如下形式：

我们需要分别估计 union=1 和 union=0 时的模型，并分别计算 grade 边际效应。方法同模型一，此处直接附上代码和图形。

1. **计算边际效应：union==1**

2. reg wage grade ttl grade_x_ttl if union==1

3. matrix b_1=e(b)

4. matrix V_1=e(V)

5. scalar b1_1=b_1[1,1]

6. scalar b3_1=b_1[1,3]

7. scalar varb1_1=V_1[1,1]

8. scalar varb3_1=V_1[3,3]

9. scalar covb1b3_1=V_1[1,3]

10. sum ttl if union==1

11. scalar ttl_mean_1=r(mean)

12. dis "Margins_1：" b1_1+b3_1*ttl_mean_1

13. dis "Std.Err_1：" sqrt(varb1_1+varb3_1*(ttl_mean_1^2)+2*covb1b3_1*ttl_mean_1)

14. gen conbx_1=b1_1+b3_1*MVZ

15. gen consx_1=sqrt(varb1_1+varb3_1*(MVZ^2)+2*covb1b3_1*MVZ)

16. gen ax_1=1.96*consx_1

17. gen upperx_1=conbx_1+ax_1

18. gen lowerx_1=conbx_1-ax_1

19. **计算边际效应：union==0**

20. reg wage grade ttl grade_x_ttl if union==0

21. matrix b_0=e(b)

22. matrix V_0=e(V)

23. scalar b1_0=b_0[1,1]

24. scalar b3_0=b_0[1,3]

25. scalar varb1_0=V_0[1,1]

26. scalar varb3_0=V_0[3,3]

27. scalar covb1b3_0=V_0[1,3]

28. sum ttl if union==0

29. scalar ttl_mean_0=r(mean)

30. dis "Margins_0：" b1_0+b3_0*ttl_mean_0

31. dis "Std.Err_0：" sqrt(varb1_0+varb3_0*(ttl_mean_0^2)+2*covb1b3_0*ttl_mean_0)

32. gen conbx_0=b1_0+b3_0*MVZ

33. gen consx_0=sqrt(varb1_0+varb3_0*(MVZ^2)+2*covb1b3_0*MVZ)

34. gen ax_0=1.96*consx_0

35. gen upperx_0=conbx_0+ax_0

36. gen lowerx_0=conbx_0-ax_0

37. **绘图**

38. graph twoway line conbx_1 MVZ, clpattern(solid) clwidth(medium) clcolor(blue) ///

39. || line upperx_1 MVZ, clpattern(dash) clwidth(thin) clcolor(blue) ///

40. || line lowerx_1 MVZ, clpattern(dash) clwidth(thin) clcolor(blue) ///

41. || line conbx_0 MVZ, clpattern(solid) clwidth(medium) clcolor(red) ///

42. || line upperx_0 MVZ, clpattern(dash) clwidth(thin) clcolor(red) ///

43. || line lowerx_0 MVZ, clpattern(dash) clwidth(thin) clcolor(red) ///

44. || line yline MVZ, clpattern(solid) clwidth(thin) clcolor(black) ///

45. || , ///

46. xlabel( 0 5 10 15 20 25 30, nogrid labsize(2)) ///

47. ylabel(-.2 0 .2 .4 .6 0.8 1 1.2 1.4 1.6 , nogrid labsize(2)) ///

48. xscale(noline) ///

49. legend(off) ///

50. xtitle("" , size(2.5)) ///

51. xsca(titlegap(2)) ///

52. ysca(titlegap(2)) ///

53. scheme(s2mono) graphregion(fcolor(white) ilcolor(white) lcolor(white))

分组绘制边际效应-代替三个交乘项

参考资料：Marginal Effect Plot for X: An Interaction Between X and Z

关于我们

• 【Stata 连享会(公众号：StataChina)】由中山大学连玉君老师团队创办，旨在定期与大家分享 Stata 应用的各种经验和技巧。

• 公众号推文同步发布于 CSDN-Stata连享会 、简书-Stata连享会 和 知乎-连玉君Stata专栏。可以在上述网站中搜索关键词Stata或Stata连享会后关注我们。

• 点击推文底部【阅读原文】可以查看推文中的链接并下载相关资料。

• Stata连享会 精彩推文1  || 精彩推文2

联系我们

• 欢迎赐稿： 欢迎将您的文章或笔记投稿至Stata连享会(公众号: StataChina)，我们会保留您的署名；录用稿件达五篇以上，即可免费获得 Stata 现场培训 (初级或高级选其一) 资格。

• 意见和资料： 欢迎您的宝贵意见，您也可以来信索取推文中提及的程序和数据。

• 招募英才： 欢迎加入我们的团队，一起学习 Stata。合作编辑或撰写稿件五篇以上，即可免费获得 Stata 现场培训 (初级或高级选其一) 资格。

• 联系邮件： [e-mail]

往期精彩推文

---

欢迎加入Stata连享会(公众号: StataChina)