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

作者: stata连享会 | 来源:发表于2019-06-09 07:07 被阅读0次

Stata: 手动计算和图示边际效应
Stata：图示连续变量的边际效应(交乘项)
Stata边际效应输出
Stata：边际效应分析
Corporate Finance (4th): Chapter
Stata：多期倍分法 (DID) 图示
Stata：图示交互效应\调节效应
Stata 概率模拟
stata中设置滞后值时需要注意的问题
建筑工程工程量计算——混凝土工程量计算（附图解及计算公式）

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

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

2020寒假Stata现场班 (北京, 1月8-17日，连玉君-江艇主讲)，「+助教招聘」

2020寒假Stata现场班

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

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

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

1. 边际效应计算原理

1.1 两个变量的交乘项

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

$Y=\beta_{0}+\beta_{1} X+\beta_{2} Z+\beta_{3} X Z + u$

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

$\frac{\partial Y}{\partial X}=\beta_{1}+\beta_{3} Z$

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

例如，假设 $\hat{\beta}_{1}=2$ ， $\hat{\beta}_{3}=0.5$ ，则当 $Z=1$ 时，「 $X$ 对 $Y$ 的边际影响」为：

$\frac{\partial Y}{\partial X}|_{ Z=1}=\beta_{1}+\beta_{3} Z = 2+0.5\times 1 = 2.5$

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

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

$\operatorname{Var}({\partial Y}/{\partial X}) = \sqrt{\operatorname{Var}\left(\beta_{1}\right)+Z^{2} \operatorname{Var}\left(\beta_{3}\right)+2 Z \operatorname{Cov}\left(\beta_{1}, \beta_{3}\right)}$

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

1.2 三个变量的交乘项

对于含有三个变量交乘项的模型（如下模型），也可以用上述原理来手动计算变量的边际效应。
$\begin{array} {c}{ Y=\beta_{0}+\beta_{1} X+\beta_{2} Z+\beta_{3} W} \\ {\qquad+\beta_{4} X Z+\beta_{5} X W+\beta_{6} Z W} \\ {+\beta_{7} X Z W+u} \end{array}$

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

$\begin{aligned} {\partial Y}/{\partial X}=\beta_{1}+\beta_{4} Z+\beta_{5} W +\beta_{7} Z W \end{aligned}$

$\begin{align} \operatorname{Var}({\partial Y}/{\partial X}) &=\operatorname{Var}\left(\beta_{1}\right)+Z^{2} \operatorname{Var}\left(\beta_{4}\right)+W^{2} 2 \operatorname{Var}\left(\beta_{5}\right) +Z^{2} W^{2} \operatorname{Var}\left(\beta_{7}\right) \\ &+2 Z \operatorname{cov}\left(\beta_{1} \beta_{4}\right) +2 W \operatorname{cov}\left(\beta_{1} \beta_{5}\right) +2 Z W \operatorname{cov}\left(\beta_{1} \beta_{7}\right) \\ &+2 Z W \operatorname{cov}\left(\beta_{4} \beta_{5}\right) +2 W Z^{2} \operatorname{cov}\left(\beta_{4} \beta_{7}\right) +2 Z W^{2} \operatorname{cov}\left(\beta_{5} \beta_{7}\right) \end{align}$

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

$Y=\beta_{0}+\beta_{1} X+\beta_{2} W+\beta_{3} D+\beta_{4} XW+\beta_{5} XD+\beta_{6} WD+\beta_{7} XWD+u$

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

$Y=\beta_{0}+\beta_{1g} X+\beta_{2g} W+\beta_{3g} XW+u, \ \ g=1,2$

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

连享会计量方法专题……

2. Stata 实现

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

2.1 两个变量的交乘项

$wage=\beta_{0}+\beta_{1} grade+\beta_{2} ttl +\beta_{3} grade\times ttl+u$

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

计算边际效应

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

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

sysuse "nlsw88.dta", clear
gen grade_x_ttl = grade*ttl
reg wage grade ttl grade_x_ttl 

      Source |       SS           df       MS      Number of obs   =     2,244
-------------+----------------------------------   F(3, 2240)      =    133.83
       Model |  11301.2662         3  3767.08872   Prob > F        =    0.0000
    Residual |  63053.0643     2,240  28.1486894   R-squared       =    0.1520
-------------+----------------------------------   Adj R-squared   =    0.1509
       Total |  74354.3305     2,243  33.1495009   Root MSE        =    5.3055

------------------------------------------------------------------------------
        wage |      Coef.   Std. Err.      t    P>|t|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
       grade |      0.252      0.132     1.91   0.056       -0.006       0.509
     ttl_exp |     -0.144      0.128    -1.12   0.264       -0.396       0.108
 grade_x_ttl |      0.032      0.010     3.21   0.001        0.013       0.052
       _cons |      0.934      1.658     0.56   0.573       -2.317       4.184
------------------------------------------------------------------------------

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

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

     matrices:
                  e(b) :  1 x 4
                  e(V) :  4 x 4

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

symmetric e(V)[4,4]
                   grade      ttl_exp  grade_x_ttl        _cons
      grade     .0173019
    ttl_exp    .01534545    .01651049
grade_x_ttl   -.00123247   -.00125844    .00009963
      _cons   -.21378964   -.19844023    .01533119    2.7477949

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

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

matrix  b = e(b)
matrix  V = e(V)
scalar b1 = b[1,1]
scalar b3 = b[1,3]
scalar varb1  = V[1,1]
scalar varb3  = V[3,3]
scalar covb1b3= V[1,3]

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

sum ttl
scalar ttl_mean = r(mean)
dis "Margins：" b1+b3*ttl_mean
dis "Std.Err：" sqrt(varb1+varb3*(ttl_mean^2)+2*covb1b3*ttl_mean)

Margins：.65359238
Std.Err：.04536131

绘制边际效应图

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

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

set obs 2875
generate MVZ=(_n+15)/100

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

gen conbx=b1+b3*MVZ 
gen consx=sqrt(varb1+varb3*(MVZ^2)+2*covb1b3*MVZ)
gen ax=1.96*consx
gen upperx=conbx+ax
gen lowerx=conbx-ax

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

gen where=-0.045
gen pipe = "|"
gen yline=0

graph twoway hist ttl, width(0.5) percent color(gs14) yaxis(2)        ///
        ||   scatter where ttl , plotr(m(b 4)) ms(none) mlabcolor(gs5) mlabel(pipe) mlabpos(6) legend(off)   ///
        ||   line conbx  MVZ,clpattern(solid) clwidth(medium) clcolor(black) yaxis(1)  ///
        ||   line upperx MVZ,clpattern(dash)  clwidth(thin)   clcolor(black)  ///
        ||   line lowerx MVZ,clpattern(dash)  clwidth(thin)   clcolor(black)  ///
        ||   line yline  MVZ,clpattern(solid) clwidth(thin)   clcolor(black)   ///
        ||   ,    ///
             xlabel( 0 5 10 15 20 25 30, nogrid labsize(2))  ///
             ylabel(-.2 0 .2 .4 .6 0.8 1 1.2 1.4 1.6 , axis(1) nogrid labsize(2))  ///
             ylabel(0 1 2 3 4 5, axis(2) nogrid labsize(2))  ///
             yscale(noline alt)  ///
             yscale(noline alt axis(2))  ///    
             xscale(noline)  ///
             legend(off)  ///
             xtitle("" , size(2.5))  ///
             ytitle("" , axis(2) size(2.5))  ///
             xsca(titlegap(2))  ///
             ysca(titlegap(2))  ///
             scheme(s2mono) graphregion(fcolor(white) ilcolor(white) lcolor(white))

image.png

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

$\begin{align} wage &= \beta_{0}+\beta_{1} grade +\beta_{2} ttl +\beta_{3} union \\ &\quad +\beta_{4} grade \times ttl +\beta_{5} grade \times union +\beta_{6} ttl \times union \\ & \quad +\beta_{7} grade \times ttl \times union + u \end{align}$

其中，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 将模型二改写成如下形式：

$wage=\beta_{0}+\beta_{1g} grade+\beta_{2g} ttl_{-} \exp +\beta_{3g} grade _{-} ttl+ u$

$\text { if } union=1, g=1 ; \text { if }union=0, g=2$

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

**计算边际效应：union==1**
reg wage grade ttl grade_x_ttl if union==1
matrix b_1=e(b)
matrix V_1=e(V)
scalar b1_1=b_1[1,1]
scalar b3_1=b_1[1,3]
scalar varb1_1=V_1[1,1]
scalar varb3_1=V_1[3,3]
scalar covb1b3_1=V_1[1,3]

sum ttl if union==1
scalar ttl_mean_1=r(mean)
dis "Margins_1：" b1_1+b3_1*ttl_mean_1
dis "Std.Err_1：" sqrt(varb1_1+varb3_1*(ttl_mean_1^2)+2*covb1b3_1*ttl_mean_1) 

gen conbx_1=b1_1+b3_1*MVZ 
gen consx_1=sqrt(varb1_1+varb3_1*(MVZ^2)+2*covb1b3_1*MVZ)
gen ax_1=1.96*consx_1
gen upperx_1=conbx_1+ax_1
gen lowerx_1=conbx_1-ax_1

**计算边际效应：union==0**
reg wage grade ttl grade_x_ttl if union==0
matrix b_0=e(b)
matrix V_0=e(V)
scalar b1_0=b_0[1,1]
scalar b3_0=b_0[1,3]
scalar varb1_0=V_0[1,1]
scalar varb3_0=V_0[3,3]
scalar covb1b3_0=V_0[1,3]

sum ttl if union==0
scalar ttl_mean_0=r(mean)
dis "Margins_0：" b1_0+b3_0*ttl_mean_0
dis "Std.Err_0：" sqrt(varb1_0+varb3_0*(ttl_mean_0^2)+2*covb1b3_0*ttl_mean_0) 

gen conbx_0=b1_0+b3_0*MVZ 
gen consx_0=sqrt(varb1_0+varb3_0*(MVZ^2)+2*covb1b3_0*MVZ)
gen ax_0=1.96*consx_0
gen upperx_0=conbx_0+ax_0
gen lowerx_0=conbx_0-ax_0

**绘图**
graph twoway line conbx_1  MVZ, clpattern(solid) clwidth(medium) clcolor(blue)  ///
        ||   line upperx_1 MVZ, clpattern(dash)  clwidth(thin)   clcolor(blue)  ///
        ||   line lowerx_1 MVZ, clpattern(dash)  clwidth(thin)   clcolor(blue)  ///
        ||   line conbx_0  MVZ, clpattern(solid) clwidth(medium) clcolor(red)   ///
        ||   line upperx_0 MVZ, clpattern(dash)  clwidth(thin)   clcolor(red)   ///
        ||   line lowerx_0 MVZ, clpattern(dash)  clwidth(thin)   clcolor(red)   ///
        ||   line yline    MVZ, clpattern(solid) clwidth(thin)   clcolor(black)    ///
        ||   ,    ///
             xlabel( 0 5 10 15 20 25 30, nogrid labsize(2))  ///
             ylabel(-.2 0 .2 .4 .6 0.8 1 1.2 1.4 1.6 , nogrid labsize(2))  ///
             xscale(noline) ///
             legend(off)  ///
             xtitle("" , size(2.5))  ///
             xsca(titlegap(2))  ///
             ysca(titlegap(2))  ///
             scheme(s2mono) graphregion(fcolor(white) ilcolor(white) lcolor(white))

image.png

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

连享会计量方法专题……

关于我们

Stata连享会 由中山大学连玉君老师团队创办，定期分享实证分析经验。
推文同步发布于 CSDN 、简书和知乎Stata专栏。可在百度中搜索关键词「Stata连享会」查看往期推文。
点击推文底部【阅读原文】可以查看推文中的链接并下载相关资料。
欢迎赐稿： 欢迎赐稿。录用稿件达三篇以上，即可免费获得一期 Stata 现场培训资格。
E-mail： StataChina@163.com
往期精彩推文：一网打尽

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

Stata: 手动计算和图示边际效应
作者：高娜娜（中南财经政法大学）连享会-交乘项 (调节效应) 专题系列 Stata：交叉项\交乘项该这么分析...
Stata：图示连续变量的边际效应(交乘项)
作者：连玉君 (知乎 | 简书 | 码云) …… Stata连享会 - 往期精彩推文 …… Stata 现场培...
Stata边际效应输出
关于OLS回归结果的输出，Stata：毕业论文大礼包 A——实证结果输出命令大比拼已经有非常详细的讲解。本文主要简...
Stata：边际效应分析
作者：连玉君 | 杨柳 ( 知乎 | 简书 | 码云 ) Stata连享会精彩推文1 || 精彩推文2 1. ...
Corporate Finance (4th): Chapter
Stata 现场培训报名中下面是公司财务课程教材中部分范例的 Stata 实现代码，以便大家对相关的概念和计算方...
Stata：多期倍分法 (DID) 图示
Stata连享会 (知乎 | 简书 | 码云) Stata 现场培训报名中 Source: mostly-harm...
Stata：图示交互效应\调节效应
作者：邓浩然 (上海财经大学)Stata 连享会：知乎 | 简书 | 码云 Stata连享会精品专题 |...
Stata 概率模拟
Stata做随机模型和Matlab类似的规则，采用分布类型和计算内容结合的方式：概率函数类型： norm ：标准...
stata中设置滞后值时需要注意的问题
在计算时间序列问题的时候，经常需要计算滞后数值。Stata 里使用 lag 函数计算，试例如下：数据库中有变量 ...
建筑工程工程量计算——混凝土工程量计算（附图解及计算公式）
一、基础混凝土计算 1．现浇混凝土基础工程量计算规则：按设计图示尺寸以体积计算。不扣除构件内钢筋、预埋铁件和伸入...