三、樣本選擇模型的Stata實(shí)操
Stata關(guān)于樣本選擇模型的官方命令是heckman，該命令能夠進(jìn)行Heckman兩步估計(jì)以及模型整體參數(shù)的極大似然估計(jì)。
3.1 語(yǔ)法介紹
關(guān)于heckman更多語(yǔ)法細(xì)節(jié)，在命令窗口鍵入如下代碼即可了解。
help heckman
常用的語(yǔ)法如下。
**# 一、樣本選擇模型
**# 1.1 語(yǔ)法介紹
/*
a. 兩步估計(jì)法
heckman y x1 x2, select( x1 x2 z1) twostep mills(newname)
等價(jià)于
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ǔ)法書(shū)寫(xiě)大致包括兩部分，選擇項(xiàng)之前的是第二階段回歸的被解釋變量y以及控制變量x1和x2，不需要寫(xiě)入IMR。選擇項(xiàng)則包括對(duì)選擇方程的具體設(shè)置，以及其他的細(xì)節(jié)。
-
select()表示寫(xiě)入選擇方程，括號(hào)內(nèi)的是選擇方程的具體變量，有兩種寫(xiě)法。- 第一種寫(xiě)法：直接寫(xiě)入控制變量（
x1和x2）和外生變量z1。此時(shí)嚴(yán)格要求原回歸的被解釋變量y存在缺失值，且缺失值不能以0作為標(biāo)記。 - 第二種寫(xiě)法：首先寫(xiě)入
y是否被觀測(cè)到的虛擬變量y_dummy（y存在缺失值的樣本y_dummy標(biāo)記為0，不存在缺失值的樣本標(biāo)記為1，且y取值為0不作為缺失值處理）作為選擇方程的被解釋變量，等號(hào)后面跟著控制變量和外生變量。若選擇這種寫(xiě)法，則要求提前生成y_dummy。
- 第一種寫(xiě)法：直接寫(xiě)入控制變量（
-
twostep表示使用兩步估計(jì)法，默認(rèn)使用MLE。 -
mills()表示生成各個(gè)樣本的IMR，并以newname作為變量名。nshazard()的作用相同。 -
nolog表示在使用MLE時(shí)回歸結(jié)果不顯示迭代過(guò)程。 -
vce()表示使用穩(wěn)健標(biāo)準(zhǔn)誤，括號(hào)內(nèi)填入robust表示使用異方差穩(wěn)健標(biāo)準(zhǔn)誤；填入cluster varname表示使用聚類(lèi)穩(wěn)健標(biāo)準(zhǔn)誤，以變量varname作為聚類(lèi)標(biāo)準(zhǔn)，根據(jù)經(jīng)驗(yàn)法則，要求聚類(lèi)數(shù)目大于30。需要注意的是，vce()不能在兩步法中使用。
3.2 實(shí)例操作
heckman的help文件附帶演示案例，下面根據(jù)演示案例中的數(shù)據(jù)和代碼，對(duì)樣本選擇模型在Stata中的具體實(shí)現(xiàn)進(jìn)行介紹。
首先調(diào)用數(shù)據(jù)集womenwk，由于該數(shù)據(jù)集并非Stata自帶，因此可能無(wú)法調(diào)用成功，這種情況下可以直接調(diào)用本次推文提供的數(shù)據(jù)。
*- Stata Version: 16 | 17
*- 定義路徑
cd "C:\Users\KEMOSABE\Desktop\heckman"
**# 1.2 實(shí)例操作
webuse womenwk, clear
/*
*- 或者直接調(diào)用已下載的數(shù)據(jù)集
use womenwk.dta, clear
*/
該例中，基準(zhǔn)回歸的被解釋變量是wage，解釋變量是educ和age；選擇方程中額外引入了兩個(gè)外生解釋變量married和children。
由于被解釋變量wage存在缺失值，因此我們懷疑基準(zhǔn)回歸方程（OLS）可能存在樣本選擇偏差，因此選用樣本選擇模型做進(jìn)一步的穩(wěn)健性檢驗(yàn)。
該數(shù)據(jù)集為截面數(shù)據(jù)，因此標(biāo)準(zhǔn)誤無(wú)法聚類(lèi)到個(gè)體層面，而同一地區(qū)的個(gè)體特征可能存在一定的相關(guān)性，因此我們嘗試將標(biāo)準(zhǔn)誤聚類(lèi)到county層面，下面判斷是否能夠使用聚類(lèi)標(biāo)準(zhǔn)誤。
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
*/
可以看到，地區(qū)數(shù)目只有10個(gè)，遠(yuǎn)無(wú)法滿足聚類(lèi)數(shù)目最低值（30個(gè)）的要求，因此選用異方差穩(wěn)健標(biāo)準(zhǔn)誤。
根據(jù)wage是否存在缺失值生成虛擬變量wage_dummy?？梢钥吹剑?,000個(gè)樣本中，wage數(shù)據(jù)缺失的有657個(gè)。
*- 查看y的缺失值情況
nmissing wage
/*
wage 657
*/
*- 生成y是否取值的虛擬變量y_dummy
gen wage_dummy = (wage != .)
然后是定義控制變量和外生變量的全局暫元。
*- 定義全局暫元
global xlist educ age
global zlist married children
第一個(gè)是基準(zhǔn)回歸，可以看到：
- 參與回歸的樣本數(shù)目為1,343個(gè)，即
wage存在缺失值的樣本（657個(gè)）在回歸時(shí)直接被drop掉。 - 基準(zhǔn)回歸中兩個(gè)解釋變量的系數(shù)均顯著為正，模型擬合程度也較好（修正R方為0.2524）。
*-（一）基準(zhǔn)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
------------------------------------------------------------------------------
*/
第二個(gè)是樣本選擇模型，使用MLE方法進(jìn)行估計(jì)，可以看到：
- 選擇方程中兩個(gè)外生變量均顯著為正，說(shuō)明外生變量的選擇是有效的。
- 在第二階段回歸中，
IMR（即lambda）的估計(jì)系數(shù)為4.2244，但顯著性未知，該值等于rho和sigma的乘積，其中：-
sigma是原方程干擾項(xiàng)的標(biāo)準(zhǔn)差。 -
rho是選擇方程干擾項(xiàng)和第二階段回歸干擾項(xiàng)的相關(guān)系數(shù)。如果rho等于0，表示第二階段回歸中IMR的系數(shù)不顯著，說(shuō)明樣本選擇偏差在原方程中不怎么嚴(yán)重，反之則需要考慮樣本選擇偏差帶來(lái)的估計(jì)偏誤?；貧w結(jié)果的末尾是LR檢驗(yàn)，檢驗(yàn)的原假設(shè)是H0: rho = 0，p值說(shuō)明至少可以在1%的水平下拒絕原假設(shè)，可以認(rèn)為rho顯著不等于0，這說(shuō)明原模型中確實(shí)存在嚴(yán)重的樣本選擇偏差，基準(zhǔn)回歸結(jié)果不可信。
-
- 第二階段回歸結(jié)果中，兩個(gè)解釋變量仍舊顯著為正，且相較于基準(zhǔn)回歸結(jié)果取值變化不大，說(shuō)明考慮到樣本選擇偏差后基準(zhǔn)回歸結(jié)果依然是穩(wěn)健的。
- 回歸結(jié)果右上角說(shuō)明總的樣本量是2,000個(gè)，參與選擇方程回歸的樣本為2,000個(gè)，參與第二階段回歸的樣本為1,343個(gè)。
*-（二）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
*/
第三個(gè)同樣是樣本選擇模型，使用MLE方法進(jìn)行估計(jì)，這里將使用異方差穩(wěn)健標(biāo)準(zhǔn)誤。可以看到，使用穩(wěn)健標(biāo)準(zhǔn)誤的回歸結(jié)果與上文使用普通標(biāo)準(zhǔn)誤的回歸結(jié)果基本保持一致。需要注意的是，在使用穩(wěn)健標(biāo)準(zhǔn)誤后，對(duì)H0的統(tǒng)計(jì)檢驗(yàn)將變成Wald檢驗(yàn)，這里Wald檢驗(yàn)同樣拒絕原假設(shè)。
*-（三）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
*/
第四個(gè)是樣本選擇模型，使用兩步法進(jìn)行估計(jì)，可以看到：
- 選擇方程中兩個(gè)工具變量的引入是有效的。
- 第二階段回歸中，
IMR的回歸系數(shù)等于4.0016，與MLE方法下的4.2244相差不大，但兩步法下IMR回歸系數(shù)可以直接進(jìn)行z檢驗(yàn)，并且統(tǒng)計(jì)結(jié)果說(shuō)明IMR回歸系數(shù)至少在1%的水平下顯著為正，這同時(shí)說(shuō)明原方程中的樣本選擇偏差問(wèn)題不可忽視。 - 第二階段回歸結(jié)果中，兩個(gè)解釋變量仍舊顯著為正，且大小與基準(zhǔn)回歸結(jié)果相比變化不大，這說(shuō)明在考慮樣本選擇偏差的情況下，基準(zhǔn)回歸結(jié)果是可信的。
- 右上角同樣說(shuō)明在兩步估計(jì)法下，參與選擇方程的樣本數(shù)是2,000個(gè)，但參與第二階段回歸的樣本數(shù)只有1,343個(gè)，減少的657個(gè)樣本與被解釋變量
wage存在缺失值的樣本完全重合。
*-（四）兩步估計(jì)法
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
------------------------------------------------------------------------------
*/
第五個(gè)是手工完成兩步估計(jì)法，但實(shí)際使用時(shí)這種方法十分不推薦（包括上文使用兩步法估計(jì)的樣本選擇模型），因?yàn)閮刹焦烙?jì)法下第一步回歸的估計(jì)誤差會(huì)帶到第二步，造成效率損失，影響第二步估計(jì)系數(shù)的顯著性。
手工兩步估計(jì)法的思路是先用控制變量和外生變量對(duì)wage_dummy做probit回歸，計(jì)算出擬合值y_hat，然后根據(jù)y_hat計(jì)算出全部樣本的IMR，最后將IMR作為額外的控制變量引入基準(zhǔn)回歸方程中做進(jìn)一步的回歸。
第一階段（選擇方程）回歸的代碼與結(jié)果如下。
*-（五）手工兩步估計(jì)法
*- 第一階段方程（選擇方程）回歸
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
------------------------------------------------------------------------------
*/
可以看到，實(shí)際參與回歸的樣本為2,000個(gè)，控制變量和兩個(gè)外生變量均顯著為正，這說(shuō)明外生變量的選擇是有效的。此外，衡量模型整體擬合程度的統(tǒng)計(jì)指標(biāo)是偽R方（Pseudo R2），偽R方為0.1889，整體擬合程度較好。
下面根據(jù)擬合值y_hat計(jì)算出額外的控制變量imr。其中，normalden()是標(biāo)準(zhǔn)正態(tài)的概率密度函數(shù)pdf，normal()是標(biāo)準(zhǔn)正態(tài)的累積分布函數(shù)cdf。
*- 計(jì)算IMR
predict y_hat, xb
gen imr = normalden(y_hat) / normal(y_hat)
獲得各個(gè)樣本的imr后，就可以進(jìn)行第二階段回歸?？梢园l(fā)現(xiàn)：
- 與樣本選擇模型的兩步估計(jì)法結(jié)果相比，手工兩步法估計(jì)結(jié)果在系數(shù)值大小方面沒(méi)有任何改變，在系數(shù)標(biāo)準(zhǔn)誤方面變化也不大，從而各個(gè)變量的系數(shù)顯著性保持高度一致。
- 這提示我們，兩步估計(jì)法（包括手工兩步法）帶來(lái)的效率損失不可忽視，雖然本例中的結(jié)果和MLE以及基準(zhǔn)回歸基本保持一致，實(shí)際使用應(yīng)盡量避免使用這種方法，除非能夠找到比較好的方法來(lái)校正誤差。
- 第二階段實(shí)際參與回歸的樣本是1,343個(gè)。
- 衡量模型整體擬合程度的統(tǒng)計(jì)指標(biāo)是修正R方，修正R方等于0.2777，整體擬合程度較高。事實(shí)上，由于該例并未規(guī)定原方程的核心解釋變量，因此對(duì)模型整體擬合程度的考察是有意義的，但在現(xiàn)實(shí)應(yīng)用中，R方即便很小問(wèn)題也不大，只要求核心解釋變量的系數(shù)大小和顯著性符合預(yù)期即可。就算對(duì)R方有具體的要求，只要我們?cè)诜匠讨幸胱銐蚨嗟目刂谱兞?，控制足夠多的“固定效?yīng)”，R方總能達(dá)到預(yù)定的要求。因此，實(shí)際應(yīng)用中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
------------------------------------------------------------------------------
*/
可以對(duì)比一下用各個(gè)方法計(jì)算的IMR。其中，imr1是用樣本選擇模型的MLE方法計(jì)算的，imr2是用樣本選擇模型的兩步估計(jì)法計(jì)算的，imr是用手工兩步法計(jì)算的。打開(kāi)Stata的數(shù)據(jù)編輯器，可以發(fā)現(xiàn)：
- 所有樣本的
IMR均有取值。 -
imr和imr2保持高度一致（imr四舍五入了），因?yàn)閮烧叨际怯脙刹焦烙?jì)法計(jì)算出來(lái)的。 -
imr1和imr2存在差距，但差別不大。
*- 對(duì)比imr、imr1和imr2
order imr1 imr2 imr wage_dummy
最后是五個(gè)模型回歸結(jié)果的對(duì)比，可以發(fā)現(xiàn)：
- 兩個(gè)解釋變量（
education和age）的回歸系數(shù)的符號(hào)、大小、顯著性在模型中基本保持一致，并且與模型
的區(qū)別不大，這說(shuō)明在考慮到樣本選擇偏差后，根據(jù)基本回歸結(jié)果得出的結(jié)論依然成立。
- 模型
缺少變量
lambda，在實(shí)際匯報(bào)時(shí)應(yīng)該根據(jù)原始回歸表進(jìn)行匯報(bào)。 - 模型
匯報(bào)的樣本量是2,000個(gè)全樣本，但應(yīng)該清楚不同回歸階段實(shí)際參與的樣本數(shù)目是不同的。
- 關(guān)于模型擬合程度的統(tǒng)計(jì)量，模型
匯報(bào)的都是修正R方，模型
的是偽R方，模型
Wald chi2，因此，在實(shí)際匯報(bào)時(shí)，應(yīng)該把Wald chi2作為這三個(gè)模型整體擬合度的匯報(bào)指標(biāo)。
*- 回歸結(jié)果輸出
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
*/
以上命令、代碼和實(shí)例均基于截面數(shù)據(jù)，現(xiàn)階段的文獻(xiàn)大多數(shù)也是將面板數(shù)據(jù)轉(zhuǎn)為截面數(shù)據(jù)處理，最多控制幾個(gè)個(gè)體虛擬變量和時(shí)間虛擬變量，并將其視為面板數(shù)據(jù)回歸。
事實(shí)上，Stata在16及以上版本中提供了估計(jì)面板數(shù)據(jù)樣本選擇模型官方命令xtheckman，但這個(gè)命令只能估計(jì)隨機(jī)效應(yīng)模型，估計(jì)固定效應(yīng)模型需安裝一個(gè)外部命令xtheckmanfe，該命令由Boston College的Fernando編寫(xiě)，可以鍵入如下代碼進(jìn)行安裝或更新。
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]
四、處理效應(yīng)模型的Stata實(shí)操
處理效應(yīng)模型在Stata中的官方命令是treatreg和etregress，其中treatreg在Stata第14版之后被etregress替代，因此下面所使用的均為etregress。
與heckman類(lèi)似，etregress同樣可以進(jìn)行處理效應(yīng)模型的兩步法估計(jì)和極大似然估計(jì)。
4.1 語(yǔ)法介紹
etregress的語(yǔ)法結(jié)構(gòu)和heckman極為相似，但還是存在幾點(diǎn)不同，具體來(lái)說(shuō)：
- 第一，選擇方程的寫(xiě)法是
treat()，而非select()。 - 第二，選擇方程中的被解釋變量是原方程中的核心解釋變量
D，該解釋變量取值為0或1，且不存在缺失值。在etregress語(yǔ)法中，變量D不可省略。 - 第三，選擇方程中的工具變量
z1直接影響是變量D，而非樣本選擇模型中的y_dummy。 - 第四，處理效應(yīng)模型中生成
IMR的選擇項(xiàng)是hazard()，而非樣本選擇模型中的mills()或nshazad()。 - 第五，雖然第二步回歸方程中沒(méi)有寫(xiě)入
D，但該命令在進(jìn)行第二部回歸時(shí)默認(rèn)自動(dòng)引入D。因此，如果在語(yǔ)法書(shū)寫(xiě)時(shí)在第二步回歸方程中手動(dòng)加入D，正式回歸時(shí)該變量會(huì)產(chǎn)生完全的多重共線性而被omitted掉。 - 第六，對(duì)于工具變量
z1的選擇，還是那句話，盡量避免使用D的滯后項(xiàng)D_lag，如果實(shí)在找不到一個(gè)良好的工具變量而選擇使用D_lag，也應(yīng)該要確保面板結(jié)構(gòu)中不同個(gè)體D取值為1的時(shí)點(diǎn)是交錯(cuò)的，這樣才能避免產(chǎn)生多重共線性問(wèn)題。當(dāng)然，如果是純粹的截面數(shù)據(jù)（非混合面板），就不存在所謂的滯后變量，在這種情況下找到一個(gè)良好的工具變量是必須的。
**# 二、處理效應(yīng)模型
**# 2.1 語(yǔ)法介紹
/*
a. 兩步估計(jì)法
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 實(shí)例操作
這里使用的是etregress命令的help文件中的示例數(shù)據(jù)集union3與代碼。
同樣，如果非內(nèi)嵌數(shù)據(jù)集調(diào)用失敗可以直接調(diào)用本次推文提供的數(shù)據(jù)。
**# 2.2 實(shí)例操作
webuse union3, clear
/*
*- 或者直接調(diào)用已下載的數(shù)據(jù)集
use union3.dta, clear
*/
由于union3為截面數(shù)據(jù)集，因此標(biāo)準(zhǔn)誤無(wú)法聚類(lèi)到個(gè)體層面，同時(shí)考慮到同一行業(yè)內(nèi)不同個(gè)體的某些特征存在一致性，因此考慮將標(biāo)準(zhǔn)誤聚類(lèi)到行業(yè)層面。
*- 判斷是否能夠使用聚類(lèi)標(biāo)準(zhǔn)誤
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
*/
可見(jiàn)，如果將標(biāo)準(zhǔn)誤聚類(lèi)到行業(yè)層面，聚類(lèi)數(shù)目少于經(jīng)驗(yàn)數(shù)目（30個(gè)），聚類(lèi)數(shù)目過(guò)少將導(dǎo)致標(biāo)準(zhǔn)誤無(wú)法收斂至真實(shí)值（該問(wèn)題在往期推文『雙重差分法(DID) | 空間DID』中有討論），因此該例還是使用異方差穩(wěn)健標(biāo)準(zhǔn)誤。
下面觀察核心解釋變量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
*/
可見(jiàn)，總共有1,693個(gè)樣本，其中union取值為1的樣本有260個(gè)，取值為0的樣本有984個(gè)，缺失值樣本有449個(gè)。剔除掉缺失值樣本，對(duì)于剩下的1,244個(gè)樣本，我們懷疑存在自選擇偏差，因此使用處理效應(yīng)模型來(lái)檢驗(yàn)?zāi)Ｐ椭惺欠翊嬖谠搯?wèn)題，并與基準(zhǔn)回歸結(jié)果對(duì)比進(jìn)行穩(wěn)健性檢驗(yàn)。
下面將定義變量的全局暫元，并對(duì)所使用的變量進(jìn)行簡(jiǎn)要說(shuō)明。
*- 定義全局暫元
global xlist1 black tenure
global xlist2 age grade smsa black tenure
global xlist3 union age grade smsa black tenure
其中，wage是被解釋變量，union是核心解釋變量，并且我們其懷疑存在自選擇問(wèn)題。age、grade、smsa、black和tenure是基準(zhǔn)回歸方程的控制變量；black和tenure是選擇方程的控制變量，south是選擇方程的工具變量。
可以看到，選擇方程的控制變量并未與基準(zhǔn)回歸的控制變量重合，其數(shù)目反而少于基準(zhǔn)回歸控制變量的數(shù)目。
第一個(gè)是基準(zhǔn)回歸，可以看到：
- 所有解釋變量的系數(shù)均至少在1%的水平下顯著，特別是我們重點(diǎn)關(guān)注的核心解釋變量
union。 - 由于樣本中各變量存在缺失值，因此最終參與回歸的樣本只有1,210個(gè)。
- 由于各樣本
union的取值可能不是隨機(jī)的，即可能存在自選擇問(wèn)題，導(dǎo)致基準(zhǔn)回歸結(jié)果存在偏誤，因此還需要做進(jìn)一步的穩(wěn)健性檢驗(yàn)才能得出令人信服的結(jié)論。
*-（一）基準(zhǔn)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
------------------------------------------------------------------------------
*/
第二個(gè)是處理效應(yīng)模型，使用MLE方法從模型整體角度對(duì)參數(shù)進(jìn)行估計(jì)?？梢钥吹剑?/p>
- 第一步回歸中
south顯著為負(fù)，說(shuō)明工具變量的選擇有效。 - 第二步回歸中核心解釋變量
union顯著為正，這一點(diǎn)在處理效應(yīng)模型回歸結(jié)果表中顯示為1.union，說(shuō)明相對(duì)于union取值為0的樣本，union取值為1的樣本的平均工資高2.9458個(gè)單位。2.9458大約是基準(zhǔn)回歸結(jié)果1.0039的三倍，這說(shuō)明在未考慮自選擇偏差的情況下，低估了union對(duì)wage的影響。 - 其余解釋變量估計(jì)系數(shù)的符號(hào)和顯著性基本沒(méi)變，但估計(jì)系數(shù)的數(shù)值改變較大，這說(shuō)明自選擇偏差帶來(lái)的估計(jì)偏誤不能忽視。
-
IMR（lambda）的估計(jì)系數(shù)為-1.1603，同時(shí)rho的LR檢驗(yàn)結(jié)果說(shuō)明至少在1%的水平下拒絕H0。這都說(shuō)明基準(zhǔn)回歸模型中確實(shí)存在嚴(yán)重的自選擇偏差。 - 第一步回歸和第二步回歸的樣本數(shù)均為1,210個(gè)。
*-（二）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
*/
第三個(gè)是處理效應(yīng)模型，使用MLE方法并使用異方差穩(wěn)健標(biāo)準(zhǔn)誤?？梢钥吹剑c使用普通標(biāo)準(zhǔn)誤的MLE方法相比，除系數(shù)的z統(tǒng)計(jì)量有細(xì)微差別，以及對(duì)rho的顯著性檢驗(yàn)使用的是Wald統(tǒng)計(jì)指標(biāo)，其他結(jié)果沒(méi)有實(shí)質(zhì)性差異。
*-（三）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
*/
第四個(gè)是用兩步法估計(jì)處理效應(yīng)模型，可以發(fā)現(xiàn)：
- 選擇方程中的工具變量選擇有效。
- 第二階段回歸中的所有解釋變量，除
tenure，均至少在1%的水平下顯著，lambda顯著為負(fù)。 - 回歸系數(shù)的大小與MLE方法下的估計(jì)系數(shù)相近，但與基準(zhǔn)回歸相差較大，同時(shí)考慮到
IMR顯著，這說(shuō)明原模型中的自選擇偏差不可忽視，在未考慮自選擇偏差的情況下，將嚴(yán)重低估union對(duì)wage的影響。
*-（四）兩步估計(jì)法
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
------------------------------------------------------------------------------
*/
第五個(gè)是手工完成兩步估計(jì)法，估計(jì)結(jié)果與直接使用etregress兩步估計(jì)法基本保持一致。值得注意的是：
-
union取值為1的樣本和union取值為0的樣本的IMR計(jì)算公式不同，因此如果混淆了公式或使用一個(gè)公式進(jìn)行計(jì)算，最后的結(jié)果必然存在偏誤，因?yàn)?code>union取值為1的樣本和union取值為0的樣本均參與了第二步回歸。 - 手工兩步法下參與回歸的樣本數(shù)目均為1,210，這點(diǎn)可以從兩張回歸結(jié)果表中直觀看出。之所以相比于全樣本缺少483個(gè)樣本，原因在于有幾個(gè)變量存在缺失值，Stata在回歸時(shí)自動(dòng)將這些變量存在缺失值的樣本剔除，這也與前面四個(gè)模型的回歸結(jié)果保持一致。
*-（五）手工兩步估計(jì)法
*- 第一階段方程（選擇方程）回歸
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
------------------------------------------------------------------------------
*/
*- 計(jì)算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
------------------------------------------------------------------------------
*/
最后是各個(gè)回歸模型下IMR的對(duì)比，以及回歸結(jié)果的匯總輸出。
*- 對(duì)比imr、imr1和imr2
order imr1 imr2 imr union
*- 回歸結(jié)果輸出
local mlist "OLS MLE MLE_r TwoStep First Second"
esttab `mlist' using 處理效應(yīng)模型.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
*/
