Heckman两步法 | 样本选择模型 & 处理效应模型

大家好，欢迎来到IT知识分享网。

这期推送简单介绍一下样本选择模型和处理效应模型，其中样本选择模型是一般意义上的Heckman两步法，后者则借鉴了Heckman两步法的构建思想，但又不完全等同于前者。模型介绍之后，将利用help文件中的示例数据与代码简单演示一下这两个模型在Stata中的具体操作，然后简单评述一下现阶段文献中对这两个模型的理解与应用情况，最后结合一篇论文的公开数据与代码进行结果复盘与二次验证。

$Note:$ 1、下划线字体为链接，可点击跳转；2、推文中的公式与代码块均可左右滑动；3、该文首发于微信公众号DMETP，欢迎关注；4、需要本次推送所使用的数据和代码的朋友，可以在公众号后台对话框内回复关键词heckman。

一、样本选择偏差与自选择偏差

上期推送『双重差分法 | PSM – DID』介绍了样本选择偏差与自选择偏差的区别，最关键的一点在于两者非随机的选择机制是不同的。

• 样本选择偏差。样本选择偏差的非随机选择机制在于对样本的选择不随机。在样本数据的采集过程中，只对某部分群体进行调查，但这部分群体与其他群体在某些方面的特征差异较大，因此根据这样的样本做回归得到的普适性结论并不可信。体现在具体的数据集中就是，数据集中只有特定群体的样本，或者，虽然有全部群体的所有解释变量数据，但除特定群体之外的其他群体的被解释变量数据缺失，在这两种情况下进行的回归，都将直接忽视其他群体的样本信息（y缺失的样本在参与回归时将被drop掉）。实质上，样本选择偏差说的就是参与回归的样本不能代表总体从而产生估计偏误的问题。
• 自选择偏差。自选择偏差的非随机选择机制在于对自变量的选择不随机。在使用DID方法评估政策效应时，一个明显的事实就是，相对于未实施政策的地区（控制组），实施政策的地区（处理组）通常情况下经济发展都较为发达、各类基础设施建设都较为完善，而所谓的“政策效果评估”也即考察政策的经济效应，因此地区是否参与政策这一行为是内生的。体现在回归方程中就是，经济指标（如，GDP、人均GDP、GDP增长率等）作为被解释变量y，地区（在某时点）是否实施该项政策的哑元变量D作为核心解释变量，但由于政策内生，因此某些影响地区是否参与决策D的（可观测或不可观测）因素也将同时影响经济指标y，由于这些因素或者无法穷尽、或者影响形式未知、或者不可测度，因此被放到随机扰动项中，造成解释变量D与扰动项ε相关，即$Cov(D,\epsilon )\ne 0$。实质上，自选择偏差说的就是实验组与控制组的先验条件存在较大差异从而导致估计偏误的问题。
• 两者的区别。非随机选择机制的不同是两者最大的区别，体现在具体回归方程中就是，样本选择偏差中被解释变量y是否被观测到或是否取值（而非取值大小）是非随机的；而自选择偏差中哑元解释变量D的取值是非随机的。陈强（2014）《高级计量经济学及Stata应用（第二版）》第539页认为，样本选择问题通常不考虑某项目或政策的效应，故个体间的差异并不在于是否得到处理，而在于是否能进入样本（即被解释变量$y_i$是否可观测），通常$D_i=1$意味着$y_i$可观测，而$D_i=0$则意味着$y_i$不可观测。而在处理效应模型中，无论$D_i=1$或$0$，结果变量$y_i$均可观测。这种说法基本概括了两者的区别，但有一个小问题，在样本选择偏差中，$D_i$的取值与$y_i$是否可观测并不存在必然的关系，因为$D_i$是一个确定并可准确测度的因素，而影响$y_i$是否可观测的却是一个不可观测的潜变量，这个潜变量由一系列控制变量与外生变量决定。

二、两个模型的估计思路

花大篇幅论述样本选择偏差与自选择偏差这两个问题，自然是为了引出解决这两个问题的具体方法。

2.1 样本选择模型

对于样本选择偏差导致的估计偏误，将使用样本选择模型（Sample Selection Model）来缓解。样本选择偏差与样本选择模型（或称Heckman两步估计法、Heckit）由诺贝尔经济学奖获得者Heckman教授于1979年提出。

[2] Heckman J J. Sample Selection Bias as a Specification Error[J]. Econometrica, 1979, 47(01): 153-161.

本质上，样本选择偏差其实是一个因遗漏变量而导致内生性的特例（具体推导请看任意一本高级计量经济学教材，如陈强（2014）《高级计量经济学及Stata应用（第二版）》第234页、Hansen（2021）《ECONOMETRICS（Version 2021）》第852页等）。

回归方程中被遗漏的变量叫做逆米尔斯比率（Inverse Mill’s Ratio，IMR或$\lambda$），也被称为风险函数（Hazard Function），计算公式为：

$$IM{R}_i(或{\lambda }_i)=\frac{\varphi ({\hat{y}}_i)}{\Phi ({\hat{y}}_i)}\text{\hspace{1em}(1)}$$

其中，${\hat{y}}_i$为第$i$个样本在第一步回归（选择方程）的拟合值，$\varphi (\cdot )$为标准正态的概率密度函数（Probability Density Function，pdf），$\Phi (\cdot )$为累积分布函数（Cumulative Distribution Function，cdf）。

因此，样本选择模型的估计思路是：首先，计算全部样本的IMR；随后，将遗漏变量IMR代入原回归方程中，具体来说：

• 第一步 ：用probit方法估计选择方程，其中原回归方程的被解释变量y是否被观测到或是否取值的虚拟变量y_dummy作为probit的被解释变量，解释变量包括原回归方程所有解释变量和至少一个外生变量，该外生变量只影响y是否取值，而不影响y的大小，即满足相关性和外生性的要求（但不是工具变量）。估计出所有变量的系数后，将样本数据代入至probit模型中，计算出拟合值$\hat{y}$，再将$\hat{y}$代入风险函数（公式$(1)$）中计算出IMR。这里有四点需要注意：
  • 第一，选择方程的被解释变量是原回归方程中被解释变量y是否被观测到或是否取值的虚拟变量，即y_dummy，当y取值不为空（包括取值为0）时，y_dummy等于1，只有当y_dummy取值为空（missing）时，y_dummy才等于0。关于这一点，现实应用中存在的问题是，即便我们十分清楚存在样本选择偏差，但由于前期数据搜集过程中直接忽视了y取值为空的样本，因此无法采用样本选择模型，因为样本选择模型第一步选择方程使用的是所有样本，包括y取值为空的样本和取值不为空的样本。由于数据搜集过程存在问题，因此许多文献使用的所谓Heckman两步法实际上是一种“伪样本选择模型”，与Heckman（1979）提出的两步估计法（Two-Step Estimation，或Heckit）完全不同，而且也不是下文将要介绍的处理效应模型。
  • 第二，选择方程的被解释变量只能是原回归方程中被解释变量y是否被观测到或是否取值的虚拟变量，而不能是其他变量，更不能是解释变量是否取值的虚拟变量。如果第一步回归的被解释变量是原回归中解释变量是否取值的虚拟变量，那么该模型就不再是样本选择模型了，而变成了下文将要介绍的处理效应模型，关于这一点，实际应用中经常被搞混。
  • 第三，第一步选择方程的解释变量必须要包括原回归中所有解释变量和至少一个外生变量，也就是说，原回归的解释变量是选择方程解释变量的真子集。如果只使用原回归中一部分的解释变量或不引入外生变量，那么就不能确保IMR与原回归的随机干扰项不相关，从而造成估计系数依然存在偏误。实际应用中，多数文献并未引入外生变量，部分文献甚至没有汇报第一步选择方程中的解释变量，这样的做法十分不推荐。此外，论文中如果引入了外生变量，就需要对相关性与外生性进行具体说明，其中相关性不能只从外生变量的回归系数显著这一个方面进行说明，还要从其他文献和从理论上进行分析；外生性的说明与之类似。
  • 第四，第一步选择方程只能使用probit模型进行回归，不能使用logit模型。在选择方程中，假设扰动项服从正态分布，从而可以推导出将IMR代入原回归方程可以缓解样本选择偏差问题，因此对于被解释变量为0-1型的虚拟变量，只能使用probit模型而不能使用logit模型，因为logit模型不具有扰动项服从正态分布的假设。但问题是，probit假设时间效应和个体效应与扰动项不相关，即第一步选择方程中只能使用随机效应模型，不能使用更一般化的固定效应模型。实际应用中，多数文献在汇报第一阶段回归结果时，在末尾加上“时间固定效应 – Yes”、“个体固定效应 – Yes”等，这样的做法是有待商榷的，因为这根本就不是固定效应模型。
• 第二步 ：将第一步回归计算得到的IMR作为控制变量引入原回归方程中。如果IMR显著，说明原回归中存在样本选择偏差，需要使用样本选择模型进行缓解，而其余变量的回归系数则是缓解样本选择偏差后更为稳健的结果；如果IMR不显著，说明原回归存在的样本选择偏差问题不是很严重，不需要使用样本选择模型，当然，使用了也没关系，因为引入控制变量的回归结果可以与原回归结果比较，作为一种形式的稳健性检验。这里有两点需要注意：
  • 第一，两步估计法中第二步回归代入的是第一步回归的结果，因此第一步回归的估计误差也将被代入第二步，造成效率损失，最终导致第二步估计系数的标准误存在偏差，影响p值进而影响系数显著性。解决方法有两种：一是对第二步回归的标准误进行校正处理，但标准误的校正方法相对复杂，因此现阶段采用这种解决方案的文献几乎没有；二是使用极大似然估计（Maximum Likelihood Estimate，MLE），直接对两阶段回归进行整体估计，这种方法在实际应用中使用较多，但存在的问题在于如果样本量太大，计算会非常耗时。因此，考虑到操作的简便性、理解的直观性以及对分布的假设更为宽松，目前国内流行使用的还是两步估计法。
  • 第二，第二步回归使用的样本数目少于第一步。假设所有的解释变量（包括第一步的外生变量）都没有缺失值，仅被解释变量y存在缺失值，那么第一步回归中使用的样本数目是全样本，因为第一步选择方程的被解释变量y_dummy设置为当y取值不为空（包括y取值为0）时y_dummy等于1，y取值为空时y_dummy等于0，故所有样本的y_dummy都有取值，因此都参与了第一步回归。而第二步回归中的被解释变量y存在缺失值，存在缺失值的样本在参与回归时将直接被剔除。因此第二步回归使用的样本数目少于第一步，这也是样本选择模型一个最直观的特征，这与下文介绍的处理效应模型形成比较。

2.2 处理效应模型

对于自选择偏差导致的估计偏误，将使用处理效应模型（Treatment Effects Model）来缓解，该模型由Maddala（1983）提出。

[4] Maddala G S. Limited-Dependent and Qualitative Variables in Econometrics[M]. Cambridge University Press, 1986.

事实上，使用处理效应模型也只是一定程度上缓解自选择偏差问题。正如『上期推文的1.3小节』所论述的，决定个体是否参与实验的因素可以分为两种：

• 一种是可观测因素，如果个体参与实验的决策依赖于可观测因素，就说明该个体的决策依可测变量选择。
• 另一种是不可观测因素，如果个体参与实验的决策依赖于不可观测因素，就说明该个体的决策依不可测变量选择。

相应地，解决自选择偏差问题的方法也大致可以分为两类：

• 解决依可测变量选择问题的方法如上期介绍的PSM，通过控制处理组与控制组协变量的取值大致相等，从而达到变量选择近似随机的目的。
• 解决依不可测变量选择问题的方法包括PSM - DID方法、断点回归方法（RDD）以及这里的处理效应模型等。需要注意的是，单纯的PSM只能解决依可测变量选择的内生问题，而将PSM和DID结合（即PSM - DID）就可以缓解一部分由不可观测因素带来的自选择偏差问题。

处理效应模型的构建基于Heckman两步法的思想，但与Heckman两步法或者样本选择模型有着本质上的区别，最明显的区别在于，样本选择模型第一阶段回归的被解释变量是第二阶段被解释变量y是否取值的虚拟变量y_dummy，并且y_dummy不参与第二阶段回归；而处理效应模型第一阶段回归的被解释变量是第二阶段的核心解释变量D，并且D的取值为0或1，不存在缺失值。

同样，自选择偏差本质上也是一个因遗漏变量而导致的内生性问题，被遗漏的变量也是IMR，但其计算公式与样本选择偏差存在区别。具体而言，存在自选择偏差的回归方程中被遗漏的IMR计算公式为：

$$IM{R}_i(或{\lambda }_i)=\{\begin{array}{llc}\frac{\varphi ({\hat{y}}_i)}{\Phi ({\hat{y}}_i)}& ,& if{D}_i=1\\ \frac{-\varphi ({\hat{y}}_i)}{\Phi (-{\hat{y}}_i)}& ,& if{D}_i=0\end{array}\text{\hspace{1em}(2)}$$

上式各字母的解释同公式$(1)$。关于IMR计算公式的更多细节，请参考Stata官方网站的回答（FAQs）。

明显可以看到，公式$(1)$说明在样本选择模型中，所有样本的IMR均用一个公式来计算；公式$(2)$说明在处理效应模型中，D取值为1的样本与D取值为0的样本的IMR计算公式不同，而且由于处理效应模型第二阶段回归中所有样本均参与了回归，因此如果混用了计量模型将直接导致变量IMR的取值错误，进而影响第二步回归的估计结果。

同样，处理效应模型的估计思路是：首先，计算全部样本的IMR；随后，将遗漏变量IMR代入原回归方程中，具体来说：

• 第一步 ：使用probit模型估计选择方程，其中选择方程的被解释变量是第二步回归中的核心解释变量D，该解释变量为虚拟变量且不存在缺失值；选择方程的解释变量包括由第二阶段回归中所有解释变量组成的控制变量集以及一个或多个外生变量组成的工具变量集Z，这里之所以直接说Z是工具变量，是因为要求Z满足相关性与外生性，而相关性说的是Z与原回归方程中的解释变量D相关，而非样本选择模型中的要求外生变量与y_dummy相关。同样，回归模型只能使用probit方法，此外也不能使用固定效应模型，在汇报时只能说是“个体效应 – Yes”或“时间效应 – Yes”。
  • 需要注意的是，选择方程中的工具变量应尽量避免使用D的滞后项D_lag，原因在于如果是普通DID，对于所有处理组来说政策实施时点都是一致的，那么在第一步回归中，D_lag会因为多重共线性而被omitted；如果是多期DID，尽管政策实施时点不固定，但总共的实施时点必然不会过多，D_lag同样也会因为多重共线性而被omitted。而对于非DID的D而言，滞后项D_lag则有可能作为一个良好的工具变量，下文第六部分『公开数据的Stata实操』就是一个非DID的例子。
• 第二步 ：将样本数据代入第一步选择方程中，得到各个样本的的拟合值${\hat{y}}_i$，再将${\hat{y}}_i$代入处理效应模型的风险函数（公式$(2)$）中，计算得到各样本的IMR，最后将IMR作为额外的控制变量引入原回归方程中，考察核心解释变量D以及IMR的估计系数。如果IMR的估计系数显著，说明自选择偏差问题不可忽视，此时核心解释变量D的系数就是考虑了自选择偏差后的估计结果，并可与基准回归结果对比构成稳健性检验；而如果IMR的估计系数不显著，则说明自选择偏差问题在原回归中不明显，基准回归结果本身就是可信的。
  • 需要注意的是，核心解释变量D在两步模型中均参与了回归，其中第一阶段回归中D作为被解释变量，在第二阶段回归中作为解释变量，并且我们假设D不存在缺失值，因此处理效应模型两步回归中的样本均是全样本，这不同于样本选择模型。

2.3 估计思路的对比

总结一下样本选择模型和处理效应模型的估计思路的异同点。

相同点在于：

1. 都是两步估计法。Heckman于1979年提出的两步估计法最开始是用于解决样本选择偏差的，即最初的Heckman两步法指的就是样本选择模型，后来有学者借鉴这种两步估计法的思想，应用于解决自选择偏差的处理效应模型。这两个模型在估计思路上是一脉相承的，而正是因为这种相似性，所以才导致各个学者对这两个模型的错误理解与错误应用，这种错误在现阶段的文献中较为常见。
2. 都可以使用MLE进行模型的整体估计。两步估计法（如2SLS、PSM – DID以及这里的样本选择模型和处理效应模型等）一个明显的缺陷是，第一步估计的误差将被带入第二步，导致效率损失。而使用MLE从整体上进行参数估计可以避免这种问题，但如果样本量过大，MLE估计耗时较长，且MLE对分布的假设较为严格，因此需要在估计的精准性、操作的简便性等方面进行权衡。
3. 第一阶段回归都需要引入外生变量，同时应包括第二阶段的所有外生解释变量。引入的外生变量需满足相关性和外生性的要求，即与选择方程中的被解释变量在理论上和统计上均具有相关性，而与第二步回归的被解释变量不具有直接的相关关系。引入外生变量的目的是确保第一步计算得到的IMR在引入原回归方程后不与干扰项相关。该外生变量在处理效应模型中可以直接称作工具变量。此外，如果核心解释变量D是DID模型的did项，那么为了防止出现多重共线性，应该尽量避免使用D的滞后项D_lag作为工具变量。事实上，如果找到了一个良好的工具变量，也完全能够使用2SLS解决内生性问题。此外，两个模型除了都需要在第一阶段引入至少一个外生变量，第一阶段回归中的其余控制变量也应该是第二阶段回归中所有的控制变量，即应该包括所有的外生解释变量，原因在于保证两阶段估计的一致性，详情请看陈强教授的推文『工具变量法（五）: 为何第一阶段回归应包括所有外生解释变量』。然而，部分文献在第一阶段并未包括第二阶段所有的外生解释变量，少部分文献甚至根本就不引入第二阶段的外生解释变量（如，考虑滞后效应，直接引入第二阶段外生解释变量的滞后项），并且在Stata处理效应模型的官方命令etregress的help文件的演示案例中，第一阶段回归也并未包括所有的外生解释变量，原因可能在于IMR是一个非线性项，因此不包含所有外生解释变量引起的内生性问题可能并没有2SLS那么严重。
4. 第一步回归都只能是probit模型。由于logit模型不具备扰动项服从正态分布的假设，如果使用logit模型估计选择方程，将直接导致IMR计算错误，因为Heckman（1979）在推导IMR时，假设选择方程的随机扰动项服从正态分布。这与PSM不同，PSM估计概率方程可以使用logit模型，也可以使用probit模型，并且实际使用中流行的是logit模型。然而，选择方程使用probit模型进行估计有一个问题不可忽视，那就是probit（包括Stata的xtprobit）不能估计固定效应模型，因此即便在回归方程中引入时间虚拟变量和个体虚拟变量，控制的也只是“时间效应”和“个体效应”，不能加入“固定”二字。

不同点在于：

1. 解决的问题不同。样本选择模型解决的是样本选择偏差导致的内生性问题，处理效应模型解决（或者“缓解”）的是依不可观测因素导致的自选择偏差问题。在实际应用中，部分文献在分析内生性问题时将样本选择偏差与自选择偏差混淆，从而使用的模型也是不恰当的。在数据搜集过程中，对被解释变量存在缺失值的样本，多数文献的做法是直接把这些样本剔除，因而即便文章中考虑到了样本选择偏差问题，我们也无法使用样本选择模型（或Heckman两步法）。事实上，囿于数据缺陷，大多数实证类论文都不具备实施Heckman两步法的条件。对于DID类的实证论文，对内生性的分析角度应该更多考虑从自选择偏差切入，而非样本选择偏差，因为各样本处理组虚拟变量D的取值本身就提供了自选择偏差分析的条件，即D取值为1的样本与D取值为0的样本在某些方面是否存在明显的特征差异？或者，是否存在某些因素影响了各样本是否实施政策的决定，而这些因素在两组间又是否存在巨大差异？同时，这些因素是否在理论与统计意义上影响我们想研究的经济指标？在这样的分析之后，就可以使用处理效应模型来缓解因自选择偏差而导致的估计偏误。
2. 变量的设置不同。在样本选择模型第一阶段回归方程中，被解释变量是原方程中的被解释变量y是否被观测到的虚拟变量y_dummy，该变量不参与第二阶段回归，同时第一阶段引入的外生变量直接影响的是y_dummy。在处理效应模型第一阶段回归方程中，被解释变量是原方程的核心解释变量D，D取值为0或1，且不存在缺失值，该变量还同时参与了第二阶段回归，此外第一阶段引入的外生变量（或称工具变量）直接影响的是D。
3. 各阶段样本参与回归的数目不同。假设除关键变量，其余变量都不存在缺失值，那么对于样本选择模型来说，第一阶段回归的解释变量均不存在缺失值，被解释变量y_dummy取值为0或1，也不存在缺失值，因此选择方程中参与回归的样本是全样本，第二阶段由于被解释变量y本身就存在缺失值，因此参与第二阶段回归的样本不是全样本，从而第一阶段的样本多于第二阶段。对于处理效应模型来说，所有变量均不存在缺失值，因此两阶段参与回归的样本是相同的，虽然在第一阶段引入滞后项D_lag作为工具变量的情况下会损失一部分样本，但由于计算出来的IMR同样也存在缺失值，从而第二阶段参与回归的样本也将与第一阶段相同。
4. IMR的计算公式不同。从公式$(1)$和公式$(2)$就可以看出，对于样本选择模型，各样本的IMR计算公式相同；对于处理效应模型来说，D取值为1的样本和D取值为0的样本IMR计算公式并不相同，并且所有样本的IMR均参与了第二步回归。所以，如果混淆了样本选择模型和处理效应模型，将直接导致变量IMR的计算错误，反而进一步造成了估计偏误。

下面推文的第三、第四部分将分别使用示例数据演示样本选择模型和处理效应模型在Stata中的规范操作。

三、样本选择模型的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. R-squared = 0.2535 -------------+---------------------------------- Adj R-squared = 0.2524 Total | 53354.8946 1,342 39. Root MSE = 5.452 ------------------------------------------------------------------------------ wage | Coef. Std. Err. t P>|t| [95% Conf. Interval] -------------+---------------------------------------------------------------- education | . .0 18.00 0.000 . . age | . .0 7.83 0.000 . . _cons | 6.084875 . 6.84 0.000 4. 7. ------------------------------------------------------------------------------ */ 

第二个是样本选择模型，使用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 | . .0 18.59 0.000 . 1.094334 age | . .0 10.34 0.000 . . _cons | . 1.077037 0.45 0.652 -1. 2.59673 -------------+---------------------------------------------------------------- wage_dummy | education | .0 .0 5.19 0.000 .0 .0 age | .0 .0041533 8.79 0.000 .0 .0 married | . .0 6.61 0.000 . . children | . .0 15.79 0.000 . . _cons | -2. . -13.16 0.000 -2. -2. -------------+---------------------------------------------------------------- /athrho | . . 8.62 0.000 . 1.072993 /lnsigma | 1. .027598 64.95 0.000 1. 1.84665 -------------+---------------------------------------------------------------- rho | . .0 . . sigma | 6.004797 . 5.68862 6. lambda | 4. . 3. 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 | . .0 18.53 0.000 . 1.094643 age | . .0 10.10 0.000 . . _cons | . 1.099121 0.44 0.659 -1. 2. -------------+---------------------------------------------------------------- wage_dummy | education | .0 .0 5.12 0.000 .034388 .0 age | .0 .0042243 8.64 0.000 .0 .0 married | . .0 6.66 0.000 . . children | . .0 16.08 0.000 . . _cons | -2. . -13.22 0.000 -2. -2. -------------+---------------------------------------------------------------- /athrho | . . 8.32 0.000 . 1.080266 /lnsigma | 1. .0 62.17 0.000 1.73605 1. -------------+---------------------------------------------------------------- rho | . .0 . . sigma | 6.004797 . 5. 6. lambda | 4. . 3. 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 | . .0 18.23 0.000 . 1.088133 age | . .0 9.61 0.000 . . _cons | . 1. 0.59 0.557 -1. 3. -------------+---------------------------------------------------------------- wage_dummy | education | .0 .0 5.32 0.000 .0 .0 age | .0 .0042293 8.21 0.000 .0 .0 married | . .074208 5.81 0.000 . . children | . .0 15.56 0.000 . . _cons | -2. . -12.81 0.000 -2. -2.089948 -------------+---------------------------------------------------------------- /mills | lambda | 4.001615 . 6.60 0.000 2. 5.19041 -------------+---------------------------------------------------------------- rho | 0.67284 sigma | 5. ------------------------------------------------------------------------------ */ 

第五个是手工完成两步估计法，但实际使用时这种方法十分不推荐（包括上文使用两步法估计的样本选择模型），因为两步估计法下第一步回归的估计误差会带到第二步，造成效率损失，影响第二步估计系数的显著性。

手工两步估计法的思路是先用控制变量和外生变量对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 | .0 .0 5.32 0.000 .0 .0 age | .0 .0042293 8.21 0.000 .0 .0 married | . .074208 5.81 0.000 . . children | . .0 15.56 0.000 . . _cons | -2. . -12.81 0.000 -2. -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. R-squared = 0.2793 -------------+---------------------------------- Adj R-squared = 0.2777 Total | 53354.8946 1,342 39. Root MSE = 5.3587 ------------------------------------------------------------------------------ wage | Coef. Std. Err. t P>|t| [95% Conf. Interval] -------------+---------------------------------------------------------------- education | . .0 19.46 0.000 . 1.08159 age | . .0 10.25 0.000 . . imr | 4.001616 . 6.93 0.000 2. 5. _cons | . 1. 0.63 0.529 -1. 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）的回归系数的符号、大小、显著性在模型$(2)(3)(4)(6)$中基本保持一致，并且与模型$(1)$的区别不大，这说明在考虑到样本选择偏差后，根据基本回归结果得出的结论依然成立。
• 模型$(2)(3)$缺少变量lambda，在实际汇报时应该根据原始回归表进行汇报。
• 模型$(2)(3)(4)$汇报的样本量是2,000个全样本，但应该清楚不同回归阶段实际参与的样本数目是不同的。
• 关于模型拟合程度的统计量，模型$(1)$和模型$(6)$汇报的都是修正R方，模型$(5)$的是伪R方，模型$(2)(3)(4)$没有汇报统计量，但观察原始回归表可以发现，模型$(2)(3)(4)$均在左上角汇报了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. Prob > F = 0.0000 Residual | 4196.70492 1,203 3. R-squared = 0.3402 -------------+---------------------------------- Adj R-squared = 0.3369 Total | 6360.19947 1,209 5. Root MSE = 1.8678 ------------------------------------------------------------------------------ wage | Coef. Std. Err. t P>|t| [95% Conf. Interval] -------------+---------------------------------------------------------------- union | 1.003913 . 7.50 0.000 . 1. age | . .0 7.56 0.000 . . grade | . .0 14.82 0.000 . . smsa | . . 7.79 0.000 . 1. black | -. .1252 -4.94 0.000 -. -. tenure | . .0 6.26 0.000 . . _cons | -4. . -8.14 0.000 -5. -3. ------------------------------------------------------------------------------ */ 

第二个是处理效应模型，使用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 | . .0 7.70 0.000 . . grade | . .0 14.33 0.000 . . smsa | . . 7.30 0.000 . 1. black | -. . -5.77 0.000 -1.05619 -. tenure | . .0 4.12 0.000 .0 . 1.union | 2. . 10.71 0.000 2.4069 3. _cons | -4. . -8.24 0.000 -5. -3. -------------+---------------------------------------------------------------- union | black | . .0 4.76 0.000 . . tenure | .0 .0 3.75 0.000 .0 . south | -. .0 -6.82 0.000 -. -. _cons | -. .0 -12.22 0.000 -1.027576 -. -------------+---------------------------------------------------------------- /athrho | -. .0 -7.19 0.000 -. -. /lnsigma | . .0 23.95 0.000 . . -------------+---------------------------------------------------------------- rho | -. .060971 -. -. sigma | 2.019151 .0 1. 2. lambda | -1.1603 . -1. -. ------------------------------------------------------------------------------ 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 | . .0 7.19 0.000 . . grade | . .0 11.13 0.000 . . smsa | . . 7.60 0.000 . 1. black | -. . -5.59 0.000 -1.064404 -. tenure | . .0 3.22 0.001 .0 . 1.union | 2. . 5.22 0.000 1. 4.051988 _cons | -4. . -6.83 0.000 -5. -3. -------------+---------------------------------------------------------------- union | black | . .093489 4.87 0.000 . . tenure | .0 .0 3.46 0.001 .0 . south | -. .0 -6.93 0.000 -. -. _cons | -. .0 -11.97 0.000 -1.030522 -. -------------+---------------------------------------------------------------- /athrho | -. . -2.82 0.005 -1. -. /lnsigma | . .0 9.27 0.000 . . -------------+---------------------------------------------------------------- rho | -. . -. -. sigma | 2.019151 . 1. 2. lambda | -1.1603 . -1. -. ------------------------------------------------------------------------------ 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 | . .0 7.92 0.000 . . grade | . .029014 14.56 0.000 . . smsa | . . 6.71 0.000 . 1. black | -. . -5.19 0.000 -1. -. tenure | . .051879 1.93 0.053 -.0013584 . union | 4. 1.006459 4.53 0.000 2. 6. _cons | -4. . -8.65 0.000 -5.72903 -3. -------------+---------------------------------------------------------------- union | black | . .0 4.52 0.000 . . tenure | .0 .0 4.22 0.000 .053396 . south | -. .0 -5.24 0.000 -. -. _cons | -. .0 -12.97 0.000 -1. -. -------------+---------------------------------------------------------------- hazard | lambda | -2.093313 . -3.61 0.000 -3. -. -------------+---------------------------------------------------------------- rho | -0.89172 sigma | 2. ------------------------------------------------------------------------------ */ 

第五个是手工完成两步估计法，估计结果与直接使用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 | . .0 4.52 0.000 . . tenure | .0 .0 4.22 0.000 .053396 . south | -. .0 -5.24 0.000 -. -. _cons | -. .0 -12.97 0.000 -1. -. ------------------------------------------------------------------------------ */ *- 计算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. Prob > F = 0.0000 Residual | 4132.88562 1,202 3. R-squared = 0.3502 -------------+---------------------------------- Adj R-squared = 0.3464 Total | 6360.19947 1,209 5. Root MSE = 1.8543 ------------------------------------------------------------------------------ wage | Coef. Std. Err. t P>|t| [95% Conf. Interval] -------------+---------------------------------------------------------------- union | 4.56386 . 5.45 0.000 2. 6. age | . .0 7.94 0.000 . . grade | . .029448 14.35 0.000 . . smsa | . .12708 6.79 0.000 . 1. black | -. . -6.45 0.000 -1. -. tenure | . .0 2.37 0.018 .0 . imr | -2.093314 . -4.31 0.000 -3.046589 -1. _cons | -4. . -8.75 0.000 -5. -3. ------------------------------------------------------------------------------ */ 

最后是各个回归模型下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 */ 

五、已有文献的简单评述

Lennox et al.（2012）总结了样本选择模型和处理效应模型在会计学领域中应用的几点问题，这几个问题也完全可以延申至其他领域。

• 第一，选择方程中没有引入排他性约束（Exclusion Restrictions）变量，也即上文所说的外生解释变量或工具变量；就算引入了外生变量，也没有对外生变量的相关性和外生性进行详细说明。
• 第二，没有汇报第一阶段的回归结果，因而无法判断是否包含外生变量，更无法判断外生变量的引入是否有效。

5.1 常见问题

除了以上几点，现阶段文献中经常出现的问题还有：

• 第一，混淆样本选择偏差和自选择偏差。如前文所述，样本选择偏差和自选择偏差有着本质上的不同，最关键的不同在于非随机选择的机制是不同的。由于数据缺陷，多数文献出现的内生性问题其实并非由样本选择偏差所导致，就算怀疑存在样本选择偏差，也囿于数据而无法进一步实施样本选择模型的回归操作。而论文中常见的做法是，通篇说的样本选择偏差，但实际出现的问题在于处理变量D的取值不随机，也即本质上是自选择偏差；或者，样本选择偏差与自选择偏差这两种说法混着用。这样的做法也确实给读者造成了困扰，笔者之前阅读过多篇论文，确实感觉不同论文这对方面的说明大相径庭。而正是由于对问题的界定不清晰，所以造成了对模型使用的偏误。
• 第二，混淆样本选择模型和处理效应模型。现阶段文献对样本选择模型和处理效应模型的实际操作主要有三种：heckman、etregress和手工两步法，其中heckman和手工两步法用的最多，etregress几乎没有。但是，正如前面所说的，样本选择模型由于数据缺陷在应用中不具有可操作性，因此heckman在多数文献中本意是解决自选择偏差，因此部分文献在heckman选择方程中设置被解释变量为D而非y_dummy就不足为奇了，但这样的做法是有问题的，除了模型使用的混淆，更重要的在于两种模型在前后两步中参与回归的样本数目是不同的，并且两种模型在计算IMR时使用的公式也不同，因此模型混用可能反而进一步导致估计偏误。此外，无论是heckman还是etregress，以及手工操作，两步估计法都不推荐使用，比较而言，MLE得出的结果更稳健。

5.2 经典论文解析

下面将解析两篇经典论文。考虑到多数文献本质上存在的是自选择偏差，而分析样本选择偏差并使用样本选择模型的（中文）文献并不多见，因此下面两篇论文均是从样本选择偏差角度来分析内生性问题。

第一篇论文是陈云松（2012）发表在《社会》的《农民工收入与村庄网络 基于多重模型识别策略的因果效应分析》，研究主题是探讨社会网络对农民工在外务工收入的影响。

[8] 陈云松. 农民工收入与村庄网络 基于多重模型识别策略的因果效应分析[J]. 社会, 2012, 32(04): 68-92.

以下是基准回归方程：

$$W_{ig}={\beta }_0+{\beta }_1{S}_g+{\beta }_2{X}_{ig}+{\beta }_3{V}_g+\epsilon \text{\hspace{1em}(3)}$$

其中，$W_{ig}$表示第$g$个村庄第$i$个农民工在城市务工的工资；$S_g$代表社会网络，用村庄在外务工的人数来表示，在论文中作为核心解释变量；$X_{ig}$表示个人层面的控制变量；$V_g$是村庄层面的控制变量。

文章怀疑农民工在外务工可能是一个选择行为，因为具有城市劳动力市场优势（男性、年轻和能力强等）的农民会更倾向于外出打工。论文进一步将这种选择性的来源分为可观测因素和不可观测因素，可观测因素包括年龄、性别等，不可观测因素包括性格、能力等。因此，在模型设置时必须要考虑样本群体是否随机和均质，即样本选择偏差问题。

在本例的数据集结构中，被解释变量$W_{ig}$存在缺失值，而由于不在外务工的农民本身就不具有在外务工的工资数据（即使他们有其他来源的收入，但不是文章研究的重点），因此这些缺失值存在的原因就是这些被调查的样本本身就不在外务工。而由于前面提到的农民是否在外务工可能是一个选择行为，即$W_{ig}$存在缺失值的样本与$W_{ig}$取值不为空的样本在某些特征因素方面本身就存在较大的差异，因此如果在回归时直接剔除这部分取值为空的样本，最后得到的结果就可能存在估计偏误，也就是说在考虑到样本选择偏差的情况下，基准回归结果可能就不再具有稳健性。

为了解决可能存在的样本选择偏差问题，作者使用了样本选择模型中的两步估计法（Heckit）。Heckit由以下方程组构成：

$$W_{ig}={\beta }_0+{\beta }_1{S}_g+{\beta }_2{X}_{ig}+{\beta }_3{V}_g+{\beta }_4{\hat{P}}_{ig}+\epsilon \text{\hspace{1em}(4)}$$

$$P_{ig}={\gamma }_0+{\gamma }_1{F}_{ig}+{\gamma }_2{S}_g+{\gamma }_3{X}_{ig}+{\gamma }_4{V}_g+\mu \text{\hspace{1em}(5)}$$

其中，方程$(5)$是第一阶段回归方程（选择方程），方程$(4)$是第二阶段回归方程；${\hat{P}}_{ig}$是逆米尔斯比率；$P_{ig}$代表样本是否外出务工的虚拟变量，即样本在外务工取值为1，否则为0；$F_{ig}$是选择方程中的外生解释变量，论文中选择的是家庭劳动力人数；方程$(4)$的所有解释变量，是方程$(5)$解释变量的严格子集（真子集）。

论文还对外生变量$F_{ig}$的相关性与外生性进行了说明，认为家庭劳动力数量对农民工的打工决策有着重要影响，而对在外务工收入的影响微乎其微，具体分析请看原文。

事实上，这里存在两个问题：

• 第一，在选择方程中，作者使用logit模型进行回归，正如前文所述，logit模型不具备干扰项服从正态分布的假设，因此根据第一步回归拟合值计算出的IMR可能存在一定程度的偏误。
• 第二，方程$(4)$和方程$(5)$对${\hat{P}}_{ig}$或$P_{ig}$的界定不清晰，因为方程$(4)$中的${\hat{P}}_{ig}$本意应该是IMR，而非字母本身所表达的“$P_{ig}$的拟合值”的含义，而IMR正是由“$P_{ig}$的拟合值”计算所得。

下表是基准OLS回归与Heckit第二步回归结果的对比，括号内为各变量估计系数的聚类稳健标准误，以村庄为聚类单位；Heckit的第一步回归结果论文并未列示，这里假定外生解释变量在第一步回归中显著且有效。

可以观察到，Heckit第二步回归结果中IMR显著为正，且数值较大（相较于其他控制变量的估计系数而言，详细结果请看原文），这说明基准OLS回归确实存在样本选择偏差，造成估计偏误，具体来说是低估了社会网络对农民工在外务工收入的影响，因为OLS模型中社会网络的估计系数仅有0.125，而Heckit模型的估计系数（0.263）是其两倍还多，且两者均至少在1%的水平下显著。至于为什么社会网络的收入促进效应在Heckit模型中高于OLS，作者在原文中给出了解释。

值得一提的是，为了解决一般性的因遗漏变量和联立方程（互为因果）导致的估计偏误问题，作者在样本选择模型的基础上进一步采用工具变量法，即采用IV - Heckit对基准OLS的稳健性进行进一步的检验，详情请看原文。

第二篇论文是祝树金和赵玉龙（2017）发表在《金融研究》的《资源错配与企业的出口行为——基于中国工业企业数据的经验研究》，主题是探讨企业资源错配对出口行为的影响。

[9] 祝树金, 赵玉龙. 资源错配与企业的出口行为——基于中国工业企业数据的经验研究[J]. 金融研究, 2017(11): 49-64.

论文一开始就考虑到了样本选择偏差问题，认为企业是否出口受制于自身条件，简单将出口企业与非出口企业同等对待将产生估计偏误，因此构建样本选择模型对这种偏误进行纠正。构建的两阶段模型如下：

$$exdu{m}_{it}={\alpha }_{1}mis{s}_{it-1}+\sum _1^J{\beta }_j{X}_{j,it-1}+{\delta }_{it}\text{\hspace{1em}(6)}$$

$$exshar{e}_{it}={\alpha }_{2}mis{s}_{it-1}+\sum _1^J{\beta }_j{Z}_{j,it-1}+{\epsilon }_{it}\text{\hspace{1em}(7)}$$

公式$(6)$是第一阶段的企业出口选择方程，其中，被解释变量$exdu{m}_{it}$代表第$t$年$i$企业的出口状态，若有出口行为，该变量记为1，否则记为0；考虑到出口滞后效应，第一、二阶段所有解释变量均滞后一期；$mis{s}_{it-1}$是论文的核心解释变量企业资源错配，分别使用企业资源错配指标$miss1$、$\tau k$和$\tau l$来表示；$X_{j,it-1}$是影响企业出口决策的第$j$个控制变量，包括一个外生解释变量$exdu{m}_{it-1}$，即企业前一期是否出口的虚拟变量，作者认为该变量满足相关性和外生性的要求，以及其他控制变量。

公式$(7)$是第二阶段的回归方程，其中，被解释变量$exshar{e}_{it}$表示企业出口强度；$Z_{j,it-1}$表示影响企业出口强度的控制变量，这些控制变量包括第一阶段的所有控制变量（除$exdu{m}_{it-1}$），以及一个根据第一阶段回归拟合值计算的IMR。

回归结果汇总如下（限于篇幅限制，这里仅展示核心解释变量为$miss1$的回归结果；括号内为标准误，具体类型未告知）：

可以发现（$Waldch{i}^2$和$\rho$的那两列结果应该合并为一列，markdown表格合并比较麻烦~），第一步选择方程（模型$(1)$）中，外生变量的估计系数显著为正，说明外生变量的选择有效；第二阶段回归（模型$(2)$）中，核心解释变量的估计系数显著为正，说明在考虑样本选择偏差的情况下，企业资源错配仍对企业出口强度产生促进作用。需要注意的是：

• 模型$(1)$和模型$(2)$引入行业虚拟变量和年份虚拟变量，并在结果汇报中称之为行业（效应）和年份（效应），而非行业固定效应与年份固定效应，这种做法是比较严谨的，因为第一阶段的probit回归不能使用固定效应模型。
• 模型$(1)$和模型$(2)$使用的样本数是不同的，模型$(1)$使用的样本数应该多于模型$(2)$。然而，结果显示两步回归均使用了177,386个样本，因此猜测论文中的样本选择模型是使用MLE方法进行估计的，因为MLE从模型整体角度进行参数估计，参与回归的样本一般汇报为样本总数。
• 模型$(2)$没有汇报IMR的估计系数，而仅仅汇报两步回归方程干扰项之间的相关系数$\rho$,因此可以基本断定该样本选择模型就是使用MLE方法进行估计的。$\rho$在1%的水平下显著为负，说明模型中存在的样本选择偏差不能忽视。

六、公开数据的Stata实操

以上论文解析都没有使用数据复现，下面将使用一篇论文的公开数据与代码简单演示一下样本选择模型与处理效应模型在Stata中的操作。

需要提前说明的是：

• 这篇论文出现的问题本质上是自选择偏差，作者误认为是样本选择偏差，并错误地使用了样本选择模型。
• 由于在Stata的etregress命令下，第一步选择方程的被解释变量D只能是0-1型的虚拟变量，如果使用连续型变量，Stata将会报错。
• 在D为连续型变量的情况下，为了使得Stata不报错，选择方程使用D是否取值的虚拟变量，第二步方程直接引入D作为核心解释变量，这样第二步回归模型中将会同时出现D以及D是否取值的虚拟变量，这种情况下可能产生多重共线性问题。
• 为了避免出现多重共线性问题，公开代码使用heckman命令而非etregress命令，因为heckman命令下第一步回归的被解释变量不会自动带入第二步回归，这或许也是论文作者选择heckman而非etregress的动机。但需要说明的是，这样的做法本质上是不严谨的，不同的模型适用于不同的问题，模型混用可能导致估计结果的严重偏误，虽然最后得到的结果可能符合预期，但不严谨的实证设计很难说服读者接受其结论。
• 为了避免模型混用导致的估计偏误，推文中的代码将使用手工两步法进行处理效应模型的估计。

这篇论文是杜勇等（2021）发表在《中国工业经济》的《共同机构所有权与企业盈余管理》，研究主题是探讨共同机构持股对企业盈余管理的影响。

[10] 杜勇, 孙帆, 邓旭. 共同机构所有权与企业盈余管理[J]. 中国工业经济, 2021(06): 155-173.

由于上文已对两个命令进行了足够详细的讲解，并且对相关问题进行了提前说明，因此下面直接放出代码与结果，代码与结果的具体解读可以参考前文。

# 【数据来源】 *- 杜勇等（2021），参见在《中国工业经济》网站（http://ciejournal.ajcass.org/Magazine/show/?id=77795） # 【参考文献】 *- 杜勇, 孙帆, 邓旭. 共同机构所有权与企业盈余管理[J]. 中国工业经济, 2021(06): 155-173. *- Stata Version: 16 | 17 *- 定义路径 cd "C:\Users\KEMOSABE\Desktop\heckman" use 数据.dta, clear *- 连续变量缩尾处理 winsor2 da1 da2 da11 rem /// coz2 coz3 coz_power coz_number /// coz22 coz33 coz2_year coz3_year /// director djg auditfees epcm /// net institution size leverage /// roa growth toptenrate independent /// magpay boardshare invrec analyst, /// cut(1 99) replace *- 设置控制变量的全局暂元 global CV institution size leverage roa /// growth toptenrate dual independent /// magpay boardshare invrec analyst /// opin aud *- 将控制变量及核心解释变量滞后一期 xtset code year foreach i of global CV { gen Lag`i' = L.`i' } gen Lagcoz1 = L.coz1 *- 设置第一阶段回归控制变量的全局暂元 global CV1 Laginstitution Lagsize Lagleverage Lagroa /// Laggrowth Lagtoptenrate Lagdual Lagindependent /// Lagmagpay Lagboardshare Laginvrec Laganalyst /// Lagopin Lagaud # 一、coz1的基准回归及稳健性检验 *- （1）基准回归 qui: reg da1 coz1 $CV i.year i.industry, r est store OLS1 *- （2）处理效应模型（MLE） qui: etregress da1 $CV i.year i.industry, /// treat(coz1 = Lagcoz1 $CV i.year i.industry) nolog vce(robust) est store et_ML1 /* Linear regression with endogenous treatment Number of obs = 20,018 Estimator: maximum likelihood Wald chi2(31) = 1485.50 Log pseudolikelihood = 18287.002 Prob > chi2 = 0.0000 ------------------------------------------------------------------------------ | Robust | Coef. Std. Err. z P>|z| [95% Conf. Interval] -------------+---------------------------------------------------------------- da1 | institution | .0001648 .0001052 1.57 0.117 -.0000414 .000371 size | -.0066258 .0008814 -7.52 0.000 -.0083534 -.0048982 leverage | .0 .0048317 4.91 0.000 .0 .0 roa | -. .0 -13.25 0.000 -. -. growth | .0 .0023365 15.11 0.000 .0 .0 toptenrate | .0001971 .0000445 4.43 0.000 .0001099 .0002844 dual | .0031043 .0015126 2.05 0.040 .0001396 .006069 independent | .0040687 .0 0.35 0.724 -.0 .0 magpay | . . 2.89 0.004 . . boardshare | -.0090979 .0037856 -2.40 0.016 -.0 -.0016784 invrec | .0 .0051727 3.54 0.000 .0081693 .0 analyst | .0005457 .000082 6.65 0.000 .0003849 .0007065 opin | -.0 .0040999 -9.50 0.000 -.047003 -.0 aud | -.0038764 .0012816 -3.02 0.002 -.0063884 -.0013644 | year | 2009 | -.0008564 .0043254 -0.20 0.843 -.0093341 .0076213 2010 | .0 .0049932 3.73 0.000 .0088619 .0 2011 | .0 .0042826 5.15 0.000 .0 .030431 2012 | -.0 .0037267 -4.26 0.000 -.0 -.0085688 2013 | -.0 .0036 -4.48 0.000 -.0 -.009063 2014 | -.0 .0035979 -4.48 0.000 -.0 -.0090528 2015 | -.0031965 .0037801 -0.85 0.398 -.0 .0042122 2016 | -.0027102 .0037549 -0.72 0.470 -.0 .0046493 2017 | -.0 .0036145 -4.24 0.000 -.0 -.0082243 2018 | -.0 .0036009 -4.30 0.000 -.0 -.0084143 2019 | -.0 .0035699 -4.99 0.000 -.0 -.0 | industry | 2 | -.0 .0 -2.82 0.005 -.0 -.0090968 3 | -.0 .0 -1.06 0.291 -.0 .0099234 4 | -.0 .0 -2.45 0.014 -.0 -.0056462 5 | -.0 .0 -4.17 0.000 -.0 -.0 6 | -.0 .0 -3.01 0.003 -.0 -.0 | 1.coz1 | -.0067129 .0019325 -3.47 0.001 -.0 -.0029253 _cons | . .0 12.94 0.000 . . -------------+---------------------------------------------------------------- coz1 | Lagcoz1 | 3. .041428 78.08 0.000 3.15344 3. institution | .0037066 .0029546 1.25 0.210 -.0020843 .0094976 size | . .0 4.50 0.000 .0 . leverage | -.0 . -0.10 0.918 -. . roa | -. . -1.67 0.095 -1. .0 growth | .083781 .037123 2.26 0.024 .0 . toptenrate | .0032139 .0012742 2.52 0.012 .0007165 .0057114 dual | -.0 .0 -1.61 0.106 -. .0 independent | -1.012848 . -2.77 0.006 -1. -. magpay | 2.53631 4.014366 0.63 0.528 -5. 10.40432 boardshare | -1.021367 . -5.88 0.000 -1. -. invrec | -. . -1.34 0.179 -. .0 analyst | -.0033357 .0025439 -1.31 0.190 -.0083216 .0016502 opin | . . 1.08 0.282 -.0 . aud | .0 .0 0.65 0.515 -.0 . | year | 2009 | .0 . 0.49 0.622 -. . 2010 | . . 1.87 0.061 -.0096534 . 2011 | -.0 . -0.27 0.788 -. . 2012 | .0080566 . 0.07 0.941 -. . 2013 | -.0 . -0.62 0.534 -. . 2014 | -.0 . -0.41 0.679 -. . 2015 | -.0 . -0.51 0.608 -. . 2016 | -.0 . -0.90 0.368 -. . 2017 | -. . -2.33 0.020 -.45645 -.0 2018 | -.0 . -0.92 0.359 -. . 2019 | .0 .0 0.26 0.797 -. . | industry | 2 | -. . -1.09 0.276 -. . 3 | -. . -1.51 0.131 -. .09736 4 | -. . -1.61 0.107 -. .0 5 | .069847 . 0.34 0.731 -. . 6 | .0 . 0.44 0.658 -. . | _cons | -4. . -7.54 0.000 -5. -3. -------------+---------------------------------------------------------------- /athrho | .0 .0 2.07 0.038 .0022294 .0 /lnsigma | -2. .0 -206.51 0.000 -2. -2. -------------+---------------------------------------------------------------- rho | .041339 .0 .0022294 .0 sigma | .0 .0010229 .0 .0 lambda | .0035624 .0017194 .0001924 .0069323 ------------------------------------------------------------------------------ Wald test of indep. eqns. (rho = 0): chi2(1) = 4.29 Prob > chi2 = 0.0383 */ *- （3）处理效应模型（TwoStep） qui: etregress da1 $CV i.year i.industry, /// treat(coz1 = Lagcoz1 $CV i.year i.industry) twostep est store et_2S1 /* Linear regression with endogenous treatment Number of obs = 20018 Estimator: two-step Wald chi2(61) = 3629.46 Prob > chi2 = 0.0000 ------------------------------------------------------------------------------ | Coef. Std. Err. z P>|z| [95% Conf. Interval] -------------+---------------------------------------------------------------- da1 | institution | .0001648 .0000966 1.71 0.088 -.0000245 .0003541 size | -.0066256 .0007952 -8.33 0.000 -.0081841 -.0050671 leverage | .0 .0038194 6.21 0.000 .0 .0 roa | -. .0098235 -22.94 0.000 -. -. growth | .0 .0011144 31.69 0.000 .0 .0 toptenrate | .0001971 .000043 4.58 0.000 .0001128 .0002814 dual | .0031042 .0014823 2.09 0.036 .0001989 .0060095 independent | .0040678 .0 0.36 0.722 -.0 .0 magpay | . . 3.68 0.000 . . boardshare | -.0090989 .0041262 -2.21 0.027 -.017186 -.0010117 invrec | .0 .0041006 4.46 0.000 .0 .0 analyst | .0005457 .0000836 6.53 0.000 .0003818 .0007096 opin | -.0 .0030457 -12.79 0.000 -.0 -.0 aud | -.0038764 .0012695 -3.05 0.002 -.0063645 -.0013882 | year | 2009 | -.0008564 .0038108 -0.22 0.822 -.0083254 .0066127 2010 | .0 .003792 4.92 0.000 .0 .0 2011 | .0 .003709 5.94 0.000 .0 .0 2012 | -.0 .0035715 -4.44 0.000 -.0 -.008873 2013 | -.0 .0035103 -4.59 0.000 -.022999 -.0092388 2014 | -.0 .0034962 -4.61 0.000 -.022957 -.0092523 2015 | -.0031967 .0035138 -0.91 0.363 -.0 .0036903 2016 | -.0027104 .0035086 -0.77 0.440 -.0095872 .0041665 2017 | -.0 .0034914 -4.38 0.000 -.0 -.0084658 2018 | -.0 .0034815 -4.44 0.000 -.0 -.0086485 2019 | -.0 .0034585 -5.15 0.000 -.0 -.0 | industry | 2 | -.0 .0066066 -4.51 0.000 -.0 -.0 3 | -.0 .0069175 -1.68 0.093 -.0 .001942 4 | -.0 .0074022 -3.82 0.000 -.0 -.0 5 | -.0 .0064804 -6.71 0.000 -.0 -.0 6 | -.0 .006855 -4.72 0.000 -.0 -.0 | coz1 | -.0067176 .0021868 -3.07 0.002 -.0 -.0024315 _cons | . .0 15.32 0.000 . . -------------+---------------------------------------------------------------- coz1 | Lagcoz1 | 3. .0 78.33 0.000 3. 3. institution | .0036154 .0029089 1.24 0.214 -.002086 .0093168 size | . .0 4.69 0.000 .0 . leverage | -.0 . -0.17 0.864 -. . roa | -. . -1.73 0.083 -1. .0 growth | .0 .0 2.89 0.004 .0 . toptenrate | .0032304 .0013442 2.40 0.016 .0005958 .0058649 dual | -.079793 .050784 -1.57 0.116 -. .0 independent | -1.019378 . -2.72 0.007 -1. -. magpay | 2.68412 3. 0.73 0.466 -4. 9. boardshare | -1.018715 . -5.87 0.000 -1. -. invrec | -. . -1.35 0.177 -. .0 analyst | -.0033943 .0025893 -1.31 0.190 -.0084692 .0016805 opin | . .0 1.17 0.241 -.0 . aud | .0 .040476 0.64 0.520 -.053281 . | year | 2009 | .0 . 0.45 0.651 -. . 2010 | . . 1.84 0.066 -.0 . 2011 | -.0 . -0.28 0.781 -. . 2012 | .0078514 . 0.07 0.942 -. . 2013 | -.0 . -0.61 0.544 -. . 2014 | -.0 . -0.42 0.674 -. . 2015 | -.0 . -0.53 0.599 -. . 2016 | -.092317 . -0.85 0.395 -. . 2017 | -. . -2.28 0.023 -. -.0 2018 | -.0 . -0.92 0.358 -. . 2019 | .0 . 0.21 0.833 -. . | industry | 2 | -. . -1.16 0.244 -. . 3 | -. . -1.53 0.126 -. .093908 4 | -. . -1.64 0.101 -. .080877 5 | .0 . 0.31 0.756 -. . 6 | .0 . 0.42 0.673 -. . | _cons | -4. . -7.92 0.000 -5. -3. -------------+---------------------------------------------------------------- hazard | lambda | .0035703 .001916 1.86 0.062 -.0001849 .0073256 -------------+---------------------------------------------------------------- rho | 0.04143 sigma | .0 ------------------------------------------------------------------------------ */ *- （4）手工处理（TwoStep） qui: probit coz1 Lagcoz1 $CV i.year i.industry est store First1 predict y1, xb gen imr1 = normalden(y1) / normal(y1) replace imr1 = -normalden(y1) / normal(-y1) if coz1 == 0 qui: reg da1 coz1 $CV i.year i.industry imr1, r est store Second1 *- 结果导出 local mlist1 "OLS1 et_ML1 et_2S1 First1 Second1" esttab `mlist1' using coz1的基准回归及稳健性检验.rtf, replace /// b(%6.4f) t(%6.4f) /// scalar(N r2_a r2_p) compress nogap /// star(* 0.1 0.05 * 0.01) /// indicate("Year FE=*.year" "Industry FE=*.industry") /// mtitle(`mlist1') /// title("Baseline Regression and Robustness Test of coz1") /* Baseline Regression and Robustness Test of coz1 --------------------------------------------------------------------------- (1) (2) (3) (4) (5) LSDV1 et_mle1 et_2s1 first1 second1 --------------------------------------------------------------------------- main coz1 -0.0048* -0.0067* -0.0067* (-2.9786) (-3.0719) (-3.4445) institut~n 0.0002 0.0002 0.0002* 0.0036 0.0002 (1.5850) (1.5665) (1.7065) (1.2429) (1.5654) size -0.0066* -0.0066* -0.0066* 0.1121* -0.0066* (-7.9796) (-7.5170) (-8.3323) (4.6877) (-7.5054) leverage 0.0212* 0.0237* 0.0237* -0.0209 0.0237* (4.5995) (4.9060) (6.2064) (-0.1712) (4.9027) roa -0.2110* -0.2254* -0.2254* -0.5579* -0.2254* (-12.8244) (-13.2518) (-22.9432) (-1.7340) (-13.2449) growth 0.0345* 0.0353* 0.0353* 0.0879* 0.0353* (15.5953) (15.1127) (31.6865) (2.8924) (15.1130) toptenrate 0.0002* 0.0002* 0.0002* 0.0032 0.0002* (5.4173) (4.4277) (4.5837) (2.4032) (4.4238) dual 0.0028* 0.0031 0.0031 -0.0798 0.0031 (1.9385) (2.0523) (2.0941) (-1.5712) (2.0507) independ~t 0.0036 0.0041 0.0041 -1.0194* 0.0041 (0.3212) (0.3527) (0.3558) (-2.7177) (0.3524) magpay 0.3756* 0.4069* 0.4070* 2.6841 0.4070* (2.9158) (2.8879) (3.6835) (0.7287) (2.8857) boardshare -0.0024 -0.0091 -0.0091 -1.0187* -0.0091 (-0.6691) (-2.4033) (-2.2052) (-5.8660) (-2.4013) invrec 0.0230* 0.0183* 0.0183* -0.1787 0.0183* (4.6611) (3.5393) (4.4646) (-1.3510) (3.5365) analyst 0.0006* 0.0005* 0.0005* -0.0034 0.0005* (7.1090) (6.6524) (6.5269) (-1.3109) (6.6451) opin -0.0409* -0.0390* -0.0390* 0.1167 -0.0390* (-10.2467) (-9.5044) (-12.7942) (1.1722) (-9.4991) aud -0.0043* -0.0039* -0.0039* 0.0261 -0.0039* (-3.5288) (-3.0246) (-3.0534) (0.6436) (-3.0220) 1.coz1 -0.0067* (-3.4737) Lagcoz1 3.2352* (78.3295) imr1 0.0036* (1.8696) _cons 0.2681* 0.2775* 0.2775* -4.3443* 0.2775* (13.4216) (12.9383) (15.3159) (-7.9174) (12.9241) --------------------------------------------------------------------------- coz1 Lagcoz1 3.2346* 3.2352* (78.0784) (78.3295) institut~n 0.0037 0.0036 (1.2545) (1.2429) size 0.1103* 0.1121* (4.4971) (4.6877) leverage -0.0127 -0.0209 (-0.1025) (-0.1712) roa -0.5487* -0.5579* (-1.6706) (-1.7340) growth 0.0838 0.0879* (2.2568) (2.8924) toptenrate 0.0032 0.0032 (2.5222) (2.4032) dual -0.0821 -0.0798 (-1.6150) (-1.5712) independ~t -1.0128* -1.0194* (-2.7677) (-2.7177) magpay 2.5363 2.6841 (0.6318) (0.7287) boardshare -1.0214* -1.0187* (-5.8833) (-5.8660) invrec -0.1794 -0.1787 (-1.3444) (-1.3510) analyst -0.0033 -0.0034 (-1.3113) (-1.3109) opin 0.1161 0.1167 (1.0766) (1.1722) aud 0.0256 0.0261 (0.6506) (0.6436) _cons -4.3163* -4.3443* (-7.5393) (-7.9174) --------------------------------------------------------------------------- / athrho 0.0414 (2.0716) lnsigma -2.4514* (-2.1e+02) --------------------------------------------------------------------------- hazard lambda 0.0036* (1.8635) Year FE Yes Yes Yes Yes Yes Industry~E Yes Yes Yes Yes Yes --------------------------------------------------------------------------- N 22591 20018 20018 20018 20018 r2_a 0.1387 0.1441 r2_p 0.6912 --------------------------------------------------------------------------- t statistics in parentheses * p<0.1, p<0.05, * p<0.01 */ 

# 二、coz2的基准回归及稳健性检验 *- （1）基准回归 qui: reg da1 coz2 $CV i.year i.industry, r est store OLS2 *- （2）样本选择模型（TwoStep） qui: heckman da1 coz2 $CV i.year i.industry, /// select( coz2 = $CV1 ) twostep est store he_2S2 /* Heckman selection model -- two-step estimates Number of obs = 20,018 (regression model with sample selection) Selected = 2,593 Nonselected = 17,425 Wald chi2(31) = 559.43 Prob > chi2 = 0.0000 -------------------------------------------------------------------------------- | Coef. Std. Err. z P>|z| [95% Conf. Interval] ---------------+---------------------------------------------------------------- da1 | coz2 | -.0 .0 -1.83 0.068 -.0 .0015293 institution | -.0001099 .0002544 -0.43 0.666 -.0006085 .0003886 size | .0065137 .0025792 2.53 0.012 .0014586 .0 leverage | .0 .0 2.48 0.013 .0060363 .0 roa | -. .0 -6.41 0.000 -. -. growth | .0 .0028657 14.19 0.000 .0 .0 toptenrate | .0005991 .0001313 4.56 0.000 .0003418 .0008563 dual | .004107 .0047744 0.86 0.390 -.0052507 .0 independent | -. .0 -3.21 0.001 -. -.0 magpay | . . 2.13 0.033 .0 1. boardshare | -. .030857 -3.74 0.000 -. -.055068 invrec | -.0 .0 -1.45 0.146 -.0 .0070255 analyst | -3.14e-06 .0002186 -0.01 0.989 -.0004316 .0004253 opin | -.0 .0083493 -4.17 0.000 -.051178 -.0 aud | .0026224 .0036027 0.73 0.467 -.0044387 .0096835 | year | 2009 | -.0 .0087203 -1.53 0.125 -.0 .0037278 2010 | -.0023983 .0085558 -0.28 0.779 -.0 .0 2011 | .007618 .0084863 0.90 0.369 -.0090149 .0 2012 | -.018136 .0082114 -2.21 0.027 -.03423 -.002042 2013 | -.0 .0082714 -4.15 0.000 -.0 -.0 2014 | -.0 .0082289 -3.26 0.001 -.0 -.0 2015 | .0043649 .0082808 0.53 0.598 -.0 .020595 2016 | -.0049588 .0082864 -0.60 0.550 -.0212 .0 2017 | -.0 .0084003 -2.61 0.009 -.0 -.0054486 2018 | -.0 .008332 -3.59 0.000 -.0 -.0 2019 | -.028234 .0081615 -3.46 0.001 -.0 -.0 | industry | 2 | -.0 .0 -3.40 0.001 -. -.0 3 | -.0 .0 -1.35 0.177 -.0 .012629 4 | -.0 .0 -0.42 0.675 -.0 .0 5 | -.0 .0 -3.41 0.001 -. -.027294 6 | -.0 .0 -2.38 0.017 -.0 -.0081345 | _cons | -. .077478 -1.37 0.172 -. .0 ---------------+---------------------------------------------------------------- coz2 | Laginstitution | .000229 .0017509 0.13 0.896 -.0032026 .0036607 Lagsize | . .0 9.43 0.000 .0 . Lagleverage | . .0 4.54 0.000 . . Lagroa | -. . -1.59 0.111 -. .0 Laggrowth | -.0 .0 -2.49 0.013 -.0 -.0 Lagtoptenrate | .0049017 .0008002 6.13 0.000 .0033334 .0064701 Lagdual | -. .030836 -3.73 0.000 -. -.0 Lagindependent | -1. . -6.37 0.000 -1.9228 -1.017622 Lagmagpay | 1.025986 2. 0.47 0.641 -3.2887 5. Lagboardshare | -1. . -15.56 0.000 -1. -1. Laginvrec | -. .0 -8.31 0.000 -. -. Laganalyst | -.0040897 .0015752 -2.60 0.009 -.0071769 -.0010024 Lagopin | .0 .0 1.33 0.184 -.0 . Lagaud | .0 .0 1.98 0.048 .0005237 .0 _cons | -3. . -12.48 0.000 -4.058582 -2. ---------------+---------------------------------------------------------------- /mills | lambda | .0 .0 5.46 0.000 .0 . ---------------+---------------------------------------------------------------- rho | 0.84225 sigma | . -------------------------------------------------------------------------------- */ *- （3）处理效应模型 - 手工处理（TwoStep） qui: probit coz1 Lagcoz1 $CV i.year i.industry est store First2 predict y2, xb gen imr2 = normalden(y2) / normal(y2) replace imr2 = -normalden(y2) / normal(-y2) if coz2 == 0 qui: reg da1 coz2 $CV i.year i.industry imr2, r est store Second2 *- 结果导出 local mlist2 "OLS2 he_2S2 First2 Second2" esttab `mlist2' using coz2的基准回归及稳健性检验.rtf, replace /// b(%6.4f) t(%6.4f) /// scalar(N r2_a r2_p) compress nogap /// star(* 0.1 0.05 * 0.01) /// indicate("Year FE=*.year" "Industry FE=*.industry") /// mtitle(`mlist2') /// title("Baseline Regression and Robustness Test of coz2") /* Baseline Regression and Robustness Test of coz2 -------------------------------------------------------------- (1) (2) (3) (4) OLS2 he_2S2 First2 Second2 -------------------------------------------------------------- main coz2 -0.0090* -0.0208* -0.0109* (-3.5991) (-1.8260) (-3.8498) institut~n 0.0002 -0.0001 0.0036 0.0002 (1.5711) (-0.4321) (1.2429) (1.5472) size -0.0066* 0.0065 0.1121* -0.0066* (-7.9283) (2.5255) (4.6877) (-7.4912) leverage 0.0213* 0.0287 -0.0209 0.0237* (4.6120) (2.4817) (-0.1712) (4.9087) roa -0.2110* -0.1696* -0.5579* -0.2253* (-12.8242) (-6.4084) (-1.7340) (-13.2397) growth 0.0344* 0.0407* 0.0879* 0.0353* (15.5849) (14.1866) (2.8924) (15.1052) toptenrate 0.0002* 0.0006* 0.0032 0.0002* (5.4342) (4.5641) (2.4032) (4.4256) dual 0.0027* 0.0041 -0.0798 0.0031 (1.9266) (0.8602) (-1.5712) (2.0504) independ~t 0.0034 -0.1094* -1.0194* 0.0040 (0.3042) (-3.2071) (-2.7177) (0.3470) magpay 0.3763* 0.7141 2.6841 0.4066* (2.9216) (2.1271) (0.7287) (2.8833) boardshare -0.0026 -0.1155* -1.0187* -0.0092 (-0.7150) (-3.7446) (-5.8660) (-2.4146) invrec 0.0229* -0.0201 -0.1787 0.0182* (4.6381) (-1.4522) (-1.3510) (3.5220) analyst 0.0006* -0.0000 -0.0034 0.0005* (7.0944) (-0.0144) (-1.3109) (6.6466) opin -0.0408* -0.0348* 0.1167 -0.0389* (-10.2360) (-4.1696) (1.1722) (-9.4915) aud -0.0043* 0.0026 0.0261 -0.0039* (-3.5238) (0.7279) (0.6436) (-3.0258) Lagcoz1 3.2352* (78.3295) imr2 0.0031* (1.6981) _cons 0.2673* -0.1058 -4.3443* 0.2773* (13.3690) (-1.3655) (-7.9174) (12.9075) -------------------------------------------------------------- coz2 Laginsti~n 0.0002 (0.1308) Lagsize 0.1220* (9.4334) Laglever~e 0.3218* (4.5355) Lagroa -0.3268 (-1.5947) Laggrowth -0.0534 (-2.4929) Lagtopte~e 0.0049* (6.1256) Lagdual -0.1151* (-3.7329) Lagindep~t -1.4702* (-6.3668) Lagmagpay 1.0260 (0.4661) Lagboard~e -1.5829* (-15.5642) Laginvrec -0.6019* (-8.3086) Laganalyst -0.0041* (-2.5964) Lagopin 0.0822 (1.3277) Lagaud 0.0487 (1.9813) _cons -3.5077* (-12.4803) -------------------------------------------------------------- /mills lambda 0.0948* (5.4618) Year FE Yes Yes Yes Yes Industry~E Yes Yes Yes Yes -------------------------------------------------------------- N 22591 20018 20018 20018 r2_a 0.1388 0.1442 r2_p 0.6912 -------------------------------------------------------------- t statistics in parentheses * p<0.1, p<0.05, * p<0.01 */ 

$coz3$的基准回归及稳健性检验的第三段代码见公众号下载附件（公众号后台回复heckman即可获取下载链接）。