Stata：图示边际效应
作者：连玉君 | 杨柳 ( 知乎 | 简书 | 码云 )
特别说明： 文中包含的链接在微信中无法生效，部分输出结果有删减。请点击本文底部左下角的【阅读原文】，转入本文【简书版】。
前期相关推文：
1. 边际效应简介
1.1 引言
考研复试结束后，你和闺蜜决定去成都旅游。当你和闺蜜正在品尝当地最有名气的麻辣火锅时，你们感觉心情非常愉快。此时，麻辣火锅将有助于有一个好心情，它对心情的边际效应是正值；当你们吃到一半时，手机上收到了一条消息，是考研复试的排名结果，打开消息后，发现你们俩都榜上有名。你们看到被录取的消息后万分高兴，吃的更 high 了，又多点了一些菜，并决定吃完后再去 KTV 庆祝下。此时，考研成功的结果极大的增加了吃火锅对心情的边际效应值。
研究生入学后，你和闺蜜都十分努力学习，认真的完成导师布置的课题任务并将课题研究内容整理成小论文投到了一个 C 刊上，但遭到拒稿。此时，努力学习对科研成果的边际效应是负值。不过你们没有就此停止努力，继续请教导师、按照意见认真修改并写成英文投到了一个 SSCI 期刊上，结果被录用了。此时，随着努力程度的增加，它对科研成果的边际效应变为正值。
总结上面的例子，我们发现吃麻辣火锅（x1）对心情（y）的边际效应受到其他变量（x2：考研成功）的调节作用，使得该边际效应值增加了。一开始，努力学习（x3）对科研成果的边际效应为负，但随着努力程度的增加（x3 值的增加），它对科研成果的边际效应变为正值。
1.2 边际效应分析的必要性
虽然 回归结果表格 中的变量的系数估计值反映了该变量对被解释变量影响作用的大小，并且一直是学者们交流回归模型结果的重要方式，但是，当回归模型中包含 类别变量、交乘项 或者 回归模型为非线性（诸如 Logit, Probit 等非线性模型）时，对系数估计值的解释就非常具有挑战性。这时，就需要计算变量的 边际效应 或者计算 预测边际值， 以探求 自变量变化对 因变量变化 的 影响作用 或分析比较不同情况时的因变量预测边际值的大小。
下面，我们就一起来学习如何在 Stata 中计算边际效应并绘制图形。
1.3 边际效应的定义
所谓 边际效应，是从已有拟合模型结果中计算出来的统计量，该数值表示 自变量的变化 对 因变量的变化 的 影响作用 的 大小。
在对模型结果进行分析时，可以计算 连续变量取某一个值 时，连续变量 对 因变量 的边际效应，也可以计算 连续变量平均值处 的边际效应，或者还可以计算 其他变量取均值 时 连续变量对 因变量 的平均边际效应。
2. 计算边际效应命令(margins)与绘图命令(marginsplot)
2.1 margins 命令
margins命令的语法如下所示：
margins [marginlist] [if] [in] [weight] [, response_options options]命令含义：
图1 margins 命令解释
Note: 有关如何使用因子变量的介绍请参见 往期推文：Stata: 因子变量全攻略。
2.2 marginsplot 绘图命令
使用
marginsplot命令可以将之前刚刚计算的边际效应的结果以图的形式展示出来。语法如下：marginsplot [, options]
图2 marginsplot 命令解释
3. margins 与 marginsplot 命令举例
3.1 基础案例
在美国的种族文化中，不可否认白人与黑人之间的差异性。例如，存在着白人与黑人在行业类别与工资方面差异的现象。于是，我们想检验种族是否为工资的显著影响因素，还想了解当行业类别相当时，不同种族的工资的平均水平分别是多少，它们之间的差别有多少。
接下来，我们使用 Stata 的自带数据 nlsw88.dta (1988年美国妇女小时工资) ，以 wage (妇女的小时工资) 作为被解释变量、以 industry (行业类别)、 race (种族类别) 作为解释变量建立线性回归模型。
race 变量为 类别变量，它包括三个类别，分别为 white、black、other。可以使用 因子变量 的语法格式，在变量前面加上前缀 i. 生成虚拟变量 (i.race)，基准组为第一个类别 white。industry 变量也同样使用 因子变量 的语法格式生成虚拟变量 (i.industry)，基准组为第一个类别 Ag/Forestry/Fisheries。 Stata 中的回归命令和结果如下所示：
. sysuse "nlsw88.dta", clear(NLSW, 1988 extract)
. reg wage i.race i.industry我们还想进一步了解当控制行业类别变量时，各个种族类别（white, black, other）的妇女的小时工资的平均水平是多少。
于是，我们使用 margins 命令附加 atmeans 的选项就可以计算当 其他变量取均值时，不同种族的妇女 的小时工资的 预测边际值。
Stata 命令如下：
. margins i.race, atmeans计算结果表明：当行业类别变量取均值时，白人妇女的小时工资的预测边际值为 8.067、黑人妇女的小时工资的预测边际值为 6.967、其他种族的妇女的小时工资的预测边际值为 8.198。
我们还想将边际效应的计算结果用图的形式表示。使用 marginsplot 命令就可以很方便的实现这个想法。Stata 命令如下所示：
. marginsplot图3 各种族类别的妇女小时工资的预测边际值
3.2 交乘项案例
下面，继续以研究妇女工资的影响因素为例对计算交乘项的边际效应的使用方法进行说明。
3.2.1 类别变量与类别变量交乘
从 3.1 节案例的回归结果中，我们已知道种族 race 是妇女工资 wage 的影响因素之一。
但是，除此之外，还有诸多因素会影响妇女工资，例如是否大学毕业 collgrad。因此，我们想探究是否大学毕业 collgrad 能否调节种族 race 对妇女工资 wage 的影响作用。于是，拟在回归模型中加入这两个变量的 交乘项 来检验是否存在调节效应。
使用 Stata 的自带数据
nlsw88.dta (1988年美国妇女小时工资)，以wage (妇女的小时工资)作为被解释变量、以industry (行业类别)、collgrad (是否大学毕业)、race (种族类别)、race (种族类别)与collgrad (是否大学毕业)的交乘项作为解释变量建立线性回归模型。使用
因子变量的语法格式，collgrad##i.race表示在模型中既包括collgrad与race变量，还包括collgrad与race变量的交乘项。Stata 中的回归命令和结果如下所示：
. reg wage i.industry collgrad##i.race
Source | SS df MS Number of obs = 2,232
--------------------+------------------------------- F(16, 2215) = 23.20
Model | 10639.3304 16 664.958149 Prob > F = 0.0000
Residual | 63477.9484 2,215 28.658216 R-squared = 0.1435
--------------------+------------------------------- Adj R-squared = 0.1374 Total | 74117.2788 2,231 33.2215503 Root MSE = 5.3533
-------------------------------------------------------------------------------
wage | Coef. Std. Err. t P>|t| [95% Conf. Interval]
--------------------+----------------------------------------------------------
collgrad |
college grad | 3.346 0.320 10.47 0.000 2.719 3.973
|
race |
black | -1.280 0.295 -4.34 0.000 -1.858 -0.701
other | -0.327 1.311 -0.25 0.803 -2.897 2.244
|
collgrad#race |
collegegrad#black | 2.099 0.659 3.18 0.001 0.806 3.391
collegegrad#other | 0.714 2.237 0.32 0.750 -3.673 5.102
|
_cons | 5.208 1.301 4.00 0.000 2.657 7.760
-------------------------------------------------------------------------------回归结果显示：大学毕业与黑人的交乘项 collgrad#black 的系数显著为正，而黑人 black的系数显著为负，表明大学毕业 collgrad 对黑人 black 与 妇女的小时工资 wage 之间的影响关系具有调节作用。
因此，我们想进一步了解 大学毕业(collgrad) 与 种族(race) 交乘项的 各个类别 对 妇女的小时工资 (wage) 的边际效应分别是多少。于是，我们使用 margins 附加 atmeans 的选项就可以计算当 其他变量取均值时，collgrad 与race交乘项的各个类别 的小时工资的 预测边际值。
Stata 命令和结果如下所示：
. margins collgrad#i.race, atmeans
Adjusted predictions Number of obs = 2,232
Model VCE : OLS
Expression : Linear prediction, predict()
at : 1.industry = .0076165 (mean)
( ……omitted…… )
0.collgrad = .7629928 (mean)
1.collgrad = .2370072 (mean)
1.race = .7289427 (mean)
2.race = .2594086 (mean)
3.race = .0116487 (mean)
--------------------------------------------------------------
| Delta-method
| Margin Std. Err. t P>|t|
------------------------+-------------------------------------
collgrad#race |
not college grad#white | 7.226 0.156 46.26 0.000
not college grad#black | 5.947 0.250 23.75 0.000
not college grad#other | 6.900 1.301 5.30 0.000
college grad#white | 10.573 0.273 38.71 0.000
college grad#black | 11.391 0.534 21.35 0.000
college grad#other | 10.960 1.789 6.13 0.000
--------------------------------------------------------------计算结果表明：当其他变量处于均值水平时，当妇女没有大学毕业时，白人妇女的小时工资的预测边际值最大，为 7.226，而当妇女为大学毕业时，黑人妇女的小时工资的预测边际值最大，为 11.391。
我们使用 marginsplot 命令来图示计算结果。Stata 命令如下所示：
. marginsplot图4 collgrad与race交乘项各个类别的妇女小时工资的预测边际值
此外，我们还可以进一步计算在 不同种族类别下(race)，大学毕业(collgrad=1) 与 非大学毕业(collgrad=0) 的 妇女小时工资 (wage) 的预测边际值的差值是多少。
此时，可以使用 margins 附加 dydx、at 与 atmeans 的选项来实现 。Stata 命令如下：
. margins, dydx(collgrad) at(race=(1 2 3)) atmeans计算结果显示：大学毕业(collgrad=1) 黑人妇女的小时工资的预测边际值与 非大学毕业(collgrad=0) 黑人妇女的小时工资的预测边际值的差值是最大的，为 5.444。简言之，上大学将在很大程度上提高黑人妇女工资。
同样的，可以使用 marginsplot 命令图示计算结果。Stata 命令如下所示：
. marginsplot图5 不同种族的大学毕业与非大学毕业的妇女小时工资的预测边际值的差别
3.2.2 类别变量与连续型变量交乘
影响妇女工资的因素较多，下面我们就来检验诸如 hours (每周工作时间) 与 union (是否工会成员) 这两个变量的 交乘项 是否会对妇女工资产生影响，其中 hours 为连续型变量， union 为类别变量。
使用 Stata 的自带数据
nlsw88.dta (1988年美国妇女小时工资)，以wage (妇女的小时工资)作为被解释变量、以industry (行业类别)、union (是否工会成员)、hours (每周工作小时数)、union (是否工会成员)与hours (每周工作小时数)的交乘项作为解释变量建立线性回归模型。使用
因子变量的语法格式，i.union##c.hours表示在模型中既包括union与hours变量，还包括union与hours变量的交乘项。Stata 中的回归命令和结果如下所示：
. reg wage i.industry c.hours##i.union
Source | SS df MS Number of obs = 1,864
-------------+---------------------------------- F(14, 1849) = 19.48
Model | 4165.55214 14 297.539439 Prob > F = 0.0000
Residual | 28235.7439 1,849 15.2708188 R-squared = 0.1286
-------------+---------------------------------- Adj R-squared = 0.1220 Total | 32401.296 1,863 17.392 Root MSE = 3.9078
-------------------------------------------------------------------------
wage | Coef. Std. Err. t P>|t| [95% CI]
----------------+--------------------------------------------------------
(output omitted)
hours | 0.057 0.010 5.42 0.000 0.036 0.077
union |
union | 3.761 0.899 4.18 0.000 1.998 5.524
union#c.hours |
union | -0.075 0.023 -3.30 0.001 -0.119 -0.030
_cons | 3.864 1.194 3.24 0.001 1.523 6.205
-------------------------------------------------------------------------回归结果显示 hours(每周工作小时数) 的系数值显著为正，为 0.056，union(是否工会成员) 的系数值显著为正，为 3.761，而 工会成员的每周工作小时(union#c.hours) 的系数值显著为负，为 -0.074，表明 hours (每周工作小时数) 对 wage (妇女的小时工资) 的边际效应会受到 union (是否工会成员) 的影响；union (是否工会成员) 对 wage (妇女的小时工资) 的边际效应也会受到 hours (每天工作小时数) 的影响。
可进一步使用 margins 命令附加 dydx 选项与 at 选项来计算当妇女为工会成员或非工会成员时，hours 对 wage 的平均边际效应分别为多少：
. margins, dydx(hour) at(union=(0 1))
Average marginal effects Number of obs = 1,864
Model VCE : OLS
Expression : Linear prediction, predict()
dy/dx w.r.t. : hours
1._at : union = 0
2._at : union = 1
------------------------------------------------
| Delta-method
| dy/dx Std. Err. t P>|t|
-------+----------------------------------------
hours |
_at |
1 | .0565515 .0104406 5.42 0.000
2 | -.0182076 .020211 -0.90 0.368
------------------------------------------------计算结果表明：当妇女为工会成员时，每周工作小时数增加 1 个单位，则小时工资下降0.018个单位，但在统计上不显著；当妇女为非工会成员时，每周工作小时增加 1 个单位，则小时工资将显著增加 0.056 个单位。
使用 marginsplot 命令将图示计算结果：
. marginsplot图6 是否工会成员的工作小时数对妇女小时工资的平均边际效应
我们还可以计算当妇女每周工作小时数不同时，
union对wage的边际效应分别为多少。可以使用margins命令附加dydx选项与at选项。在计算之前，我们需要事先知道hours变量的取值范围，可使用sum命令查看。
. keep if e(sample)
(382 observations deleted)
. sum hours
Variable | Obs Mean Std. Dev. Min Max-------------+---------------------------------------------------------
hours | 1,864 37.62071 9.959845 1 80
. margins, dydx(union) at(hours=(1(5)80))使用 marginsplot 命令进行绘图：
图7 工会成员与非工会成员的工作小时对妇女小时工资的边际效应的差别
3.2.3 连续型变量与连续型变量交乘
在实证研究中常常会分析两个连续型变量的交乘项的影响作用和变量的调节作用。
例如，当车辆重量 weight 与每加仑汽油行驶的距离 mpg 增加时，汽车价格 price 会有所增加。
现在，我们想进一步了解每加仑汽油行驶距离 mpg 能否调节车辆重量 weight 与汽车价格 price 之间的影响关系。于是，拟在回归模型中加入这两个 连续变量 的 交乘项，然后再计算当每加仑汽油行驶距离 mpg 取不同的数值时，车辆重量 weight 对汽车价格 price 的边际效应。
使用 Stata 的自带数据
auto.dta (1978年美国汽车数据)，以price (汽车价格)作为被解释变量、以foreign (是否进口车)、mpg (每加仑汽油能够行驶的英里数)、weight (汽车重量)、mpg (每加仑汽油能够行驶的英里数)与weight (汽车重量)的交乘项作为解释变量建立线性回归模型。使用
因子变量的语法格式c.mpg##c.weight表示在模型中既包括mpg与weight变量，还包括mpg与weight变量的交乘项。Stata 中的回归命令和结果如下所示：
. sysuse "auto.dta", clear. reg price foreign c.mpg##c.weight
. (output omitted)
price | Coef. Std. Err. t P>|t|
---------------+-----------------------------------------
foreign | 3369.814 691.422 4.87 0.000
mpg | 292.830 162.298 1.80 0.076
weight | 5.383 1.199 4.49 0.000
|
c.mpg#c.weight | -0.119 0.064 -1.87 0.066
|
_cons | -1.01e+04 4023.204 -2.51 0.014
---------------------------------------------------------可见，mpg 与 weight 的系数值显著为正，而 c.mpg#c.weight 的系数值显著为负 -0.118，表明 mpg (每加仑汽油能够行驶的英里数) 能够调节 weight (汽车重量) 对 汽车价格 (price) 的边际影响。
下面，我们计算当 mpg (每加仑汽油能够行驶的英里数) 取不同数值时，weight (汽车重量) 对 price (汽车价格) 的边际效应。使用 margins 命令附加 dydx 选项与 at 选项来计算。
首先，需要知道 mpg (每加仑汽油能够行驶的英里数) 的取值范围，因此，先使用 sum 命令查看该变量的基本统计量，再使用 margins 命令附加 dydx 选项与 at 选项。Stata 中的命令和结果如下所示：
. keep if e(sample)
(0 observations deleted)
. sum mpg //查看 mpg 的基本统计量
Variable | Obs Mean Std. Dev. Min Max-------------+---------------------------------------------------------
mpg | 74 21.2973 5.785503 12 41
. margins, dydx(weight) at(mpg=(12(2)41))接下来，使用 marginsplot 命令进行绘图：
图8 mpg取值不同时，weight对price的平均边际效应
计算结果表明：当 mpg 小于28时，随着 mpg 增加，weight 对 price 的边际效应逐渐减小且在统计上显著；当 mpg 大于28时，weight 对 price 的边际效应逐渐减小在统计上不显著。
3.3 非线性模型案例 (Logit Model)
由于在非线性模型的回归结果中，例如 Logit Model，自变量的系数值不能直接代表该变量对因变量的边际效应值，因此，我们需要借助 margins 命令来计算边际效应。
使用 Stata 的自带数据 auto.dta (1978年美国汽车数据)，以 foreign (是否进口车) 作为被解释变量、以 mpg (每加仑汽油能够行驶的英里数)、weight (汽车重量) 作为解释变量建立 Logit 模型。Stata 中的回归命令和结果如下所示：
. sysuse "auto.dta", clear. logit foreign mpg weight, nologLogistic regression Number of obs = 74
LR chi2(2) = 35.72 Prob > chi2 = 0.0000Log likelihood = -27.175156 Pseudo R2 = 0.3966
---------------------------------------------------------------
foreign | Coef. Std. Err. z P>|z| [95% Conf. Interval]
---------+-----------------------------------------------------
mpg | -0.169 0.092 -1.83 0.067 -0.349 0.012
weight | -0.004 0.001 -3.86 0.000 -0.006 -0.002
_cons | 13.708 4.519 3.03 0.002 4.852 22.565
---------------------------------------------------------------回归结果显示 mpg 与 weight 的系数值显著为负，表明当车辆的 mpg 与 weight 增加时，该车辆是进口车的概率减小。但从这两个变量的系数值无法直接看出 mpg 与 weight 对 车辆是进口车的概率 的边际效应。于是，我们可以使用 margins 命令附加 dydx 选项来进行计算。Stata 中的命令和结果如下所示：
Average marginal effects
Expression : Pr(foreign), predict()
dy/dx w.r.t. : mpg
----------------------------------------------
| Delta-method
| dy/dx Std. Err. z P>|z|
-----+----------------------------------------
mpg | -0.020 0.010 -2.03 0.042
----------------------------------------------
. margins, dydx(weight)
(output omitted)
----------------------------------------------
| Delta-method
| dy/dx Std. Err. z P>|z|
--------+-------------------------------------
weight | -0.000 0.000 -8.01 0.000
----------------------------------------------结果表明：当 mpg 增加 1 个单位时，车辆为进口车的概率减少 1.97%；当 weight 增加 1个单位时，车辆为进口车的概率减少 0.04%。
4. 相关阅读资料
Marginal means, adjusted predictions, and marginal effects
The Stata Blog » margins
The Stata Blog » marginal effects
Margins plots
Trouble at the margins
In the spotlight: Estimating, graphing, and interpreting interactions using margins
In the spotlight: Visualizing continuous-by-continuous interactions with margins and twoway contour
Effects of nonlinear models with interactions of discrete and continuous variables: Estimating, graphing, and interpreting
Ben Jann, 2013, PPT, Predictive Margins and Marginal Effects in Stata
特别说明： 文中包含的链接在微信中无法生效。请点击本文底部左下角的【阅读原文】，转入本文【简书版】。
关于我们
【Stata 连享会(公众号：StataChina)】由中山大学连玉君老师团队创办，旨在定期与大家分享 Stata 应用的各种经验和技巧。
公众号推文同步发布于 【简书-Stata连享会】 和 【知乎-连玉君Stata专栏】。可以在简书和知乎中搜索关键词
Stata或Stata连享会后关注我们。点击推文底部【阅读原文】可以查看推文中的链接并下载相关资料。
联系我们
欢迎赐稿： 欢迎将您的文章或笔记投稿至
Stata连享会(公众号: StataChina)，我们会保留您的署名；录用稿件达五篇以上，即可免费获得 Stata 现场培训 (初级或高级选其一) 资格。意见和资料： 欢迎您的宝贵意见，您也可以来信索取推文中提及的程序和数据。
招募英才： 欢迎加入我们的团队，一起学习 Stata。合作编辑或撰写稿件五篇以上，即可免费获得 Stata 现场培训 (初级或高级选其一) 资格。
联系邮件： StataChina@163.com
往期精彩推文
Stata连享会推文列表