本 Notebook 配合第三章「谁进入样本，谁被我们观察到」使用。代码使用 Stata 自带的 mroz.dta 数据，并通过 nbstata 在 Jupyter 中执行。
本 Notebook 的目标不是做完整论文复现，而是演示四件事：
结果变量 wage 只在劳动参与者中可见；
可以用选择方程解释 inlf；
Heckman 模型同时报告结果方程和选择方程；
rho = 0 的检验用于判断选择方程与结果方程的误差项是否独立。
运行环境
需要先在本地配置好 nbstata。进入 Notebook 后，选择 nbstata kernel。数据无需下载到本地，Stata 会通过 webuse 在线读入。
若无法联网，可在 Stata 中先执行：
webuse mroz.dta, clear save data /mroz.dta, replace 再把后续 webuse 改成 use data/mroz.dta, clear。
读入数据并查看变量
webuse "mroz.dta" , clear des wage educ exper expersq city age faminc kidslt6 inlf
Running D:\stata19/profile.do ...
Variable Storage Display Value
name type format label Variable label
-------------------------------------------------------------------------------
wage float %9.0g estimated wage from earns., hours
educ byte %9.0g years of schooling
exper byte %9.0g actual labor mkt exper
expersq int %9.0g exper^2
city byte %9.0g =1 if live in SMSA
age byte %9.0g woman's age in yrs
faminc float %9.0g family income, 1975
kidslt6 byte %9.0g # kids < 6 years
inlf byte %9.0g =1 if in labor force, 1975
看看工资变量为何存在样本选择
wage 只在 inlf = 1 的女性中被观察到。先检查样本量和工资缺失结构。
tab inlfcount if !missing (wage)count if missing (wage)
=1 if in |
labor |
force, 1975 | Freq. Percent Cum.
------------+-----------------------------------
0 | 325 43.16 43.16
1 | 428 56.84 100.00
------------+-----------------------------------
Total | 753 100.00
428
325
"== Summary statistics: inlf==1 (in the labor force)" if inlf == 1"== Summary statistics: inlf==0 (not in the labor force)" if inlf == 0
== Summary statistics: inlf==1 (in the labor force)
Variable | Obs Mean Std. dev. Min Max
-------------+---------------------------------------------------------
wage | 428 4.177682 3.310282 .1282 25
educ | 428 12.65888 2.285376 5 17
exper | 428 13.03738 8.055923 0 38
expersq | 428 234.7196 270.0434 0 1444
city | 428 .6401869 .4805071 0 1
-------------+---------------------------------------------------------
age | 428 41.97196 7.721084 30 60
faminc | 428 24130.42 11671.26 2400 91044
kidslt6 | 428 .1401869 .3919231 0 2
== Summary statistics: inlf==0 (not in the labor force)
Variable | Obs Mean Std. dev. Min Max
-------------+---------------------------------------------------------
educ | 325 11.79692 2.181995 5 17
exper | 325 7.461538 6.918567 0 45
expersq | 325 103.3938 196.7098 0 2025
city | 325 .6461538 .4788995 0 1
age | 325 43.28308 8.467796 30 60
-------------+---------------------------------------------------------
faminc | 325 21698.05 12728.15 1500 96000
kidslt6 | 325 .3661538 .6368995 0 3
直接在有工资样本中做 OLS
这一步只使用工资可见的女性样本。它可以作为 Heckman 结果的对照，但不能自动解释为总体中的教育回报。
reg wage educ exper expersq city if inlf == 1
Source | SS df MS Number of obs = 428
-------------+---------------------------------- F(4, 423) = 15.08
Model | 583.846708 4 145.961677 Prob > F = 0.0000
Residual | 4095.20623 423 9.6813386 R-squared = 0.1248
-------------+---------------------------------- Adj R-squared = 0.1165
Total | 4679.05293 427 10.9579694 Root MSE = 3.1115
------------------------------------------------------------------------------
wage | Coefficient Std. err. t P>|t| [95% conf. interval]
-------------+----------------------------------------------------------------
educ | 0.481 0.067 7.19 0.000 0.350 0.612
exper | 0.032 0.062 0.53 0.598 -0.089 0.154
expersq | -0.000 0.002 -0.14 0.887 -0.004 0.003
city | 0.449 0.318 1.41 0.158 -0.175 1.074
_cons | -2.561 0.929 -2.76 0.006 -4.387 -0.735
------------------------------------------------------------------------------
Heckman 极大似然估计
下面使用第三章正文中的设定。结果方程解释 wage，选择方程解释 inlf。选择方程中包含年龄、年龄平方、家庭收入、6 岁以下子女数量和教育。
heckman wage educ exper expersq city, /// /// nolog
Heckman selection model Number of obs = 753
(regression model with sample selection) Selected = 428
Nonselected = 325
Wald chi2(4) = 88.79
Log likelihood = -1474.133 Prob > chi2 = 0.0000
------------------------------------------------------------------------------
| Coefficient Std. err. z P>|z| [95% conf. interval]
-------------+----------------------------------------------------------------
wage |
educ | 0.671 0.076 8.86 0.000 0.523 0.820
exper | 0.072 0.033 2.17 0.030 0.007 0.137
expersq | -0.001 0.001 -1.02 0.306 -0.003 0.001
city | -0.079 0.195 -0.40 0.686 -0.461 0.304
_cons | -7.517 0.993 -7.57 0.000 -9.464 -5.570
-------------+----------------------------------------------------------------
inlf |
age | -0.032 0.033 -0.98 0.328 -0.097 0.032
|
c.age#c.age | 0.000 0.000 0.89 0.374 -0.000 0.001
|
faminc | -0.000 0.000 -3.19 0.001 -0.000 -0.000
kidslt6 | -0.162 0.049 -3.30 0.001 -0.259 -0.066
educ | 0.161 0.020 8.11 0.000 0.122 0.200
_cons | -0.941 0.749 -1.26 0.209 -2.409 0.528
-------------+----------------------------------------------------------------
/athrho | 2.766 0.229 12.10 0.000 2.318 3.214
/lnsigma | 1.426 0.040 35.74 0.000 1.347 1.504
-------------+----------------------------------------------------------------
rho | 0.992 0.004 0.981 0.997
sigma | 4.160 0.166 3.848 4.499
lambda | 4.128 0.171 3.792 4.464
------------------------------------------------------------------------------
LR test of indep. eqns. (rho = 0): chi2(1) = 164.38 Prob > chi2 = 0.0000
Heckman 两步法估计
两步法更接近教材中先估计选择方程、再构造 IMR 的叙述方式。实证论文中可以把两步法和极大似然法并列作为稳健性比较。
heckman wage educ exper expersq city, /// ///
Heckman selection model -- two-step estimates Number of obs = 753
(regression model with sample selection) Selected = 428
Nonselected = 325
Wald chi2(4) = 46.82
Prob > chi2 = 0.0000
------------------------------------------------------------------------------
| Coefficient Std. err. z P>|z| [95% conf. interval]
-------------+----------------------------------------------------------------
wage |
educ | 0.489 0.080 6.15 0.000 0.333 0.645
exper | 0.033 0.061 0.54 0.586 -0.087 0.154
expersq | -0.000 0.002 -0.17 0.867 -0.004 0.003
city | 0.447 0.316 1.41 0.158 -0.173 1.067
_cons | -2.741 1.368 -2.00 0.045 -5.423 -0.060
-------------+----------------------------------------------------------------
inlf |
age | -0.012 0.069 -0.17 0.864 -0.147 0.123
|
c.age#c.age | -0.000 0.001 -0.33 0.742 -0.002 0.001
|
faminc | 0.000 0.000 0.75 0.452 -0.000 0.000
kidslt6 | -0.853 0.114 -7.46 0.000 -1.077 -0.629
educ | 0.117 0.024 4.93 0.000 0.070 0.163
_cons | -0.140 1.514 -0.09 0.926 -3.108 2.828
-------------+----------------------------------------------------------------
/mills |
lambda | 0.132 0.740 0.18 0.858 -1.318 1.582
-------------+----------------------------------------------------------------
rho | 0.04272
sigma | 3.0946906
------------------------------------------------------------------------------
手动生成 IMR
手动计算 IMR 的核心是先估计选择方程，再利用选择方程的线性预测值构造逆米尔斯比率。以 Heckman 两步法为例，第一阶段通过 Probit 模型估计：
\[
P(inlf_i=1|Z_i)=\Phi(Z_i\gamma)
\]
其中，\(Z_i\gamma\) 是 Probit 模型的线性预测值。根据 Heckman (1979) 的推导，对于进入第二阶段回归的样本 (\(inlf_i=1\) )，需要加入逆米尔斯比率：
\[
IMR_i=\lambda_i=\frac{\phi(Z_i\hat{\gamma})}{\Phi(Z_i\hat{\gamma})}
\]
其中，\(\phi(\cdot)\) 和 \(\Phi(\cdot)\) 分别表示标准正态分布的密度函数和累计分布函数。
因此，手动计算 IMR 的步骤为：
使用 probit 估计选择方程；
使用 predict, xb 获取线性预测值 \(\widehat{Z_i\gamma}\) ；
根据上述公式计算 IMR；
将 IMR 加入结果方程进行第二阶段 OLS 回归。
对应的 Stata 代码如下：
webuse "mroz.dta" , clear "== 第一步：估计选择方程 Probit ==" probit inlf age c.age#c.age faminc kidslt6 educpredict xb , xb gen imr = normalden (xb ) / normal (xb )"== 查看 IMR ==" if inlf == 1"== 第三步：结果方程，加入 IMR ==" reg wage educ exper expersq city imr if inlf == 1
== 第一步：估计选择方程 Probit ==
Iteration 0: Log likelihood = -514.8732
Iteration 1: Log likelihood = -465.74979
Iteration 2: Log likelihood = -465.70761
Iteration 3: Log likelihood = -465.70761
Probit regression Number of obs = 753
LR chi2(5) = 98.33
Prob > chi2 = 0.0000
Log likelihood = -465.70761 Pseudo R2 = 0.0955
------------------------------------------------------------------------------
inlf | Coefficient Std. err. z P>|z| [95% conf. interval]
-------------+----------------------------------------------------------------
age | -0.012 0.069 -0.17 0.864 -0.147 0.123
|
c.age#c.age | -0.000 0.001 -0.33 0.742 -0.002 0.001
|
faminc | 0.000 0.000 0.75 0.452 -0.000 0.000
kidslt6 | -0.853 0.114 -7.46 0.000 -1.077 -0.629
educ | 0.117 0.024 4.93 0.000 0.070 0.163
_cons | -0.140 1.514 -0.09 0.926 -3.108 2.828
------------------------------------------------------------------------------
== 查看 IMR ==
Variable | Obs Mean Std. dev. Min Max
-------------+---------------------------------------------------------
imr | 428 .6207169 .2484473 .1687163 1.792734
== 第三步：结果方程，加入 IMR ==
Source | SS df MS Number of obs = 428
-------------+---------------------------------- F(5, 422) = 12.04
Model | 584.152429 5 116.830486 Prob > F = 0.0000
Residual | 4094.90051 422 9.7035557 R-squared = 0.1248
-------------+---------------------------------- Adj R-squared = 0.1145
Total | 4679.05293 427 10.9579694 Root MSE = 3.1151
------------------------------------------------------------------------------
wage | Coefficient Std. err. t P>|t| [95% conf. interval]
-------------+----------------------------------------------------------------
educ | 0.489 0.080 6.11 0.000 0.331 0.646
exper | 0.033 0.062 0.54 0.589 -0.088 0.155
expersq | -0.000 0.002 -0.17 0.868 -0.004 0.003
city | 0.447 0.318 1.40 0.161 -0.179 1.073
imr | 0.132 0.745 0.18 0.859 -1.332 1.596
_cons | -2.741 1.377 -1.99 0.047 -5.449 -0.034
------------------------------------------------------------------------------
该代码得到的 imr 变量以及第二阶段回归结果，应与使用 heckman 并附加 twostep 选项的两步法估计结果一致：
heckman wage educ exper expersq city, /// ///
Heckman selection model -- two-step estimates Number of obs = 753
(regression model with sample selection) Selected = 428
Nonselected = 325
Wald chi2(4) = 46.82
Prob > chi2 = 0.0000
------------------------------------------------------------------------------
| Coefficient Std. err. z P>|z| [95% conf. interval]
-------------+----------------------------------------------------------------
wage |
educ | 0.489 0.080 6.15 0.000 0.333 0.645
exper | 0.033 0.061 0.54 0.586 -0.087 0.154
expersq | -0.000 0.002 -0.17 0.867 -0.004 0.003
city | 0.447 0.316 1.41 0.158 -0.173 1.067
_cons | -2.741 1.368 -2.00 0.045 -5.423 -0.060
-------------+----------------------------------------------------------------
inlf |
age | -0.012 0.069 -0.17 0.864 -0.147 0.123
|
c.age#c.age | -0.000 0.001 -0.33 0.742 -0.002 0.001
|
faminc | 0.000 0.000 0.75 0.452 -0.000 0.000
kidslt6 | -0.853 0.114 -7.46 0.000 -1.077 -0.629
educ | 0.117 0.024 4.93 0.000 0.070 0.163
_cons | -0.140 1.514 -0.09 0.926 -3.108 2.828
-------------+----------------------------------------------------------------
/mills |
lambda | 0.132 0.740 0.18 0.858 -1.318 1.582
-------------+----------------------------------------------------------------
rho | 0.04272
sigma | 3.0946906
------------------------------------------------------------------------------
报告时至少说明什么
在论文或讲义中报告 Heckman 模型时，至少要说明：
为什么结果变量只在选择样本中可见；
选择方程中的变量为什么能解释进入样本；
哪些变量只进入选择方程，为什么可以作为排他性变量；
rho = 0 或 IMR 系数是否显著；OLS 和 Heckman 结果是否有实质差异；
Heckman 的使用是否对应样本选择问题，而不是把它当成泛化的内生性修正工具。