查看原文
其他

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

Stata连享会 连享会 2022-12-31

作者:高娜娜(中南财经政法大学)

2019暑期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 表示个体的工作时间,而 gradettl 的乘积。

  • 计算边际效应

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

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

  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β3var(β1)var(β3)cov(β1 β3),分别赋给标量b1b3varb1varb3covb1b3

  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)含有三个变量的交乘项

其中,wagegradettlgrade_x_ttl 的含义同模型一。此外,该模型中加入了 union 变量,该变量表示该妇女是否为工会成员,分别用 1 和 0 表示。并加入 了 gradeunionttlunion 的交乘项 grade_unionttl_union,以及 gradettlunion 三个变量的交乘项 grade_x_ttl_union

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

我们需要分别估计 union=1union=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专栏。可以在上述网站中搜索关键词StataStata连享会后关注我们。

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

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

联系我们

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

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

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

  • 联系邮件: StataChina@163.com

往期精彩推文


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

一起学空间计量……

您可能也对以下帖子感兴趣

文章有问题?点此查看未经处理的缓存