三、样本选择模型的Stata实操
Stata关于样本选择模型的官方命令是heckman，该命令能够进行Heckman两步估计以及模型整体参数的极大似然估计。
3.1 语法介绍
关于heckman更多语法细节，在命令窗口键入如下代码即可了解。
help heckman
常用的语法如下。
**# 一、样本选择模型
**# 1.1 语法介绍
/*
a. 两步估计法
heckman y x1 x2, select( x1 x2 z1) twostep mills(newname)
等价于
heckman y x1 x2, select(y_dummy = x1 x2 z1) twostep
b. MLE
heckman y x1 x2, select(y_dummy = x1 x2 z1) nolog
+ robust/cluster SE
heckman y x1 x2, select(y_dummy = x1 x2 z1) nolog vce(cluster varname)
*/
语法书写大致包括两部分，选择项之前的是第二阶段回归的被解释变量y以及控制变量x1和x2，不需要写入IMR。选择项则包括对选择方程的具体设置，以及其他的细节。
-
select()表示写入选择方程，括号内的是选择方程的具体变量，有两种写法。- 第一种写法：直接写入控制变量（
x1和x2）和外生变量z1。此时严格要求原回归的被解释变量y存在缺失值，且缺失值不能以0作为标记。 - 第二种写法：首先写入
y是否被观测到的虚拟变量y_dummy（y存在缺失值的样本y_dummy标记为0，不存在缺失值的样本标记为1，且y取值为0不作为缺失值处理）作为选择方程的被解释变量，等号后面跟着控制变量和外生变量。若选择这种写法，则要求提前生成y_dummy。
- 第一种写法：直接写入控制变量（
-
twostep表示使用两步估计法，默认使用MLE。 -
mills()表示生成各个样本的IMR，并以newname作为变量名。nshazard()的作用相同。 -
nolog表示在使用MLE时回归结果不显示迭代过程。 -
vce()表示使用稳健标准误，括号内填入robust表示使用异方差稳健标准误；填入cluster varname表示使用聚类稳健标准误，以变量varname作为聚类标准，根据经验法则，要求聚类数目大于30。需要注意的是，vce()不能在两步法中使用。
3.2 实例操作
heckman的help文件附带演示案例，下面根据演示案例中的数据和代码，对样本选择模型在Stata中的具体实现进行介绍。
首先调用数据集womenwk，由于该数据集并非Stata自带，因此可能无法调用成功，这种情况下可以直接调用本次推文提供的数据。
*- Stata Version: 16 | 17
*- 定义路径
cd "C:\Users\KEMOSABE\Desktop\heckman"
**# 1.2 实例操作
webuse womenwk, clear
/*
*- 或者直接调用已下载的数据集
use womenwk.dta, clear
*/
该例中，基准回归的被解释变量是wage，解释变量是educ和age；选择方程中额外引入了两个外生解释变量married和children。
由于被解释变量wage存在缺失值，因此我们怀疑基准回归方程（OLS）可能存在样本选择偏差，因此选用样本选择模型做进一步的稳健性检验。
该数据集为截面数据，因此标准误无法聚类到个体层面，而同一地区的个体特征可能存在一定的相关性，因此我们尝试将标准误聚类到county层面，下面判断是否能够使用聚类标准误。
tabulate county
/*
county of |
residence | Freq. Percent Cum.
------------+-----------------------------------
0 | 200 10.00 10.00
1 | 200 10.00 20.00
2 | 200 10.00 30.00
3 | 200 10.00 40.00
4 | 200 10.00 50.00
5 | 200 10.00 60.00
6 | 200 10.00 70.00
7 | 200 10.00 80.00
8 | 200 10.00 90.00
9 | 200 10.00 100.00
------------+-----------------------------------
Total | 2,000 100.00
*/
可以看到，地区数目只有10个，远无法满足聚类数目最低值（30个）的要求，因此选用异方差稳健标准误。
根据wage是否存在缺失值生成虚拟变量wage_dummy。可以看到，2,000个样本中，wage数据缺失的有657个。
*- 查看y的缺失值情况
nmissing wage
/*
wage 657
*/
*- 生成y是否取值的虚拟变量y_dummy
gen wage_dummy = (wage != .)
然后是定义控制变量和外生变量的全局暂元。
*- 定义全局暂元
global xlist educ age
global zlist married children
第一个是基准回归，可以看到：
- 参与回归的样本数目为1,343个，即
wage存在缺失值的样本（657个）在回归时直接被drop掉。 - 基准回归中两个解释变量的系数均显著为正，模型拟合程度也较好（修正R方为0.2524）。
*-（一）基准OLS
reg wage $xlist
est store OLS
/*
Source | SS df MS Number of obs = 1,343
-------------+---------------------------------- F(2, 1340) = 227.49
Model | 13524.0337 2 6762.01687 Prob > F = 0.0000
Residual | 39830.8609 1,340 29.7245231 R-squared = 0.2535
-------------+---------------------------------- Adj R-squared = 0.2524
Total | 53354.8946 1,342 39.7577456 Root MSE = 5.452
------------------------------------------------------------------------------
wage | Coef. Std. Err. t P>|t| [95% Conf. Interval]
-------------+----------------------------------------------------------------
education | .8965829 .0498061 18.00 0.000 .7988765 .9942893
age | .1465739 .0187135 7.83 0.000 .109863 .1832848
_cons | 6.084875 .8896182 6.84 0.000 4.339679 7.830071
------------------------------------------------------------------------------
*/
第二个是样本选择模型，使用MLE方法进行估计，可以看到：
- 选择方程中两个外生变量均显著为正，说明外生变量的选择是有效的。
- 在第二阶段回归中，
IMR（即lambda）的估计系数为4.2244，但显著性未知，该值等于rho和sigma的乘积，其中：-
sigma是原方程干扰项的标准差。 -
rho是选择方程干扰项和第二阶段回归干扰项的相关系数。如果rho等于0，表示第二阶段回归中IMR的系数不显著，说明样本选择偏差在原方程中不怎么严重，反之则需要考虑样本选择偏差带来的估计偏误。回归结果的末尾是LR检验，检验的原假设是H0: rho = 0，p值说明至少可以在1%的水平下拒绝原假设，可以认为rho显著不等于0，这说明原模型中确实存在严重的样本选择偏差，基准回归结果不可信。
-
- 第二阶段回归结果中，两个解释变量仍旧显著为正，且相较于基准回归结果取值变化不大，说明考虑到样本选择偏差后基准回归结果依然是稳健的。
- 回归结果右上角说明总的样本量是2,000个，参与选择方程回归的样本为2,000个，参与第二阶段回归的样本为1,343个。
*-（二）MLE
heckman wage $xlist , select( wage_dummy = $xlist $zlist ) nolog mills(imr1)
est store MLE
/*
Heckman selection model Number of obs = 2,000
(regression model with sample selection) Selected = 1,343
Nonselected = 657
Wald chi2(2) = 508.44
Log likelihood = -5178.304 Prob > chi2 = 0.0000
------------------------------------------------------------------------------
| Coef. Std. Err. z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
wage |
education | .9899537 .0532565 18.59 0.000 .8855729 1.094334
age | .2131294 .0206031 10.34 0.000 .1727481 .2535108
_cons | .4857752 1.077037 0.45 0.652 -1.625179 2.59673
-------------+----------------------------------------------------------------
wage_dummy |
education | .0557318 .0107349 5.19 0.000 .0346917 .0767718
age | .0365098 .0041533 8.79 0.000 .0283694 .0446502
married | .4451721 .0673954 6.61 0.000 .3130794 .5772647
children | .4387068 .0277828 15.79 0.000 .3842534 .4931601
_cons | -2.491015 .1893402 -13.16 0.000 -2.862115 -2.119915
-------------+----------------------------------------------------------------
/athrho | .8742086 .1014225 8.62 0.000 .6754241 1.072993
/lnsigma | 1.792559 .027598 64.95 0.000 1.738468 1.84665
-------------+----------------------------------------------------------------
rho | .7035061 .0512264 .5885365 .7905862
sigma | 6.004797 .1657202 5.68862 6.338548
lambda | 4.224412 .3992265 3.441942 5.006881
------------------------------------------------------------------------------
LR test of indep. eqns. (rho = 0): chi2(1) = 61.20 Prob > chi2 = 0.0000
*/
第三个同样是样本选择模型，使用MLE方法进行估计，这里将使用异方差稳健标准误。可以看到，使用稳健标准误的回归结果与上文使用普通标准误的回归结果基本保持一致。需要注意的是，在使用稳健标准误后，对H0的统计检验将变成Wald检验，这里Wald检验同样拒绝原假设。
*-（三）MLE + robust SE
heckman wage $xlist , select( wage_dummy = $xlist $zlist ) nolog vce(robust)
est store MLE_r
/*
Heckman selection model Number of obs = 2,000
(regression model with sample selection) Selected = 1,343
Nonselected = 657
Wald chi2(2) = 497.82
Log pseudolikelihood = -5178.304 Prob > chi2 = 0.0000
------------------------------------------------------------------------------
| Robust
| Coef. Std. Err. z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
wage |
education | .9899537 .0534141 18.53 0.000 .8852641 1.094643
age | .2131294 .0211095 10.10 0.000 .1717555 .2545034
_cons | .4857752 1.099121 0.44 0.659 -1.668462 2.640013
-------------+----------------------------------------------------------------
wage_dummy |
education | .0557318 .0108899 5.12 0.000 .034388 .0770755
age | .0365098 .0042243 8.64 0.000 .0282303 .0447893
married | .4451721 .0668243 6.66 0.000 .3141988 .5761453
children | .4387068 .0272779 16.08 0.000 .385243 .4921705
_cons | -2.491015 .1884227 -13.22 0.000 -2.860317 -2.121713
-------------+----------------------------------------------------------------
/athrho | .8742086 .1051331 8.32 0.000 .6681514 1.080266
/lnsigma | 1.792559 .0288316 62.17 0.000 1.73605 1.849068
-------------+----------------------------------------------------------------
rho | .7035061 .0531006 .5837626 .7932976
sigma | 6.004797 .1731281 5.674882 6.353893
lambda | 4.224412 .4172197 3.406676 5.042147
------------------------------------------------------------------------------
Wald test of indep. eqns. (rho = 0): chi2(1) = 69.14 Prob > chi2 = 0.0000
*/
第四个是样本选择模型，使用两步法进行估计，可以看到：
- 选择方程中两个工具变量的引入是有效的。
- 第二阶段回归中，
IMR的回归系数等于4.0016，与MLE方法下的4.2244相差不大，但两步法下IMR回归系数可以直接进行z检验，并且统计结果说明IMR回归系数至少在1%的水平下显著为正，这同时说明原方程中的样本选择偏差问题不可忽视。 - 第二阶段回归结果中，两个解释变量仍旧显著为正，且大小与基准回归结果相比变化不大，这说明在考虑样本选择偏差的情况下，基准回归结果是可信的。
- 右上角同样说明在两步估计法下，参与选择方程的样本数是2,000个，但参与第二阶段回归的样本数只有1,343个，减少的657个样本与被解释变量
wage存在缺失值的样本完全重合。
*-（四）两步估计法
heckman wage $xlist , select( wage_dummy = $xlist $zlist ) twostep mills(imr2)
est store TwoStep
/*
Heckman selection model -- two-step estimates Number of obs = 2,000
(regression model with sample selection) Selected = 1,343
Nonselected = 657
Wald chi2(2) = 442.54
Prob > chi2 = 0.0000
------------------------------------------------------------------------------
| Coef. Std. Err. z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
wage |
education | .9825259 .0538821 18.23 0.000 .8769189 1.088133
age | .2118695 .0220511 9.61 0.000 .1686502 .2550888
_cons | .7340391 1.248331 0.59 0.557 -1.712645 3.180723
-------------+----------------------------------------------------------------
wage_dummy |
education | .0583645 .0109742 5.32 0.000 .0368555 .0798735
age | .0347211 .0042293 8.21 0.000 .0264318 .0430105
married | .4308575 .074208 5.81 0.000 .2854125 .5763025
children | .4473249 .0287417 15.56 0.000 .3909922 .5036576
_cons | -2.467365 .1925635 -12.81 0.000 -2.844782 -2.089948
-------------+----------------------------------------------------------------
/mills |
lambda | 4.001615 .6065388 6.60 0.000 2.812821 5.19041
-------------+----------------------------------------------------------------
rho | 0.67284
sigma | 5.9473529
------------------------------------------------------------------------------
*/
第五个是手工完成两步估计法，但实际使用时这种方法十分不推荐（包括上文使用两步法估计的样本选择模型），因为两步估计法下第一步回归的估计误差会带到第二步，造成效率损失，影响第二步估计系数的显著性。
手工两步估计法的思路是先用控制变量和外生变量对wage_dummy做probit回归，计算出拟合值y_hat，然后根据y_hat计算出全部样本的IMR，最后将IMR作为额外的控制变量引入基准回归方程中做进一步的回归。
第一阶段（选择方程）回归的代码与结果如下。
*-（五）手工两步估计法
*- 第一阶段方程（选择方程）回归
probit wage_dummy $xlist $zlist , nolog
est store First
/*
Probit regression Number of obs = 2,000
LR chi2(4) = 478.32
Prob > chi2 = 0.0000
Log likelihood = -1027.0616 Pseudo R2 = 0.1889
------------------------------------------------------------------------------
wage_dummy | Coef. Std. Err. z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
education | .0583645 .0109742 5.32 0.000 .0368555 .0798735
age | .0347211 .0042293 8.21 0.000 .0264318 .0430105
married | .4308575 .074208 5.81 0.000 .2854125 .5763025
children | .4473249 .0287417 15.56 0.000 .3909922 .5036576
_cons | -2.467365 .1925635 -12.81 0.000 -2.844782 -2.089948
------------------------------------------------------------------------------
*/
可以看到，实际参与回归的样本为2,000个，控制变量和两个外生变量均显著为正，这说明外生变量的选择是有效的。此外，衡量模型整体拟合程度的统计指标是伪R方（Pseudo R2），伪R方为0.1889，整体拟合程度较好。
下面根据拟合值y_hat计算出额外的控制变量imr。其中，normalden()是标准正态的概率密度函数pdf，normal()是标准正态的累积分布函数cdf。
*- 计算IMR
predict y_hat, xb
gen imr = normalden(y_hat) / normal(y_hat)
获得各个样本的imr后，就可以进行第二阶段回归。可以发现：
- 与样本选择模型的两步估计法结果相比，手工两步法估计结果在系数值大小方面没有任何改变，在系数标准误方面变化也不大，从而各个变量的系数显著性保持高度一致。
- 这提示我们，两步估计法（包括手工两步法）带来的效率损失不可忽视，虽然本例中的结果和MLE以及基准回归基本保持一致，实际使用应尽量避免使用这种方法，除非能够找到比较好的方法来校正误差。
- 第二阶段实际参与回归的样本是1,343个。
- 衡量模型整体拟合程度的统计指标是修正R方，修正R方等于0.2777，整体拟合程度较高。事实上，由于该例并未规定原方程的核心解释变量，因此对模型整体拟合程度的考察是有意义的，但在现实应用中，R方即便很小问题也不大，只要求核心解释变量的系数大小和显著性符合预期即可。就算对R方有具体的要求，只要我们在方程中引入足够多的控制变量，控制足够多的“固定效应”，R方总能达到预定的要求。因此，实际应用中R方意义不大。
*- 第二阶段方程回归
reg wage $xlist imr
est store Second
/*
Source | SS df MS Number of obs = 1,343
-------------+---------------------------------- F(3, 1339) = 173.01
Model | 14904.6806 3 4968.22688 Prob > F = 0.0000
Residual | 38450.214 1,339 28.7156191 R-squared = 0.2793
-------------+---------------------------------- Adj R-squared = 0.2777
Total | 53354.8946 1,342 39.7577456 Root MSE = 5.3587
------------------------------------------------------------------------------
wage | Coef. Std. Err. t P>|t| [95% Conf. Interval]
-------------+----------------------------------------------------------------
education | .9825259 .0504982 19.46 0.000 .8834616 1.08159
age | .2118695 .0206636 10.25 0.000 .171333 .252406
imr | 4.001616 .5771027 6.93 0.000 2.869492 5.133739
_cons | .7340391 1.166214 0.63 0.529 -1.553766 3.021844
------------------------------------------------------------------------------
*/
可以对比一下用各个方法计算的IMR。其中，imr1是用样本选择模型的MLE方法计算的，imr2是用样本选择模型的两步估计法计算的，imr是用手工两步法计算的。打开Stata的数据编辑器，可以发现：
- 所有样本的
IMR均有取值。 -
imr和imr2保持高度一致（imr四舍五入了），因为两者都是用两步估计法计算出来的。 -
imr1和imr2存在差距，但差别不大。
*- 对比imr、imr1和imr2
order imr1 imr2 imr wage_dummy
最后是五个模型回归结果的对比，可以发现：
- 两个解释变量（
education和age）的回归系数的符号、大小、显著性在模型中基本保持一致，并且与模型
的区别不大，这说明在考虑到样本选择偏差后，根据基本回归结果得出的结论依然成立。
- 模型
缺少变量
lambda，在实际汇报时应该根据原始回归表进行汇报。 - 模型
汇报的样本量是2,000个全样本，但应该清楚不同回归阶段实际参与的样本数目是不同的。
- 关于模型拟合程度的统计量，模型
汇报的都是修正R方，模型
的是伪R方，模型
Wald chi2，因此，在实际汇报时，应该把Wald chi2作为这三个模型整体拟合度的汇报指标。
*- 回归结果输出
local mlist "OLS MLE MLE_r TwoStep First Second"
esttab `mlist' using 样本选择模型.rtf, replace ///
b(%6.4f) t(%6.4f) ///
scalar(N r2_a r2_p) compress nogap ///
star(* 0.1 ** 0.05 *** 0.01) ///
mtitle(`mlist') ///
title("Sample Selection Model")
/*
Sample Selection Model
----------------------------------------------------------------------------------------
(1) (2) (3) (4) (5) (6)
OLS MLE MLE_r TwoStep First Second
----------------------------------------------------------------------------------------
main
education 0.8966*** 0.9900*** 0.9900*** 0.9825*** 0.0584*** 0.9825***
(18.0015) (18.5884) (18.5336) (18.2347) (5.3184) (19.4566)
age 0.1466*** 0.2131*** 0.2131*** 0.2119*** 0.0347*** 0.2119***
(7.8325) (10.3445) (10.0964) (9.6081) (8.2096) (10.2533)
married 0.4309***
(5.8061)
children 0.4473***
(15.5636)
imr 4.0016***
(6.9340)
_cons 6.0849*** 0.4858 0.4858 0.7340 -2.4674*** 0.7340
(6.8399) (0.4510) (0.4420) (0.5880) (-12.8133) (0.6294)
----------------------------------------------------------------------------------------
wage_dummy
education 0.0557*** 0.0557*** 0.0584***
(5.1916) (5.1178) (5.3184)
age 0.0365*** 0.0365*** 0.0347***
(8.7905) (8.6428) (8.2096)
married 0.4452*** 0.4452*** 0.4309***
(6.6054) (6.6618) (5.8061)
children 0.4387*** 0.4387*** 0.4473***
(15.7906) (16.0828) (15.5636)
_cons -2.4910*** -2.4910*** -2.4674***
(-13.1563) (-13.2204) (-12.8133)
----------------------------------------------------------------------------------------
/
athrho 0.8742*** 0.8742***
(8.6195) (8.3153)
lnsigma 1.7926*** 1.7926***
(64.9526) (62.1733)
----------------------------------------------------------------------------------------
/mills
lambda 4.0016***
(6.5975)
----------------------------------------------------------------------------------------
N 1343 2000 2000 2000 2000 1343
r2_a 0.2524 0.2777
r2_p 0.1889
----------------------------------------------------------------------------------------
t statistics in parentheses
* p<0.1, ** p<0.05, *** p<0.01
*/
以上命令、代码和实例均基于截面数据，现阶段的文献大多数也是将面板数据转为截面数据处理，最多控制几个个体虚拟变量和时间虚拟变量，并将其视为面板数据回归。
事实上，Stata在16及以上版本中提供了估计面板数据样本选择模型官方命令xtheckman，但这个命令只能估计随机效应模型，估计固定效应模型需安装一个外部命令xtheckmanfe，该命令由Boston College的Fernando编写，可以键入如下代码进行安装或更新。
ssc install xtheckmanfe, replace
[6] Fernando R-A. XTHECKMANFE Stata module to fit panel data models in the presence of endogeneity and selection: Boston College Department of Economics, 2020. [Program]
四、处理效应模型的Stata实操
处理效应模型在Stata中的官方命令是treatreg和etregress，其中treatreg在Stata第14版之后被etregress替代，因此下面所使用的均为etregress。
与heckman类似，etregress同样可以进行处理效应模型的两步法估计和极大似然估计。
4.1 语法介绍
etregress的语法结构和heckman极为相似，但还是存在几点不同，具体来说：
- 第一，选择方程的写法是
treat()，而非select()。 - 第二，选择方程中的被解释变量是原方程中的核心解释变量
D，该解释变量取值为0或1，且不存在缺失值。在etregress语法中，变量D不可省略。 - 第三，选择方程中的工具变量
z1直接影响是变量D，而非样本选择模型中的y_dummy。 - 第四，处理效应模型中生成
IMR的选择项是hazard()，而非样本选择模型中的mills()或nshazad()。 - 第五，虽然第二步回归方程中没有写入
D，但该命令在进行第二部回归时默认自动引入D。因此，如果在语法书写时在第二步回归方程中手动加入D，正式回归时该变量会产生完全的多重共线性而被omitted掉。 - 第六，对于工具变量
z1的选择，还是那句话，尽量避免使用D的滞后项D_lag，如果实在找不到一个良好的工具变量而选择使用D_lag，也应该要确保面板结构中不同个体D取值为1的时点是交错的，这样才能避免产生多重共线性问题。当然，如果是纯粹的截面数据（非混合面板），就不存在所谓的滞后变量，在这种情况下找到一个良好的工具变量是必须的。
**# 二、处理效应模型
**# 2.1 语法介绍
/*
a. 两步估计法
etregress y x1 x2, treat(D = x1 x2 z1) twostep hazard(newname)
b. MLE
etregress y x1 x2, treat(D = x1 x2 z1) nolog
+ robust/cluster SE
etregress y x1 x2, treat(D = x1 x2 z1) nolog vce(cluster varname)
*/
4.2 实例操作
这里使用的是etregress命令的help文件中的示例数据集union3与代码。
同样，如果非内嵌数据集调用失败可以直接调用本次推文提供的数据。
**# 2.2 实例操作
webuse union3, clear
/*
*- 或者直接调用已下载的数据集
use union3.dta, clear
*/
由于union3为截面数据集，因此标准误无法聚类到个体层面，同时考虑到同一行业内不同个体的某些特征存在一致性，因此考虑将标准误聚类到行业层面。
*- 判断是否能够使用聚类标准误
tabulate ind_code
/*
industry of |
employment | Freq. Percent Cum.
------------+-----------------------------------
1 | 13 0.79 0.79
2 | 2 0.12 0.91
3 | 12 0.73 1.64
4 | 331 20.13 21.78
5 | 97 5.90 27.68
6 | 346 21.05 48.72
7 | 158 9.61 58.33
8 | 61 3.71 62.04
9 | 136 8.27 70.32
10 | 13 0.79 71.11
11 | 374 22.75 93.86
12 | 101 6.14 100.00
------------+-----------------------------------
Total | 1,644 100.00
*/
可见，如果将标准误聚类到行业层面，聚类数目少于经验数目（30个），聚类数目过少将导致标准误无法收敛至真实值（该问题在往期推文『双重差分法(DID) | 空间DID』中有讨论），因此该例还是使用异方差稳健标准误。
下面观察核心解释变量union的取值情况。
*- 查看D的取值情况
tabulate union
/*
1 if union | Freq. Percent Cum.
------------+-----------------------------------
0 | 984 79.10 79.10
1 | 260 20.90 100.00
------------+-----------------------------------
Total | 1,244 100.00
*/
nmissing union
/*
union 449
*/
可见，总共有1,693个样本，其中union取值为1的样本有260个，取值为0的样本有984个，缺失值样本有449个。剔除掉缺失值样本，对于剩下的1,244个样本，我们怀疑存在自选择偏差，因此使用处理效应模型来检验模型中是否存在该问题，并与基准回归结果对比进行稳健性检验。
下面将定义变量的全局暂元，并对所使用的变量进行简要说明。
*- 定义全局暂元
global xlist1 black tenure
global xlist2 age grade smsa black tenure
global xlist3 union age grade smsa black tenure
其中，wage是被解释变量，union是核心解释变量，并且我们其怀疑存在自选择问题。age、grade、smsa、black和tenure是基准回归方程的控制变量；black和tenure是选择方程的控制变量，south是选择方程的工具变量。
可以看到，选择方程的控制变量并未与基准回归的控制变量重合，其数目反而少于基准回归控制变量的数目。
第一个是基准回归，可以看到：
- 所有解释变量的系数均至少在1%的水平下显著，特别是我们重点关注的核心解释变量
union。 - 由于样本中各变量存在缺失值，因此最终参与回归的样本只有1,210个。
- 由于各样本
union的取值可能不是随机的，即可能存在自选择问题，导致基准回归结果存在偏误，因此还需要做进一步的稳健性检验才能得出令人信服的结论。
*-（一）基准OLS
reg wage $xlist3
est store OLS
/*
Source | SS df MS Number of obs = 1,210
-------------+---------------------------------- F(6, 1203) = 103.36
Model | 2163.49455 6 360.582425 Prob > F = 0.0000
Residual | 4196.70492 1,203 3.48853277 R-squared = 0.3402
-------------+---------------------------------- Adj R-squared = 0.3369
Total | 6360.19947 1,209 5.26071089 Root MSE = 1.8678
------------------------------------------------------------------------------
wage | Coef. Std. Err. t P>|t| [95% Conf. Interval]
-------------+----------------------------------------------------------------
union | 1.003913 .1338083 7.50 0.000 .7413894 1.266437
age | .1475526 .0195242 7.56 0.000 .1092474 .1858578
grade | .4368545 .0294718 14.82 0.000 .3790327 .4946764
smsa | .9754817 .1252669 7.79 0.000 .7297159 1.221248
black | -.6183346 .1252 -4.94 0.000 -.8639693 -.3726999
tenure | .2118016 .0338481 6.26 0.000 .1453938 .2782094
_cons | -4.326028 .5315474 -8.14 0.000 -5.368891 -3.283165
------------------------------------------------------------------------------
*/
第二个是处理效应模型，使用MLE方法从模型整体角度对参数进行估计。可以看到：
- 第一步回归中
south显著为负，说明工具变量的选择有效。 - 第二步回归中核心解释变量
union显著为正，这一点在处理效应模型回归结果表中显示为1.union，说明相对于union取值为0的样本，union取值为1的样本的平均工资高2.9458个单位。2.9458大约是基准回归结果1.0039的三倍，这说明在未考虑自选择偏差的情况下，低估了union对wage的影响。 - 其余解释变量估计系数的符号和显著性基本没变，但估计系数的数值改变较大，这说明自选择偏差带来的估计偏误不能忽视。
-
IMR（lambda）的估计系数为-1.1603，同时rho的LR检验结果说明至少在1%的水平下拒绝H0。这都说明基准回归模型中确实存在严重的自选择偏差。 - 第一步回归和第二步回归的样本数均为1,210个。
*-（二）MLE
etregress wage $xlist2 , treat(union = $xlist1 south) nolog hazard(imr1)
est store MLE
/*
Linear regression with endogenous treatment Number of obs = 1,210
Estimator: maximum likelihood Wald chi2(6) = 681.89
Log likelihood = -3051.575 Prob > chi2 = 0.0000
------------------------------------------------------------------------------
| Coef. Std. Err. z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
wage |
age | .1487409 .0193291 7.70 0.000 .1108566 .1866252
grade | .4205658 .0293577 14.33 0.000 .3630258 .4781058
smsa | .9117044 .1249041 7.30 0.000 .6668969 1.156512
black | -.7882471 .1367078 -5.77 0.000 -1.05619 -.5203048
tenure | .1524015 .0369596 4.12 0.000 .0799621 .2248409
1.union | 2.945815 .2749621 10.71 0.000 2.4069 3.484731
_cons | -4.351572 .5283952 -8.24 0.000 -5.387208 -3.315936
-------------+----------------------------------------------------------------
union |
black | .4557499 .0958042 4.76 0.000 .2679771 .6435226
tenure | .0871536 .0232483 3.75 0.000 .0415878 .1327195
south | -.5807419 .0851111 -6.82 0.000 -.7475566 -.4139271
_cons | -.8855758 .0724506 -12.22 0.000 -1.027576 -.7435753
-------------+----------------------------------------------------------------
/athrho | -.6544347 .0910314 -7.19 0.000 -.832853 -.4760164
/lnsigma | .7026769 .0293372 23.95 0.000 .645177 .7601767
-------------+----------------------------------------------------------------
rho | -.5746478 .060971 -.682005 -.4430476
sigma | 2.019151 .0592362 1.906325 2.138654
lambda | -1.1603 .1495097 -1.453334 -.8672668
------------------------------------------------------------------------------
LR test of indep. eqns. (rho = 0): chi2(1) = 19.84 Prob > chi2 = 0.0000
*/
第三个是处理效应模型，使用MLE方法并使用异方差稳健标准误。可以看到，与使用普通标准误的MLE方法相比，除系数的z统计量有细微差别，以及对rho的显著性检验使用的是Wald统计指标，其他结果没有实质性差异。
*-（三）MLE + robust SE
etregress wage $xlist2 , treat(union = $xlist1 south) nolog vce(robust)
est store MLE_r
/*
Linear regression with endogenous treatment Number of obs = 1,210
Estimator: maximum likelihood Wald chi2(6) = 487.30
Log pseudolikelihood = -3051.575 Prob > chi2 = 0.0000
------------------------------------------------------------------------------
| Robust
| Coef. Std. Err. z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
wage |
age | .1487409 .0206909 7.19 0.000 .1081875 .1892943
grade | .4205658 .0377922 11.13 0.000 .3464945 .4946371
smsa | .9117044 .1199377 7.60 0.000 .6766308 1.146778
black | -.7882471 .1408988 -5.59 0.000 -1.064404 -.5120906
tenure | .1524015 .0473782 3.22 0.001 .0595419 .2452611
1.union | 2.945815 .5643841 5.22 0.000 1.839643 4.051988
_cons | -4.351572 .6371071 -6.83 0.000 -5.600279 -3.102865
-------------+----------------------------------------------------------------
union |
black | .4557499 .093489 4.87 0.000 .2725148 .6389849
tenure | .0871536 .0251543 3.46 0.001 .0378521 .1364552
south | -.5807419 .0837513 -6.93 0.000 -.7448914 -.4165924
_cons | -.8855758 .0739533 -11.97 0.000 -1.030522 -.7406301
-------------+----------------------------------------------------------------
/athrho | -.6544347 .2322442 -2.82 0.005 -1.109625 -.1992445
/lnsigma | .7026769 .0758152 9.27 0.000 .5540819 .8512719
-------------+----------------------------------------------------------------
rho | -.5746478 .1555525 -.8039298 -.1966492
sigma | 2.019151 .1530822 1.740342 2.342624
lambda | -1.1603 .3823277 -1.909649 -.4109519
------------------------------------------------------------------------------
Wald test of indep. eqns. (rho = 0): chi2(1) = 7.94 Prob > chi2 = 0.0048
*/
第四个是用两步法估计处理效应模型，可以发现：
- 选择方程中的工具变量选择有效。
- 第二阶段回归中的所有解释变量，除
tenure，均至少在1%的水平下显著，lambda显著为负。 - 回归系数的大小与MLE方法下的估计系数相近，但与基准回归相差较大，同时考虑到
IMR显著，这说明原模型中的自选择偏差不可忽视，在未考虑自选择偏差的情况下，将严重低估union对wage的影响。
*-（四）两步估计法
etregress wage $xlist2 , treat(union = $xlist1 south) twostep hazard(imr2)
est store TwoStep
/*
Linear regression with endogenous treatment Number of obs = 1210
Estimator: two-step Wald chi2(8) = 566.56
Prob > chi2 = 0.0000
------------------------------------------------------------------------------
| Coef. Std. Err. z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
wage |
age | .1543231 .0194903 7.92 0.000 .1161227 .1925234
grade | .4225025 .029014 14.56 0.000 .3656362 .4793689
smsa | .8628628 .1285907 6.71 0.000 .6108297 1.114896
black | -.9206944 .1774617 -5.19 0.000 -1.268513 -.572876
tenure | .1003226 .051879 1.93 0.053 -.0013584 .2020037
union | 4.563859 1.006459 4.53 0.000 2.591236 6.536483
_cons | -4.670352 .5401517 -8.65 0.000 -5.72903 -3.611674
-------------+----------------------------------------------------------------
union |
black | .4397974 .0972261 4.52 0.000 .2492377 .6303572
tenure | .0997638 .0236575 4.22 0.000 .053396 .1461317
south | -.4895032 .0933276 -5.24 0.000 -.6724221 -.3065844
_cons | -.9679795 .0746464 -12.97 0.000 -1.114284 -.8216753
-------------+----------------------------------------------------------------
hazard |
lambda | -2.093313 .5801968 -3.61 0.000 -3.230478 -.9561486
-------------+----------------------------------------------------------------
rho | -0.89172
sigma | 2.3475104
------------------------------------------------------------------------------
*/
第五个是手工完成两步估计法，估计结果与直接使用etregress两步估计法基本保持一致。值得注意的是：
-
union取值为1的样本和union取值为0的样本的IMR计算公式不同，因此如果混淆了公式或使用一个公式进行计算，最后的结果必然存在偏误，因为union取值为1的样本和union取值为0的样本均参与了第二步回归。 - 手工两步法下参与回归的样本数目均为1,210，这点可以从两张回归结果表中直观看出。之所以相比于全样本缺少483个样本，原因在于有几个变量存在缺失值，Stata在回归时自动将这些变量存在缺失值的样本剔除，这也与前面四个模型的回归结果保持一致。
*-（五）手工两步估计法
*- 第一阶段方程（选择方程）回归
probit union $xlist1 south, nolog
est store First
/*
Probit regression Number of obs = 1,210
LR chi2(3) = 56.54
Prob > chi2 = 0.0000
Log likelihood = -592.15536 Pseudo R2 = 0.0456
------------------------------------------------------------------------------
union | Coef. Std. Err. z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
black | .4397974 .0972261 4.52 0.000 .2492377 .6303572
tenure | .0997638 .0236575 4.22 0.000 .053396 .1461317
south | -.4895032 .0933276 -5.24 0.000 -.6724221 -.3065844
_cons | -.9679795 .0746464 -12.97 0.000 -1.114284 -.8216753
------------------------------------------------------------------------------
*/
*- 计算IMR
predict y_hat, xb
gen imr = normalden(y_hat) / normal( y_hat) if union != .
replace imr = -normalden(y_hat) / normal(-y_hat) if union == 0
*- 第二阶段方程回归
reg wage $xlist3 imr
est store Second
/*
Source | SS df MS Number of obs = 1,210
-------------+---------------------------------- F(7, 1202) = 92.54
Model | 2227.31385 7 318.187693 Prob > F = 0.0000
Residual | 4132.88562 1,202 3.43834078 R-squared = 0.3502
-------------+---------------------------------- Adj R-squared = 0.3464
Total | 6360.19947 1,209 5.26071089 Root MSE = 1.8543
------------------------------------------------------------------------------
wage | Coef. Std. Err. t P>|t| [95% Conf. Interval]
-------------+----------------------------------------------------------------
union | 4.56386 .836918 5.45 0.000 2.921877 6.205842
age | .1543231 .0194468 7.94 0.000 .1161696 .1924765
grade | .4225025 .029448 14.35 0.000 .3647272 .4802778
smsa | .8628628 .12708 6.79 0.000 .6135395 1.112186
black | -.9206944 .1427409 -6.45 0.000 -1.200743 -.6406455
tenure | .1003226 .0424118 2.37 0.018 .0171133 .1835319
imr | -2.093314 .4858841 -4.31 0.000 -3.046589 -1.140038
_cons | -4.670352 .5337275 -8.75 0.000 -5.717493 -3.623211
------------------------------------------------------------------------------
*/
最后是各个回归模型下IMR的对比，以及回归结果的汇总输出。
*- 对比imr、imr1和imr2
order imr1 imr2 imr union
*- 回归结果输出
local mlist "OLS MLE MLE_r TwoStep First Second"
esttab `mlist' using 处理效应模型.rtf, replace ///
b(%6.4f) t(%6.4f) ///
scalar(N r2_a r2_p) compress nogap ///
star(* 0.1 ** 0.05 *** 0.01) ///
mtitle(`mlist') ///
title("Treatment Effects Model")
/*
Treatment Effects Model
----------------------------------------------------------------------------------------
(1) (2) (3) (4) (5) (6)
OLS MLE MLE_r TwoStep First Second
----------------------------------------------------------------------------------------
main
union 1.0039*** 4.5639*** 4.5639***
(7.5026) (4.5346) (5.4532)
age 0.1476*** 0.1487*** 0.1487*** 0.1543*** 0.1543***
(7.5574) (7.6952) (7.1887) (7.9179) (7.9357)
grade 0.4369*** 0.4206*** 0.4206*** 0.4225*** 0.4225***
(14.8228) (14.3256) (11.1284) (14.5620) (14.3474)
smsa 0.9755*** 0.9117*** 0.9117*** 0.8629*** 0.8629***
(7.7872) (7.2992) (7.6015) (6.7102) (6.7899)
black -0.6183*** -0.7882*** -0.7882*** -0.9207*** 0.4398*** -0.9207***
(-4.9388) (-5.7659) (-5.5944) (-5.1881) (4.5234) (-6.4501)
tenure 0.2118*** 0.1524*** 0.1524*** 0.1003* 0.0998*** 0.1003**
(6.2574) (4.1235) (3.2167) (1.9338) (4.2170) (2.3654)
1.union 2.9458*** 2.9458***
(10.7135) (5.2195)
south -0.4895***
(-5.2450)
imr -2.0933***
(-4.3083)
_cons -4.3260*** -4.3516*** -4.3516*** -4.6704*** -0.9680*** -4.6704***
(-8.1386) (-8.2354) (-6.8302) (-8.6464) (-12.9675) (-8.7504)
----------------------------------------------------------------------------------------
union
black 0.4557*** 0.4557*** 0.4398***
(4.7571) (4.8749) (4.5234)
tenure 0.0872*** 0.0872*** 0.0998***
(3.7488) (3.4648) (4.2170)
south -0.5807*** -0.5807*** -0.4895***
(-6.8233) (-6.9341) (-5.2450)
_cons -0.8856*** -0.8856*** -0.9680***
(-12.2232) (-11.9748) (-12.9675)
----------------------------------------------------------------------------------------
/
athrho -0.6544*** -0.6544***
(-7.1891) (-2.8179)
lnsigma 0.7027*** 0.7027***
(23.9517) (9.2683)
----------------------------------------------------------------------------------------
hazard
lambda -2.0933***
(-3.6079)
----------------------------------------------------------------------------------------
N 1210 1210 1210 1210 1210 1210
r2_a 0.3369 0.3464
r2_p 0.0456
----------------------------------------------------------------------------------------
t statistics in parentheses
* p<0.1, ** p<0.05, *** p<0.01
*/
