-
2021-01-15 17:48:44
考虑一类混合效应模型yij=ai+xiβ+eij,j=1,2,…,mi,i=1,2,…,n,其中Eai=0,Eai2=σa2,Eeij=0,Eeij2=σe2.给出了参数β,σa2和σe2的估计,并证明了这些估计量的强相合性.还讨论了随机效应ai的方差为零的检验问题,给出了检验方案.
随机效应模型 random effects
models随机效应模型(random effects models)是经典的线性模型的一种推广，就是把原来(固定)的回归系数看作是随机变量，一般都是假设是来自正态分布。如果模型里一部分系数是随机的，另外一些是固定的，一般就叫做混合模型(mixed models)。
虽然定义很简单，对线性混合模型的研究与应用也已经比较成熟了，但是如果从不同的侧面来看，可以把很多的统计思想方法综合联系起来。概括地来说，这个模型是频率派和贝叶斯模型的结合，是经典的参数统计到高维数据分析的先驱，是拟合具有一定相关结构的观测的典型工具。
随机效应最直观的用处就是把固定效应推广到随机效应。注意，这时随机效应是一个群体概念，代表了一个分布的信息
or 特征，而对固定效应而言，我们所做的推断仅限于那几个固定的(未知的)参数。例如，如果要研究一些水稻的品种是否与产量有影响，如果用于分析的品种是从一个很大的品种集合里随机选取的，那么这时用随机效应模型分析就可以推断所有品种构成的整体的一些信息。这里，就体现了经典的频率派的思想-任何样本都来源于一个无限的群体(population)。
同时，引入随机效应就可以使个体观测之间就有一定的相关性，所以就可以用来拟合非独立观测的数据。经典的就有重复观测的数据，多时间点的记录等等，很多时候就叫做纵向数据(
更多相关内容 -
stata︱政策处理效应模型sata基本命令汇总
2020-12-29 01:06:39本文来源经管之家论坛，由坛友cuifengbao归纳Use ,文件名.dta,clearSsc installpamatch2，replace一、首先做一元回归reg 结果变量 处理变量，r二、直接引入协变量，再做多元回归reg 结果变量 处理变量 协变量1 协...展开全文本文来源经管之家论坛，由坛友cuifengbao归纳
Use ,文件名.dta,clear
Ssc installpamatch2，replace
一、首先做一元回归
reg 结果变量 处理变量，r
二、直接引入协变量，再做多元回归
reg 结果变量 处理变量 协变量1 协变量2 协变量3……，r
三、接下来进行倾向得分匹配
1.将数据随机排序
set seed 10101
gen ranorder = runiform()
sort ranorder
2.一对一匹配，又放回匹配，允许并列
Psmatch2 处理变量 协变量1 协变量2 …,outcome(结果变量) n(1) ate ties logit common
解释：n()表示一对几匹配
ate 表示同时汇报ATE ATU ATT。如果默认表示仅汇报ATT
ties 表示并列个体
logit表示对数单位模型
common表示仅对共同取值范围内个体进行匹配。如果默认表示对所有个体进匹配。
3.使用引导程序显示全部处理效应结果
Set seed 10101
Bootstrap r(att) r(atu) r(ate),reps(500):psmatch2处理变量 协变量1 协变量2 …,outcome(结果变量) n(1) ate ties logit common
4.使用命令pstest考察匹配结果是否较好地平衡了数据
Quietly psmatch2 处理变量 协变量1 协变量2 …，outcome(结果变量) n(1) ate ties Logit co
-
『Stata』政策处理效应PSM模型基本命令汇总
2020-12-19 08:55:08use ,文件名.dta,clearssc install pamatch2，replace一、首先做一元回归reg 结果变量 处理变量，r二、直接引入协变量，再做多元回归reg 结果变量 处理变量 协变量1 协变量2 协变量3……，r三、接下来进行倾向得分...展开全文use ,文件名.dta,clear
ssc install pamatch2，replace
一、首先做一元回归
reg 结果变量 处理变量，r
二、直接引入协变量，再做多元回归
reg 结果变量 处理变量 协变量1 协变量2 协变量3……，r
三、接下来进行倾向得分匹配
1.将数据随机排序
set seed 10101
gen ranorder = runiform()
sort ranorder
2.一对一匹配，又放回匹配，允许并列
Psmatch2 处理变量 协变量1 协变量2 …,outcome(结果变量) n(1) ate ties logit common
解释：n()表示一对几匹配
ate 表示同时汇报ATE ATU ATT。如果默认表示仅汇报ATT
ties 表示并列个体
logit表示对数单位模型
common表示仅对共同取值范围内个体进行匹配。如果默认表示对所有个体进匹配。
3.使用引导程序显示全部处理效应结果
Set seed 10101
Bootstrap r(att) r(atu) r(ate),reps(500):psmatch2处理变量 协变量1 协变量2 …,outcome(结果变量) n(1) ate ties logit common
4.使用命令pstest考察匹配结果是否较好地平衡了数据
Quietly psmatch2 处理变量 协变量1 协变量2 …，outcome(结果变量) n(1) ate ties Logit common
pstest协变量1 协变量2 …，both graph
5.用条形图显示倾向得分的共同取值范围
psgraph
6.K进行近邻匹配，令k=4，为节省空间，可使用选择项“quietly”略去对倾向得分估计结果的汇报。
psmatch2处理变量 协变量1 协变量2 …，outcome(结果变量) n(4) ate ties Logit common quietly
7.进行卡尺内一对四匹配。计算倾向得分的标准差，然后乘0.25
sum _psscore
dis 0.25*r(sd)
01979237
psmatch2处理变量 协变量1 协变量2 …，outcome(结果变量) n(4)cal(0.01) ate ties Logit common quietly
8.进行半径(卡尺)匹配
psmatch2处理变量 协变量1 协变量2 …，outcome(结果变量) radius cal(0.01) ate ties Logitcommon quietly
9.进行核匹配(使用默认的核函数与带宽)
psmatch2处理变量 协变量1 协变量2 …，outcome(结果变量) kernel ate ties Logit common quietly
10.进行局部线性回归匹配(使用默认的核函数与带宽)
psmatch2处理变量 协变量1 协变量2 …，outcome(结果变量) llr ate ties Logit common quietly
然后可以使用自助法得到自助标准差
Set seed 10101
Bootstrap r(att) r(atu) r(ate),reps(500):psmatch2处理变量 协变量1 协变量2 …,outcome(结果变量) llr ate ties logit common
11.进行样条匹配和马氏匹配省略不介绍。
-
Heckman两步法 | 样本选择模型 & 处理效应模型
2021-08-07 19:24:13这期推送简单介绍一下样本选择模型和处理效应模型，其中样本选择模型是一般意义上的Heckman两步法，后者则借鉴了Heckman两步法的构建思想，但又不完全等同于前者。模型介绍之后，将利用help文件中的示例数据与代码...展开全文这期推送简单介绍一下样本选择模型和处理效应模型，其中样本选择模型是一般意义上的
Heckman两步法，后者则借鉴了Heckman两步法的构建思想，但又不完全等同于前者。模型介绍之后，将利用help文件中的示例数据与代码简单演示一下这两个模型在Stata中的具体操作，然后简单评述一下现阶段文献中对这两个模型的理解与应用情况，最后结合一篇论文的公开数据与代码进行结果复盘与二次验证。N o t e : Note: Note: 1、下划线字体为链接，可点击跳转；2、推文中的公式与代码块均可左右滑动；3、该文首发于微信公众号
DMETP，欢迎关注；4、需要本次推送所使用的数据和代码的朋友，可以在公众号后台对话框内回复关键词heckman。文章目录
一、样本选择偏差与自选择偏差
上期推送『双重差分法 | PSM - DID』介绍了样本选择偏差与自选择偏差的区别，最关键的一点在于两者非随机的选择机制是不同的。
-
样本选择偏差。样本选择偏差的非随机选择机制在于对样本的选择不随机。在样本数据的采集过程中，只对某部分群体进行调查，但这部分群体与其他群体在某些方面的特征差异较大，因此根据这样的样本做回归得到的普适性结论并不可信。体现在具体的数据集中就是，数据集中只有特定群体的样本，或者，虽然有全部群体的所有解释变量数据，但除特定群体之外的其他群体的被解释变量数据缺失，在这两种情况下进行的回归，都将直接忽视其他群体的样本信息（
y缺失的样本在参与回归时将被drop掉）。实质上，样本选择偏差说的就是参与回归的样本不能代表总体从而产生估计偏误的问题。 -
自选择偏差。自选择偏差的非随机选择机制在于对自变量的选择不随机。在使用DID方法评估政策效应时，一个明显的事实就是，相对于未实施政策的地区（控制组），实施政策的地区（处理组）通常情况下经济发展都较为发达、各类基础设施建设都较为完善，而所谓的“政策效果评估”也即考察政策的经济效应，因此地区是否参与政策这一行为是内生的。体现在回归方程中就是，经济指标（如，
GDP、人均GDP、GDP增长率等）作为被解释变量y，地区（在某时点）是否实施该项政策的哑元变量D作为核心解释变量，但由于政策内生，因此某些影响地区是否参与决策D的（可观测或不可观测）因素也将同时影响经济指标y，由于这些因素或者无法穷尽、或者影响形式未知、或者不可测度，因此被放到随机扰动项中，造成解释变量D与扰动项ε相关，即 C o v ( D , ε ) ≠ 0 Cov(D, \varepsilon)\neq 0 Cov(D,ε)=0。实质上，自选择偏差说的就是实验组与控制组的先验条件存在较大差异从而导致估计偏误的问题。 -
两者的区别。非随机选择机制的不同是两者最大的区别，体现在具体回归方程中就是，样本选择偏差中被解释变量
y是否被观测到或是否取值（而非取值大小）是非随机的；而自选择偏差中哑元解释变量D的取值是非随机的。陈强（2014）《高级计量经济学及Stata应用（第二版）》第539页认为，样本选择问题通常不考虑某项目或政策的效应，故个体间的差异并不在于是否得到处理，而在于是否能进入样本（即被解释变量 y i y_i yi是否可观测），通常 D i = 1 D_i=1 Di=1意味着 y i y_i yi可观测，而 D i = 0 D_i=0 Di=0则意味着 y i y_i yi不可观测。而在处理效应模型中，无论 D i = 1 D_i=1 Di=1或 0 0 0，结果变量 y i y_i yi均可观测。这种说法基本概括了两者的区别，但有一个小问题，在样本选择偏差中， D i D_i Di的取值与 y i y_i yi是否可观测并不存在必然的关系，因为 D i D_i Di是一个确定并可准确测度的因素，而影响 y i y_i yi是否可观测的却是一个不可观测的潜变量，这个潜变量由一系列控制变量与外生变量决定。
二、两个模型的估计思路
花大篇幅论述样本选择偏差与自选择偏差这两个问题，自然是为了引出解决这两个问题的具体方法。
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），计算公式为：I M R i ( 或 λ i ) = ϕ ( y ^ i ) Φ ( y ^ i ) (1) IMR_i(或\lambda_i)=\frac{\phi(\hat y_i)}{\Phi(\hat y_i)} \tag{1} IMRi(或λi)=Φ(y^i)ϕ(y^i)(1)
其中， y ^ i \hat y_i y^i为第 i i i个样本在第一步回归（选择方程）的拟合值， ϕ ( ⋅ ) \phi(·) ϕ(⋅)为标准正态的概率密度函数（Probability Density Function，
pdf）， Φ ( ⋅ ) \Phi(·) Φ(⋅)为累积分布函数（Cumulative Distribution Function，cdf）。因此，样本选择模型的估计思路是：首先，计算全部样本的
IMR；随后，将遗漏变量IMR代入原回归方程中，具体来说：-
第一步 ：用probit方法估计选择方程，其中原回归方程的被解释变量
y是否被观测到或是否取值的虚拟变量y_dummy作为probit的被解释变量，解释变量包括原回归方程所有解释变量和至少一个外生变量，该外生变量只影响y是否取值，而不影响y的大小，即满足相关性和外生性的要求（但不是工具变量）。估计出所有变量的系数后，将样本数据代入至probit模型中，计算出拟合值 y ^ \hat y y^，再将 y ^ \hat y y^代入风险函数（公式 ( 1 ) (1) (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计算公式为：I M R i ( 或 λ i ) = { ϕ ( y ^ i ) Φ ( y ^ i ) , i f D i = 1 − ϕ ( y ^ i ) Φ ( − y ^ i ) , i f D i = 0 (2) IMR_i(或\lambda_i)= \begin{cases} \frac{\phi(\hat y_i)}{\Phi(\hat y_i)}&,&if~D_i=1\\ \frac{-\phi(\hat y_i)}{\Phi(-\hat y_i)}&,&if~D_i=0\\ \end{cases}\tag{2} IMRi(或λi)={Φ(y^i)ϕ(y^i)Φ(−y^i)−ϕ(y^i),,if Di=1if Di=0(2)
上式各字母的解释同公式 ( 1 ) (1) (1)。关于
IMR计算公式的更多细节，请参考Stata官方网站的回答（FAQs）。明显可以看到，公式 ( 1 ) (1) (1)说明在样本选择模型中，所有样本的
IMR均用一个公式来计算；公式 ( 2 ) (2) (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的例子。
- 需要注意的是，选择方程中的工具变量应尽量避免使用
-
第二步 ：将样本数据代入第一步选择方程中，得到各个样本的的拟合值 y ^ i \hat y_i y^i，再将 y ^ i \hat y_i y^i代入处理效应模型的风险函数（公式 ( 2 ) (2) (2)）中，计算得到各样本的
IMR，最后将IMR作为额外的控制变量引入原回归方程中，考察核心解释变量D以及IMR的估计系数。如果IMR的估计系数显著，说明自选择偏差问题不可忽视，此时核心解释变量D的系数就是考虑了自选择偏差后的估计结果，并可与基准回归结果对比构成稳健性检验；而如果IMR的估计系数不显著，则说明自选择偏差问题在原回归中不明显，基准回归结果本身就是可信的。- 需要注意的是，核心解释变量
D在两步模型中均参与了回归，其中第一阶段回归中D作为被解释变量，在第二阶段回归中作为解释变量，并且我们假设D不存在缺失值，因此处理效应模型两步回归中的样本均是全样本，这不同于样本选择模型。
- 需要注意的是，核心解释变量
2.3 估计思路的对比
总结一下样本选择模型和处理效应模型的估计思路的异同点。
相同点在于：
-
都是两步估计法。Heckman于1979年提出的两步估计法最开始是用于解决样本选择偏差的，即最初的Heckman两步法指的就是样本选择模型，后来有学者借鉴这种两步估计法的思想，应用于解决自选择偏差的处理效应模型。这两个模型在估计思路上是一脉相承的，而正是因为这种相似性，所以才导致各个学者对这两个模型的错误理解与错误应用，这种错误在现阶段的文献中较为常见。
-
都可以使用MLE进行模型的整体估计。两步估计法（如2SLS、PSM - DID以及这里的样本选择模型和处理效应模型等）一个明显的缺陷是，第一步估计的误差将被带入第二步，导致效率损失。而使用MLE从整体上进行参数估计可以避免这种问题，但如果样本量过大，MLE估计耗时较长，且MLE对分布的假设较为严格，因此需要在估计的精准性、操作的简便性等方面进行权衡。
-
第一阶段回归都需要引入外生变量，同时应包括第二阶段的所有外生解释变量。引入的外生变量需满足相关性和外生性的要求，即与选择方程中的被解释变量在理论上和统计上均具有相关性，而与第二步回归的被解释变量不具有直接的相关关系。引入外生变量的目的是确保第一步计算得到的
IMR在引入原回归方程后不与干扰项相关。该外生变量在处理效应模型中可以直接称作工具变量。此外，如果核心解释变量D是DID模型的did项，那么为了防止出现多重共线性，应该尽量避免使用D的滞后项D_lag作为工具变量。事实上，如果找到了一个良好的工具变量，也完全能够使用2SLS解决内生性问题。此外，两个模型除了都需要在第一阶段引入至少一个外生变量，第一阶段回归中的其余控制变量也应该是第二阶段回归中所有的控制变量，即应该包括所有的外生解释变量，原因在于保证两阶段估计的一致性，详情请看陈强教授的推文『工具变量法（五）: 为何第一阶段回归应包括所有外生解释变量』。然而，部分文献在第一阶段并未包括第二阶段所有的外生解释变量，少部分文献甚至根本就不引入第二阶段的外生解释变量（如，考虑滞后效应，直接引入第二阶段外生解释变量的滞后项），并且在Stata处理效应模型的官方命令etregress的help文件的演示案例中，第一阶段回归也并未包括所有的外生解释变量，原因可能在于IMR是一个非线性项，因此不包含所有外生解释变量引起的内生性问题可能并没有2SLS那么严重。 -
第一步回归都只能是probit模型。由于logit模型不具备扰动项服从正态分布的假设，如果使用logit模型估计选择方程，将直接导致
IMR计算错误，因为Heckman（1979）在推导IMR时，假设选择方程的随机扰动项服从正态分布。这与PSM不同，PSM估计概率方程可以使用logit模型，也可以使用probit模型，并且实际使用中流行的是logit模型。然而，选择方程使用probit模型进行估计有一个问题不可忽视，那就是probit（包括Stata的xtprobit）不能估计固定效应模型，因此即便在回归方程中引入时间虚拟变量和个体虚拟变量，控制的也只是“时间效应”和“个体效应”，不能加入“固定”二字。
不同点在于：
-
解决的问题不同。样本选择模型解决的是样本选择偏差导致的内生性问题，处理效应模型解决（或者“缓解”）的是依不可观测因素导致的自选择偏差问题。在实际应用中，部分文献在分析内生性问题时将样本选择偏差与自选择偏差混淆，从而使用的模型也是不恰当的。在数据搜集过程中，对被解释变量存在缺失值的样本，多数文献的做法是直接把这些样本剔除，因而即便文章中考虑到了样本选择偏差问题，我们也无法使用样本选择模型（或Heckman两步法）。事实上，囿于数据缺陷，大多数实证类论文都不具备实施Heckman两步法的条件。对于DID类的实证论文，对内生性的分析角度应该更多考虑从自选择偏差切入，而非样本选择偏差，因为各样本处理组虚拟变量
D的取值本身就提供了自选择偏差分析的条件，即D取值为1的样本与D取值为0的样本在某些方面是否存在明显的特征差异？或者，是否存在某些因素影响了各样本是否实施政策的决定，而这些因素在两组间又是否存在巨大差异？同时，这些因素是否在理论与统计意义上影响我们想研究的经济指标？在这样的分析之后，就可以使用处理效应模型来缓解因自选择偏差而导致的估计偏误。 -
变量的设置不同。在样本选择模型第一阶段回归方程中，被解释变量是原方程中的被解释变量
y是否被观测到的虚拟变量y_dummy，该变量不参与第二阶段回归，同时第一阶段引入的外生变量直接影响的是y_dummy。在处理效应模型第一阶段回归方程中，被解释变量是原方程的核心解释变量D，D取值为0或1，且不存在缺失值，该变量还同时参与了第二阶段回归，此外第一阶段引入的外生变量（或称工具变量）直接影响的是D。 -
各阶段样本参与回归的数目不同。假设除关键变量，其余变量都不存在缺失值，那么对于样本选择模型来说，第一阶段回归的解释变量均不存在缺失值，被解释变量
y_dummy取值为0或1，也不存在缺失值，因此选择方程中参与回归的样本是全样本，第二阶段由于被解释变量y本身就存在缺失值，因此参与第二阶段回归的样本不是全样本，从而第一阶段的样本多于第二阶段。对于处理效应模型来说，所有变量均不存在缺失值，因此两阶段参与回归的样本是相同的，虽然在第一阶段引入滞后项D_lag作为工具变量的情况下会损失一部分样本，但由于计算出来的IMR同样也存在缺失值，从而第二阶段参与回归的样本也将与第一阶段相同。 -
IMR的计算公式不同。从公式 ( 1 ) (1) (1)和公式 ( 2 ) (2) (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.7245231 R-squared = 0.2535 -------------+---------------------------------- Adj R-squared = 0.2524 Total | 53354.8946 1,342 39.7577456 Root MSE = 5.452 ------------------------------------------------------------------------------ wage | Coef. Std. Err. t P>|t| [95% Conf. Interval] -------------+---------------------------------------------------------------- education | .8965829 .0498061 18.00 0.000 .7988765 .9942893 age | .1465739 .0187135 7.83 0.000 .109863 .1832848 _cons | 6.084875 .8896182 6.84 0.000 4.339679 7.830071 ------------------------------------------------------------------------------ */第二个是样本选择模型，使用MLE方法进行估计，可以看到：
- 选择方程中两个外生变量均显著为正，说明外生变量的选择是有效的。
- 在第二阶段回归中，
IMR（即lambda）的估计系数为4.2244，但显著性未知，该值等于rho和sigma的乘积，其中：sigma是原方程干扰项的标准差。rho是选择方程干扰项和第二阶段回归干扰项的相关系数。如果rho等于0，表示第二阶段回归中IMR的系数不显著，说明样本选择偏差在原方程中不怎么严重，反之则需要考虑样本选择偏差带来的估计偏误。回归结果的末尾是LR检验，检验的原假设是H0: rho = 0，p值说明至少可以在1%的水平下拒绝原假设，可以认为rho显著不等于0，这说明原模型中确实存在严重的样本选择偏差，基准回归结果不可信。
- 第二阶段回归结果中，两个解释变量仍旧显著为正，且相较于基准回归结果取值变化不大，说明考虑到样本选择偏差后基准回归结果依然是稳健的。
- 回归结果右上角说明总的样本量是2,000个，参与选择方程回归的样本为2,000个，参与第二阶段回归的样本为1,343个。
*-（二）MLE heckman wage $xlist , select( wage_dummy = $xlist $zlist ) nolog mills(imr1) est store MLE /* Heckman selection model Number of obs = 2,000 (regression model with sample selection) Selected = 1,343 Nonselected = 657 Wald chi2(2) = 508.44 Log likelihood = -5178.304 Prob > chi2 = 0.0000 ------------------------------------------------------------------------------ | Coef. Std. Err. z P>|z| [95% Conf. Interval] -------------+---------------------------------------------------------------- wage | education | .9899537 .0532565 18.59 0.000 .8855729 1.094334 age | .2131294 .0206031 10.34 0.000 .1727481 .2535108 _cons | .4857752 1.077037 0.45 0.652 -1.625179 2.59673 -------------+---------------------------------------------------------------- wage_dummy | education | .0557318 .0107349 5.19 0.000 .0346917 .0767718 age | .0365098 .0041533 8.79 0.000 .0283694 .0446502 married | .4451721 .0673954 6.61 0.000 .3130794 .5772647 children | .4387068 .0277828 15.79 0.000 .3842534 .4931601 _cons | -2.491015 .1893402 -13.16 0.000 -2.862115 -2.119915 -------------+---------------------------------------------------------------- /athrho | .8742086 .1014225 8.62 0.000 .6754241 1.072993 /lnsigma | 1.792559 .027598 64.95 0.000 1.738468 1.84665 -------------+---------------------------------------------------------------- rho | .7035061 .0512264 .5885365 .7905862 sigma | 6.004797 .1657202 5.68862 6.338548 lambda | 4.224412 .3992265 3.441942 5.006881 ------------------------------------------------------------------------------ LR test of indep. eqns. (rho = 0): chi2(1) = 61.20 Prob > chi2 = 0.0000 */第三个同样是样本选择模型，使用MLE方法进行估计，这里将使用异方差稳健标准误。可以看到，使用稳健标准误的回归结果与上文使用普通标准误的回归结果基本保持一致。需要注意的是，在使用稳健标准误后，对
H0的统计检验将变成Wald检验，这里Wald检验同样拒绝原假设。*-（三）MLE + robust SE heckman wage $xlist , select( wage_dummy = $xlist $zlist ) nolog vce(robust) est store MLE_r /* Heckman selection model Number of obs = 2,000 (regression model with sample selection) Selected = 1,343 Nonselected = 657 Wald chi2(2) = 497.82 Log pseudolikelihood = -5178.304 Prob > chi2 = 0.0000 ------------------------------------------------------------------------------ | Robust | Coef. Std. Err. z P>|z| [95% Conf. Interval] -------------+---------------------------------------------------------------- wage | education | .9899537 .0534141 18.53 0.000 .8852641 1.094643 age | .2131294 .0211095 10.10 0.000 .1717555 .2545034 _cons | .4857752 1.099121 0.44 0.659 -1.668462 2.640013 -------------+---------------------------------------------------------------- wage_dummy | education | .0557318 .0108899 5.12 0.000 .034388 .0770755 age | .0365098 .0042243 8.64 0.000 .0282303 .0447893 married | .4451721 .0668243 6.66 0.000 .3141988 .5761453 children | .4387068 .0272779 16.08 0.000 .385243 .4921705 _cons | -2.491015 .1884227 -13.22 0.000 -2.860317 -2.121713 -------------+---------------------------------------------------------------- /athrho | .8742086 .1051331 8.32 0.000 .6681514 1.080266 /lnsigma | 1.792559 .0288316 62.17 0.000 1.73605 1.849068 -------------+---------------------------------------------------------------- rho | .7035061 .0531006 .5837626 .7932976 sigma | 6.004797 .1731281 5.674882 6.353893 lambda | 4.224412 .4172197 3.406676 5.042147 ------------------------------------------------------------------------------ Wald test of indep. eqns. (rho = 0): chi2(1) = 69.14 Prob > chi2 = 0.0000 */第四个是样本选择模型，使用两步法进行估计，可以看到：
- 选择方程中两个工具变量的引入是有效的。
- 第二阶段回归中，
IMR的回归系数等于4.0016，与MLE方法下的4.2244相差不大，但两步法下IMR回归系数可以直接进行z检验，并且统计结果说明IMR回归系数至少在1%的水平下显著为正，这同时说明原方程中的样本选择偏差问题不可忽视。 - 第二阶段回归结果中，两个解释变量仍旧显著为正，且大小与基准回归结果相比变化不大，这说明在考虑样本选择偏差的情况下，基准回归结果是可信的。
- 右上角同样说明在两步估计法下，参与选择方程的样本数是2,000个，但参与第二阶段回归的样本数只有1,343个，减少的657个样本与被解释变量
wage存在缺失值的样本完全重合。
*-（四）两步估计法 heckman wage $xlist , select( wage_dummy = $xlist $zlist ) twostep mills(imr2) est store TwoStep /* Heckman selection model -- two-step estimates Number of obs = 2,000 (regression model with sample selection) Selected = 1,343 Nonselected = 657 Wald chi2(2) = 442.54 Prob > chi2 = 0.0000 ------------------------------------------------------------------------------ | Coef. Std. Err. z P>|z| [95% Conf. Interval] -------------+---------------------------------------------------------------- wage | education | .9825259 .0538821 18.23 0.000 .8769189 1.088133 age | .2118695 .0220511 9.61 0.000 .1686502 .2550888 _cons | .7340391 1.248331 0.59 0.557 -1.712645 3.180723 -------------+---------------------------------------------------------------- wage_dummy | education | .0583645 .0109742 5.32 0.000 .0368555 .0798735 age | .0347211 .0042293 8.21 0.000 .0264318 .0430105 married | .4308575 .074208 5.81 0.000 .2854125 .5763025 children | .4473249 .0287417 15.56 0.000 .3909922 .5036576 _cons | -2.467365 .1925635 -12.81 0.000 -2.844782 -2.089948 -------------+---------------------------------------------------------------- /mills | lambda | 4.001615 .6065388 6.60 0.000 2.812821 5.19041 -------------+---------------------------------------------------------------- rho | 0.67284 sigma | 5.9473529 ------------------------------------------------------------------------------ */第五个是手工完成两步估计法，但实际使用时这种方法十分不推荐（包括上文使用两步法估计的样本选择模型），因为两步估计法下第一步回归的估计误差会带到第二步，造成效率损失，影响第二步估计系数的显著性。
手工两步估计法的思路是先用控制变量和外生变量对
wage_dummy做probit回归，计算出拟合值y_hat，然后根据y_hat计算出全部样本的IMR，最后将IMR作为额外的控制变量引入基准回归方程中做进一步的回归。第一阶段（选择方程）回归的代码与结果如下。
*-（五）手工两步估计法 *- 第一阶段方程（选择方程）回归 probit wage_dummy $xlist $zlist , nolog est store First /* Probit regression Number of obs = 2,000 LR chi2(4) = 478.32 Prob > chi2 = 0.0000 Log likelihood = -1027.0616 Pseudo R2 = 0.1889 ------------------------------------------------------------------------------ wage_dummy | Coef. Std. Err. z P>|z| [95% Conf. Interval] -------------+---------------------------------------------------------------- education | .0583645 .0109742 5.32 0.000 .0368555 .0798735 age | .0347211 .0042293 8.21 0.000 .0264318 .0430105 married | .4308575 .074208 5.81 0.000 .2854125 .5763025 children | .4473249 .0287417 15.56 0.000 .3909922 .5036576 _cons | -2.467365 .1925635 -12.81 0.000 -2.844782 -2.089948 ------------------------------------------------------------------------------ */可以看到，实际参与回归的样本为2,000个，控制变量和两个外生变量均显著为正，这说明外生变量的选择是有效的。此外，衡量模型整体拟合程度的统计指标是伪R方（
Pseudo R2），伪R方为0.1889，整体拟合程度较好。下面根据拟合值
y_hat计算出额外的控制变量imr。其中，normalden()是标准正态的概率密度函数pdf，normal()是标准正态的累积分布函数cdf。*- 计算IMR predict y_hat, xb gen imr = normalden(y_hat) / normal(y_hat)获得各个样本的
imr后，就可以进行第二阶段回归。可以发现：- 与样本选择模型的两步估计法结果相比，手工两步法估计结果在系数值大小方面没有任何改变，在系数标准误方面变化也不大，从而各个变量的系数显著性保持高度一致。
- 这提示我们，两步估计法（包括手工两步法）带来的效率损失不可忽视，虽然本例中的结果和MLE以及基准回归基本保持一致，实际使用应尽量避免使用这种方法，除非能够找到比较好的方法来校正误差。
- 第二阶段实际参与回归的样本是1,343个。
- 衡量模型整体拟合程度的统计指标是修正R方，修正R方等于0.2777，整体拟合程度较高。事实上，由于该例并未规定原方程的核心解释变量，因此对模型整体拟合程度的考察是有意义的，但在现实应用中，R方即便很小问题也不大，只要求核心解释变量的系数大小和显著性符合预期即可。就算对R方有具体的要求，只要我们在方程中引入足够多的控制变量，控制足够多的“固定效应”，R方总能达到预定的要求。因此，实际应用中R方意义不大。
*- 第二阶段方程回归 reg wage $xlist imr est store Second /* Source | SS df MS Number of obs = 1,343 -------------+---------------------------------- F(3, 1339) = 173.01 Model | 14904.6806 3 4968.22688 Prob > F = 0.0000 Residual | 38450.214 1,339 28.7156191 R-squared = 0.2793 -------------+---------------------------------- Adj R-squared = 0.2777 Total | 53354.8946 1,342 39.7577456 Root MSE = 5.3587 ------------------------------------------------------------------------------ wage | Coef. Std. Err. t P>|t| [95% Conf. Interval] -------------+---------------------------------------------------------------- education | .9825259 .0504982 19.46 0.000 .8834616 1.08159 age | .2118695 .0206636 10.25 0.000 .171333 .252406 imr | 4.001616 .5771027 6.93 0.000 2.869492 5.133739 _cons | .7340391 1.166214 0.63 0.529 -1.553766 3.021844 ------------------------------------------------------------------------------ */可以对比一下用各个方法计算的
IMR。其中，imr1是用样本选择模型的MLE方法计算的，imr2是用样本选择模型的两步估计法计算的，imr是用手工两步法计算的。打开Stata的数据编辑器，可以发现：- 所有样本的
IMR均有取值。 imr和imr2保持高度一致（imr四舍五入了），因为两者都是用两步估计法计算出来的。imr1和imr2存在差距，但差别不大。
*- 对比imr、imr1和imr2 order imr1 imr2 imr wage_dummy最后是五个模型回归结果的对比，可以发现：
- 两个解释变量（
education和age）的回归系数的符号、大小、显著性在模型 ( 2 ) ( 3 ) ( 4 ) ( 6 ) (2)(3)(4)(6) (2)(3)(4)(6)中基本保持一致，并且与模型 ( 1 ) (1) (1)的区别不大，这说明在考虑到样本选择偏差后，根据基本回归结果得出的结论依然成立。 - 模型
(
2
)
(
3
)
(2)(3)
(2)(3)缺少变量
lambda，在实际汇报时应该根据原始回归表进行汇报。 - 模型 ( 2 ) ( 3 ) ( 4 ) (2)(3)(4) (2)(3)(4)汇报的样本量是2,000个全样本，但应该清楚不同回归阶段实际参与的样本数目是不同的。
- 关于模型拟合程度的统计量，模型
(
1
)
(1)
(1)和模型
(
6
)
(6)
(6)汇报的都是修正R方，模型
(
5
)
(5)
(5)的是伪R方，模型
(
2
)
(
3
)
(
4
)
(2)(3)(4)
(2)(3)(4)没有汇报统计量，但观察原始回归表可以发现，模型
(
2
)
(
3
)
(
4
)
(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.582425 Prob > F = 0.0000 Residual | 4196.70492 1,203 3.48853277 R-squared = 0.3402 -------------+---------------------------------- Adj R-squared = 0.3369 Total | 6360.19947 1,209 5.26071089 Root MSE = 1.8678 ------------------------------------------------------------------------------ wage | Coef. Std. Err. t P>|t| [95% Conf. Interval] -------------+---------------------------------------------------------------- union | 1.003913 .1338083 7.50 0.000 .7413894 1.266437 age | .1475526 .0195242 7.56 0.000 .1092474 .1858578 grade | .4368545 .0294718 14.82 0.000 .3790327 .4946764 smsa | .9754817 .1252669 7.79 0.000 .7297159 1.221248 black | -.6183346 .1252 -4.94 0.000 -.8639693 -.3726999 tenure | .2118016 .0338481 6.26 0.000 .1453938 .2782094 _cons | -4.326028 .5315474 -8.14 0.000 -5.368891 -3.283165 ------------------------------------------------------------------------------ */第二个是处理效应模型，使用MLE方法从模型整体角度对参数进行估计。可以看到：
- 第一步回归中
south显著为负，说明工具变量的选择有效。 - 第二步回归中核心解释变量
union显著为正，这一点在处理效应模型回归结果表中显示为1.union，说明相对于union取值为0的样本，union取值为1的样本的平均工资高2.9458个单位。2.9458大约是基准回归结果1.0039的三倍，这说明在未考虑自选择偏差的情况下，低估了union对wage的影响。 - 其余解释变量估计系数的符号和显著性基本没变，但估计系数的数值改变较大，这说明自选择偏差带来的估计偏误不能忽视。
IMR（lambda）的估计系数为-1.1603，同时rho的LR检验结果说明至少在1%的水平下拒绝H0。这都说明基准回归模型中确实存在严重的自选择偏差。- 第一步回归和第二步回归的样本数均为1,210个。
*-（二）MLE etregress wage $xlist2 , treat(union = $xlist1 south) nolog hazard(imr1) est store MLE /* Linear regression with endogenous treatment Number of obs = 1,210 Estimator: maximum likelihood Wald chi2(6) = 681.89 Log likelihood = -3051.575 Prob > chi2 = 0.0000 ------------------------------------------------------------------------------ | Coef. Std. Err. z P>|z| [95% Conf. Interval] -------------+---------------------------------------------------------------- wage | age | .1487409 .0193291 7.70 0.000 .1108566 .1866252 grade | .4205658 .0293577 14.33 0.000 .3630258 .4781058 smsa | .9117044 .1249041 7.30 0.000 .6668969 1.156512 black | -.7882471 .1367078 -5.77 0.000 -1.05619 -.5203048 tenure | .1524015 .0369596 4.12 0.000 .0799621 .2248409 1.union | 2.945815 .2749621 10.71 0.000 2.4069 3.484731 _cons | -4.351572 .5283952 -8.24 0.000 -5.387208 -3.315936 -------------+---------------------------------------------------------------- union | black | .4557499 .0958042 4.76 0.000 .2679771 .6435226 tenure | .0871536 .0232483 3.75 0.000 .0415878 .1327195 south | -.5807419 .0851111 -6.82 0.000 -.7475566 -.4139271 _cons | -.8855758 .0724506 -12.22 0.000 -1.027576 -.7435753 -------------+---------------------------------------------------------------- /athrho | -.6544347 .0910314 -7.19 0.000 -.832853 -.4760164 /lnsigma | .7026769 .0293372 23.95 0.000 .645177 .7601767 -------------+---------------------------------------------------------------- rho | -.5746478 .060971 -.682005 -.4430476 sigma | 2.019151 .0592362 1.906325 2.138654 lambda | -1.1603 .1495097 -1.453334 -.8672668 ------------------------------------------------------------------------------ LR test of indep. eqns. (rho = 0): chi2(1) = 19.84 Prob > chi2 = 0.0000 */第三个是处理效应模型，使用MLE方法并使用异方差稳健标准误。可以看到，与使用普通标准误的MLE方法相比，除系数的z统计量有细微差别，以及对
rho的显著性检验使用的是Wald统计指标，其他结果没有实质性差异。*-（三）MLE + robust SE etregress wage $xlist2 , treat(union = $xlist1 south) nolog vce(robust) est store MLE_r /* Linear regression with endogenous treatment Number of obs = 1,210 Estimator: maximum likelihood Wald chi2(6) = 487.30 Log pseudolikelihood = -3051.575 Prob > chi2 = 0.0000 ------------------------------------------------------------------------------ | Robust | Coef. Std. Err. z P>|z| [95% Conf. Interval] -------------+---------------------------------------------------------------- wage | age | .1487409 .0206909 7.19 0.000 .1081875 .1892943 grade | .4205658 .0377922 11.13 0.000 .3464945 .4946371 smsa | .9117044 .1199377 7.60 0.000 .6766308 1.146778 black | -.7882471 .1408988 -5.59 0.000 -1.064404 -.5120906 tenure | .1524015 .0473782 3.22 0.001 .0595419 .2452611 1.union | 2.945815 .5643841 5.22 0.000 1.839643 4.051988 _cons | -4.351572 .6371071 -6.83 0.000 -5.600279 -3.102865 -------------+---------------------------------------------------------------- union | black | .4557499 .093489 4.87 0.000 .2725148 .6389849 tenure | .0871536 .0251543 3.46 0.001 .0378521 .1364552 south | -.5807419 .0837513 -6.93 0.000 -.7448914 -.4165924 _cons | -.8855758 .0739533 -11.97 0.000 -1.030522 -.7406301 -------------+---------------------------------------------------------------- /athrho | -.6544347 .2322442 -2.82 0.005 -1.109625 -.1992445 /lnsigma | .7026769 .0758152 9.27 0.000 .5540819 .8512719 -------------+---------------------------------------------------------------- rho | -.5746478 .1555525 -.8039298 -.1966492 sigma | 2.019151 .1530822 1.740342 2.342624 lambda | -1.1603 .3823277 -1.909649 -.4109519 ------------------------------------------------------------------------------ Wald test of indep. eqns. (rho = 0): chi2(1) = 7.94 Prob > chi2 = 0.0048 */第四个是用两步法估计处理效应模型，可以发现：
- 选择方程中的工具变量选择有效。
- 第二阶段回归中的所有解释变量，除
tenure，均至少在1%的水平下显著，lambda显著为负。 - 回归系数的大小与MLE方法下的估计系数相近，但与基准回归相差较大，同时考虑到
IMR显著，这说明原模型中的自选择偏差不可忽视，在未考虑自选择偏差的情况下，将严重低估union对wage的影响。
*-（四）两步估计法 etregress wage $xlist2 , treat(union = $xlist1 south) twostep hazard(imr2) est store TwoStep /* Linear regression with endogenous treatment Number of obs = 1210 Estimator: two-step Wald chi2(8) = 566.56 Prob > chi2 = 0.0000 ------------------------------------------------------------------------------ | Coef. Std. Err. z P>|z| [95% Conf. Interval] -------------+---------------------------------------------------------------- wage | age | .1543231 .0194903 7.92 0.000 .1161227 .1925234 grade | .4225025 .029014 14.56 0.000 .3656362 .4793689 smsa | .8628628 .1285907 6.71 0.000 .6108297 1.114896 black | -.9206944 .1774617 -5.19 0.000 -1.268513 -.572876 tenure | .1003226 .051879 1.93 0.053 -.0013584 .2020037 union | 4.563859 1.006459 4.53 0.000 2.591236 6.536483 _cons | -4.670352 .5401517 -8.65 0.000 -5.72903 -3.611674 -------------+---------------------------------------------------------------- union | black | .4397974 .0972261 4.52 0.000 .2492377 .6303572 tenure | .0997638 .0236575 4.22 0.000 .053396 .1461317 south | -.4895032 .0933276 -5.24 0.000 -.6724221 -.3065844 _cons | -.9679795 .0746464 -12.97 0.000 -1.114284 -.8216753 -------------+---------------------------------------------------------------- hazard | lambda | -2.093313 .5801968 -3.61 0.000 -3.230478 -.9561486 -------------+---------------------------------------------------------------- rho | -0.89172 sigma | 2.3475104 ------------------------------------------------------------------------------ */第五个是手工完成两步估计法，估计结果与直接使用
etregress两步估计法基本保持一致。值得注意的是：union取值为1的样本和union取值为0的样本的IMR计算公式不同，因此如果混淆了公式或使用一个公式进行计算，最后的结果必然存在偏误，因为union取值为1的样本和union取值为0的样本均参与了第二步回归。- 手工两步法下参与回归的样本数目均为1,210，这点可以从两张回归结果表中直观看出。之所以相比于全样本缺少483个样本，原因在于有几个变量存在缺失值，Stata在回归时自动将这些变量存在缺失值的样本剔除，这也与前面四个模型的回归结果保持一致。
*-（五）手工两步估计法 *- 第一阶段方程（选择方程）回归 probit union $xlist1 south, nolog est store First /* Probit regression Number of obs = 1,210 LR chi2(3) = 56.54 Prob > chi2 = 0.0000 Log likelihood = -592.15536 Pseudo R2 = 0.0456 ------------------------------------------------------------------------------ union | Coef. Std. Err. z P>|z| [95% Conf. Interval] -------------+---------------------------------------------------------------- black | .4397974 .0972261 4.52 0.000 .2492377 .6303572 tenure | .0997638 .0236575 4.22 0.000 .053396 .1461317 south | -.4895032 .0933276 -5.24 0.000 -.6724221 -.3065844 _cons | -.9679795 .0746464 -12.97 0.000 -1.114284 -.8216753 ------------------------------------------------------------------------------ */ *- 计算IMR predict y_hat, xb gen imr = normalden(y_hat) / normal( y_hat) if union != . replace imr = -normalden(y_hat) / normal(-y_hat) if union == 0 *- 第二阶段方程回归 reg wage $xlist3 imr est store Second /* Source | SS df MS Number of obs = 1,210 -------------+---------------------------------- F(7, 1202) = 92.54 Model | 2227.31385 7 318.187693 Prob > F = 0.0000 Residual | 4132.88562 1,202 3.43834078 R-squared = 0.3502 -------------+---------------------------------- Adj R-squared = 0.3464 Total | 6360.19947 1,209 5.26071089 Root MSE = 1.8543 ------------------------------------------------------------------------------ wage | Coef. Std. Err. t P>|t| [95% Conf. Interval] -------------+---------------------------------------------------------------- union | 4.56386 .836918 5.45 0.000 2.921877 6.205842 age | .1543231 .0194468 7.94 0.000 .1161696 .1924765 grade | .4225025 .029448 14.35 0.000 .3647272 .4802778 smsa | .8628628 .12708 6.79 0.000 .6135395 1.112186 black | -.9206944 .1427409 -6.45 0.000 -1.200743 -.6406455 tenure | .1003226 .0424118 2.37 0.018 .0171133 .1835319 imr | -2.093314 .4858841 -4.31 0.000 -3.046589 -1.140038 _cons | -4.670352 .5337275 -8.75 0.000 -5.717493 -3.623211 ------------------------------------------------------------------------------ */最后是各个回归模型下
IMR的对比，以及回归结果的汇总输出。*- 对比imr、imr1和imr2 order imr1 imr2 imr union *- 回归结果输出 local mlist "OLS MLE MLE_r TwoStep First Second" esttab `mlist' using 处理效应模型.rtf, replace /// b(%6.4f) t(%6.4f) /// scalar(N r2_a r2_p) compress nogap /// star(* 0.1 ** 0.05 *** 0.01) /// mtitle(`mlist') /// title("Treatment Effects Model") /* Treatment Effects Model ---------------------------------------------------------------------------------------- (1) (2) (3) (4) (5) (6) OLS MLE MLE_r TwoStep First Second ---------------------------------------------------------------------------------------- main union 1.0039*** 4.5639*** 4.5639*** (7.5026) (4.5346) (5.4532) age 0.1476*** 0.1487*** 0.1487*** 0.1543*** 0.1543*** (7.5574) (7.6952) (7.1887) (7.9179) (7.9357) grade 0.4369*** 0.4206*** 0.4206*** 0.4225*** 0.4225*** (14.8228) (14.3256) (11.1284) (14.5620) (14.3474) smsa 0.9755*** 0.9117*** 0.9117*** 0.8629*** 0.8629*** (7.7872) (7.2992) (7.6015) (6.7102) (6.7899) black -0.6183*** -0.7882*** -0.7882*** -0.9207*** 0.4398*** -0.9207*** (-4.9388) (-5.7659) (-5.5944) (-5.1881) (4.5234) (-6.4501) tenure 0.2118*** 0.1524*** 0.1524*** 0.1003* 0.0998*** 0.1003** (6.2574) (4.1235) (3.2167) (1.9338) (4.2170) (2.3654) 1.union 2.9458*** 2.9458*** (10.7135) (5.2195) south -0.4895*** (-5.2450) imr -2.0933*** (-4.3083) _cons -4.3260*** -4.3516*** -4.3516*** -4.6704*** -0.9680*** -4.6704*** (-8.1386) (-8.2354) (-6.8302) (-8.6464) (-12.9675) (-8.7504) ---------------------------------------------------------------------------------------- union black 0.4557*** 0.4557*** 0.4398*** (4.7571) (4.8749) (4.5234) tenure 0.0872*** 0.0872*** 0.0998*** (3.7488) (3.4648) (4.2170) south -0.5807*** -0.5807*** -0.4895*** (-6.8233) (-6.9341) (-5.2450) _cons -0.8856*** -0.8856*** -0.9680*** (-12.2232) (-11.9748) (-12.9675) ---------------------------------------------------------------------------------------- / athrho -0.6544*** -0.6544*** (-7.1891) (-2.8179) lnsigma 0.7027*** 0.7027*** (23.9517) (9.2683) ---------------------------------------------------------------------------------------- hazard lambda -2.0933*** (-3.6079) ---------------------------------------------------------------------------------------- N 1210 1210 1210 1210 1210 1210 r2_a 0.3369 0.3464 r2_p 0.0456 ---------------------------------------------------------------------------------------- t statistics in parentheses * p<0.1, ** p<0.05, *** p<0.01 */五、已有文献的简单评述
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 i g = β 0 + β 1 S g + β 2 X i g + β 3 V g + ε (3) W_{ig}=\beta _0+\beta _1S_g+\beta _2X_{ig}+\beta _3V_g+\varepsilon \tag{3} Wig=β0+β1Sg+β2Xig+β3Vg+ε(3)
其中， W i g W_{ig} Wig表示第 g g g个村庄第 i i i个农民工在城市务工的工资； S g S_g Sg代表社会网络，用村庄在外务工的人数来表示，在论文中作为核心解释变量； X i g X_{ig} Xig表示个人层面的控制变量； V g V_g Vg是村庄层面的控制变量。
文章怀疑农民工在外务工可能是一个选择行为，因为具有城市劳动力市场优势（男性、年轻和能力强等）的农民会更倾向于外出打工。论文进一步将这种选择性的来源分为可观测因素和不可观测因素，可观测因素包括年龄、性别等，不可观测因素包括性格、能力等。因此，在模型设置时必须要考虑样本群体是否随机和均质，即样本选择偏差问题。
在本例的数据集结构中，被解释变量 W i g W_{ig} Wig存在缺失值，而由于不在外务工的农民本身就不具有在外务工的工资数据（即使他们有其他来源的收入，但不是文章研究的重点），因此这些缺失值存在的原因就是这些被调查的样本本身就不在外务工。而由于前面提到的农民是否在外务工可能是一个选择行为，即 W i g W_{ig} Wig存在缺失值的样本与 W i g W_{ig} Wig取值不为空的样本在某些特征因素方面本身就存在较大的差异，因此如果在回归时直接剔除这部分取值为空的样本，最后得到的结果就可能存在估计偏误，也就是说在考虑到样本选择偏差的情况下，基准回归结果可能就不再具有稳健性。
为了解决可能存在的样本选择偏差问题，作者使用了样本选择模型中的两步估计法（
Heckit）。Heckit由以下方程组构成：W i g = β 0 + β 1 S g + β 2 X i g + β 3 V g + β 4 P ^ i g + ε (4) W_{ig}=\beta _0+\beta _1S_g+\beta _2X_{ig}+\beta _3V_g+\beta _4\hat P_{ig}+\varepsilon \tag{4} Wig=β0+β1Sg+β2Xig+β3Vg+β4P^ig+ε(4)
P i g = γ 0 + γ 1 F i g + γ 2 S g + γ 3 X i g + γ 4 V g + μ (5) P_{ig}=\gamma _0+\gamma _1F_{ig}+\gamma _2S_g+\gamma _3X_{ig}+\gamma _4V_g+\mu \tag{5} Pig=γ0+γ1Fig+γ2Sg+γ3Xig+γ4Vg+μ(5)
其中，方程 ( 5 ) (5) (5)是第一阶段回归方程（选择方程），方程 ( 4 ) (4) (4)是第二阶段回归方程； P ^ i g \hat P_{ig} P^ig是逆米尔斯比率； P i g P_{ig} Pig代表样本是否外出务工的虚拟变量，即样本在外务工取值为1，否则为0； F i g F_{ig} Fig是选择方程中的外生解释变量，论文中选择的是家庭劳动力人数；方程 ( 4 ) (4) (4)的所有解释变量，是方程 ( 5 ) (5) (5)解释变量的严格子集（真子集）。
论文还对外生变量 F i g F_{ig} Fig的相关性与外生性进行了说明，认为家庭劳动力数量对农民工的打工决策有着重要影响，而对在外务工收入的影响微乎其微，具体分析请看原文。
事实上，这里存在两个问题：
- 第一，在选择方程中，作者使用logit模型进行回归，正如前文所述，logit模型不具备干扰项服从正态分布的假设，因此根据第一步回归拟合值计算出的
IMR可能存在一定程度的偏误。 - 第二，方程
(
4
)
(4)
(4)和方程
(
5
)
(5)
(5)对
P
^
i
g
\hat P_{ig}
P^ig或
P
i
g
P_{ig}
Pig的界定不清晰，因为方程
(
4
)
(4)
(4)中的
P
^
i
g
\hat P_{ig}
P^ig本意应该是
IMR，而非字母本身所表达的“ P i g P_{ig} Pig的拟合值”的含义，而IMR正是由“ P i g P_{ig} Pig的拟合值”计算所得。
下表是基准OLS回归与
Heckit第二步回归结果的对比，括号内为各变量估计系数的聚类稳健标准误，以村庄为聚类单位；Heckit的第一步回归结果论文并未列示，这里假定外生解释变量在第一步回归中显著且有效。(1)
OLS(2)
Heckit社会网络 0.125***
(0.0349)0.263***
(0.0760)个体、村庄
控制变量YES YES IMR — 0.754**
(0.3790)可以观察到，
Heckit第二步回归结果中IMR显著为正，且数值较大（相较于其他控制变量的估计系数而言，详细结果请看原文），这说明基准OLS回归确实存在样本选择偏差，造成估计偏误，具体来说是低估了社会网络对农民工在外务工收入的影响，因为OLS模型中社会网络的估计系数仅有0.125，而Heckit模型的估计系数（0.263）是其两倍还多，且两者均至少在1%的水平下显著。至于为什么社会网络的收入促进效应在Heckit模型中高于OLS，作者在原文中给出了解释。值得一提的是，为了解决一般性的因遗漏变量和联立方程（互为因果）导致的估计偏误问题，作者在样本选择模型的基础上进一步采用工具变量法，即采用
IV - Heckit对基准OLS的稳健性进行进一步的检验，详情请看原文。第二篇论文是祝树金和赵玉龙（2017）发表在《金融研究》的《资源错配与企业的出口行为——基于中国工业企业数据的经验研究》，主题是探讨企业资源错配对出口行为的影响。
[9] 祝树金, 赵玉龙. 资源错配与企业的出口行为——基于中国工业企业数据的经验研究[J]. 金融研究, 2017(11): 49-64.
论文一开始就考虑到了样本选择偏差问题，认为企业是否出口受制于自身条件，简单将出口企业与非出口企业同等对待将产生估计偏误，因此构建样本选择模型对这种偏误进行纠正。构建的两阶段模型如下：
e x d u m i t = α 1 m i s s i t − 1 + ∑ 1 J β j X j , i t − 1 + δ i t (6) exdum_{it}=\alpha _1miss_{it-1}+\sum _1^J\beta _jX_{j,it-1}+\delta_{it} \tag{6} exdumit=α1missit−1+1∑JβjXj,it−1+δit(6)
e x s h a r e i t = α 2 m i s s i t − 1 + ∑ 1 J β j Z j , i t − 1 + ε i t (7) exshare_{it}=\alpha _2miss_{it-1}+\sum _1^J\beta _jZ_{j,it-1}+\varepsilon_{it} \tag{7} exshareit=α2missit−1+1∑JβjZj,it−1+εit(7)
公式 ( 6 ) (6) (6)是第一阶段的企业出口选择方程，其中，被解释变量 e x d u m i t exdum_{it} exdumit代表第 t t t年 i i i企业的出口状态，若有出口行为，该变量记为1，否则记为0；考虑到出口滞后效应，第一、二阶段所有解释变量均滞后一期； m i s s i t − 1 miss_{it-1} missit−1是论文的核心解释变量企业资源错配，分别使用企业资源错配指标 m i s s 1 miss1 miss1、 τ k \tau k τk和 τ l \tau l τl来表示； X j , i t − 1 X_{j,it-1} Xj,it−1是影响企业出口决策的第 j j j个控制变量，包括一个外生解释变量 e x d u m i t − 1 exdum_{it-1} exdumit−1，即企业前一期是否出口的虚拟变量，作者认为该变量满足相关性和外生性的要求，以及其他控制变量。
公式 ( 7 ) (7) (7)是第二阶段的回归方程，其中，被解释变量 e x s h a r e i t exshare_{it} exshareit表示企业出口强度； Z j , i t − 1 Z_{j,it-1} Zj,it−1表示影响企业出口强度的控制变量，这些控制变量包括第一阶段的所有控制变量（除 e x d u m i t − 1 exdum_{it-1} exdumit−1），以及一个根据第一阶段回归拟合值计算的
IMR。回归结果汇总如下（限于篇幅限制，这里仅展示核心解释变量为 m i s s 1 miss1 miss1的回归结果；括号内为标准误，具体类型未告知）：
(1)
e x d u m exdum exdum(2)
e x s h a r e exshare exsharem i s s 1 t − 1 miss1_{t-1} miss1t−1 0.012***
(0.0030)0.008***
(0.0009)e x d u m t − 1 exdum_{t-1} exdumt−1 2.791***
(0.0099)— 控制变量 YES YES 行业、年份 YES YES W a l d c h i 2 Wald~chi^2 Wald chi2 11,949.25***(右同) ρ \rho ρ -0.3858***(右同) N N N 177,386 177,386 可以发现（ W a l d c h i 2 Wald~chi^2 Wald chi2和 ρ \rho ρ的那两列结果应该合并为一列，
markdown表格合并比较麻烦~），第一步选择方程（模型 ( 1 ) (1) (1)）中，外生变量的估计系数显著为正，说明外生变量的选择有效；第二阶段回归（模型 ( 2 ) (2) (2)）中，核心解释变量的估计系数显著为正，说明在考虑样本选择偏差的情况下，企业资源错配仍对企业出口强度产生促进作用。需要注意的是：- 模型 ( 1 ) (1) (1)和模型 ( 2 ) (2) (2)引入行业虚拟变量和年份虚拟变量，并在结果汇报中称之为行业（效应）和年份（效应），而非行业固定效应与年份固定效应，这种做法是比较严谨的，因为第一阶段的probit回归不能使用固定效应模型。
- 模型 ( 1 ) (1) (1)和模型 ( 2 ) (2) (2)使用的样本数是不同的，模型 ( 1 ) (1) (1)使用的样本数应该多于模型 ( 2 ) (2) (2)。然而，结果显示两步回归均使用了177,386个样本，因此猜测论文中的样本选择模型是使用MLE方法进行估计的，因为MLE从模型整体角度进行参数估计，参与回归的样本一般汇报为样本总数。
- 模型
(
2
)
(2)
(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 | .0237043 .0048317 4.91 0.000 .0142344 .0331742 roa | -.2253829 .0170078 -13.25 0.000 -.2587175 -.1920482 growth | .0353113 .0023365 15.11 0.000 .0307318 .0398909 toptenrate | .0001971 .0000445 4.43 0.000 .0001099 .0002844 dual | .0031043 .0015126 2.05 0.040 .0001396 .006069 independent | .0040687 .0115344 0.35 0.724 -.0185383 .0266758 magpay | .4069487 .1409151 2.89 0.004 .1307602 .6831373 boardshare | -.0090979 .0037856 -2.40 0.016 -.0165174 -.0016784 invrec | .0183075 .0051727 3.54 0.000 .0081693 .0284458 analyst | .0005457 .000082 6.65 0.000 .0003849 .0007065 opin | -.0389673 .0040999 -9.50 0.000 -.047003 -.0309316 aud | -.0038764 .0012816 -3.02 0.002 -.0063884 -.0013644 | year | 2009 | -.0008564 .0043254 -0.20 0.843 -.0093341 .0076213 2010 | .0186483 .0049932 3.73 0.000 .0088619 .0284348 2011 | .0220373 .0042826 5.15 0.000 .0136435 .030431 2012 | -.0158729 .0037267 -4.26 0.000 -.0231771 -.0085688 2013 | -.0161189 .0036 -4.48 0.000 -.0231748 -.009063 2014 | -.0161046 .0035979 -4.48 0.000 -.0231564 -.0090528 2015 | -.0031965 .0037801 -0.85 0.398 -.0106053 .0042122 2016 | -.0027102 .0037549 -0.72 0.470 -.0100697 .0046493 2017 | -.0153086 .0036145 -4.24 0.000 -.0223928 -.0082243 2018 | -.0154719 .0036009 -4.30 0.000 -.0225295 -.0084143 2019 | -.0178076 .0035699 -4.99 0.000 -.0248046 -.0108107 | industry | 2 | -.0297887 .0105573 -2.82 0.005 -.0504806 -.0090968 3 | -.0116159 .0109897 -1.06 0.291 -.0331553 .0099234 4 | -.0282726 .0115443 -2.45 0.014 -.0508991 -.0056462 5 | -.0434586 .0104201 -4.17 0.000 -.0638816 -.0230356 6 | -.0323764 .0107643 -3.01 0.003 -.0534741 -.0112787 | 1.coz1 | -.0067129 .0019325 -3.47 0.001 -.0105006 -.0029253 _cons | .2775419 .0214512 12.94 0.000 .2354982 .3195855 -------------+---------------------------------------------------------------- coz1 | Lagcoz1 | 3.234637 .041428 78.08 0.000 3.15344 3.315835 institution | .0037066 .0029546 1.25 0.210 -.0020843 .0094976 size | .1103075 .0245284 4.50 0.000 .0622327 .1583822 leverage | -.0127087 .1239879 -0.10 0.918 -.2557205 .2303031 roa | -.5487199 .328454 -1.67 0.095 -1.192478 .0950382 growth | .083781 .037123 2.26 0.024 .0110212 .1565408 toptenrate | .0032139 .0012742 2.52 0.012 .0007165 .0057114 dual | -.0820718 .0508186 -1.61 0.106 -.1816744 .0175308 independent | -1.012848 .3659513 -2.77 0.006 -1.730099 -.2955967 magpay | 2.53631 4.014366 0.63 0.528 -5.331703 10.40432 boardshare | -1.021367 .1736045 -5.88 0.000 -1.361625 -.6811082 invrec | -.1793844 .1334323 -1.34 0.179 -.440907 .0821381 analyst | -.0033357 .0025439 -1.31 0.190 -.0083216 .0016502 opin | .1161474 .1078862 1.08 0.282 -.0953057 .3276005 aud | .0256347 .0394008 0.65 0.515 -.0515895 .1028589 | year | 2009 | .0545774 .1108197 0.49 0.622 -.1626253 .2717801 2010 | .2064666 .1102673 1.87 0.061 -.0096534 .4225866 2011 | -.0294345 .1096026 -0.27 0.788 -.2442517 .1853827 2012 | .0080566 .1085939 0.07 0.941 -.2047835 .2208966 2013 | -.0643272 .1033437 -0.62 0.534 -.2668772 .1382227 2014 | -.0427743 .103273 -0.41 0.679 -.2451857 .1596371 2015 | -.0567423 .1105243 -0.51 0.608 -.273366 .1598813 2016 | -.0928139 .1030316 -0.90 0.368 -.2947522 .1091244 2017 | -.2477218 .106496 -2.33 0.020 -.45645 -.0389935 2018 | -.0962982 .1050334 -0.92 0.359 -.3021598 .1095634 2019 | .0256564 .0998903 0.26 0.797 -.1701251 .2214378 | industry | 2 | -.2286483 .2099405 -1.09 0.276 -.6401241 .1828274 3 | -.3286038 .2173325 -1.51 0.131 -.7545676 .09736 4 | -.4177377 .259037 -1.61 0.107 -.9254409 .0899654 5 | .069847 .2032102 0.34 0.731 -.3284377 .4681318 6 | .0932101 .210735 0.44 0.658 -.3198229 .5062431 | _cons | -4.316258 .5725048 -7.54 0.000 -5.438347 -3.194169 -------------+---------------------------------------------------------------- /athrho | .0413626 .0199663 2.07 0.038 .0022294 .0804958 /lnsigma | -2.451377 .0118703 -206.51 0.000 -2.474643 -2.428112 -------------+---------------------------------------------------------------- rho | .041339 .0199322 .0022294 .0803224 sigma | .0861748 .0010229 .0841931 .0882032 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 | .0237045 .0038194 6.21 0.000 .0162187 .0311904 roa | -.2253832 .0098235 -22.94 0.000 -.2446369 -.2061294 growth | .0353113 .0011144 31.69 0.000 .0331271 .0374955 toptenrate | .0001971 .000043 4.58 0.000 .0001128 .0002814 dual | .0031042 .0014823 2.09 0.036 .0001989 .0060095 independent | .0040678 .0114332 0.36 0.722 -.0183409 .0264765 magpay | .4069558 .1104819 3.68 0.000 .1904153 .6234964 boardshare | -.0090989 .0041262 -2.21 0.027 -.017186 -.0010117 invrec | .0183072 .0041006 4.46 0.000 .0102702 .0263441 analyst | .0005457 .0000836 6.53 0.000 .0003818 .0007096 opin | -.0389672 .0030457 -12.79 0.000 -.0449367 -.0329978 aud | -.0038764 .0012695 -3.05 0.002 -.0063645 -.0013882 | year | 2009 | -.0008564 .0038108 -0.22 0.822 -.0083254 .0066127 2010 | .0186484 .003792 4.92 0.000 .0112163 .0260805 2011 | .0220373 .003709 5.94 0.000 .0147678 .0293068 2012 | -.0158729 .0035715 -4.44 0.000 -.0228729 -.008873 2013 | -.0161189 .0035103 -4.59 0.000 -.022999 -.0092388 2014 | -.0161047 .0034962 -4.61 0.000 -.022957 -.0092523 2015 | -.0031967 .0035138 -0.91 0.363 -.0100836 .0036903 2016 | -.0027104 .0035086 -0.77 0.440 -.0095872 .0041665 2017 | -.0153088 .0034914 -4.38 0.000 -.0221519 -.0084658 2018 | -.0154721 .0034815 -4.44 0.000 -.0222958 -.0086485 2019 | -.0178078 .0034585 -5.15 0.000 -.0245864 -.0110293 | industry | 2 | -.0297886 .0066066 -4.51 0.000 -.0427372 -.0168399 3 | -.0116162 .0069175 -1.68 0.093 -.0251743 .001942 4 | -.0282727 .0074022 -3.82 0.000 -.0427807 -.0137646 5 | -.0434582 .0064804 -6.71 0.000 -.0561595 -.0307569 6 | -.0323759 .006855 -4.72 0.000 -.0458115 -.0189404 | coz1 | -.0067176 .0021868 -3.07 0.002 -.0110037 -.0024315 _cons | .2775378 .0181209 15.32 0.000 .2420215 .3130541 -------------+---------------------------------------------------------------- coz1 | Lagcoz1 | 3.235202 .0413025 78.33 0.000 3.154251 3.316153 institution | .0036154 .0029089 1.24 0.214 -.002086 .0093168 size | .1121405 .0239223 4.69 0.000 .0652536 .1590274 leverage | -.0208505 .1217953 -0.17 0.864 -.2595648 .2178638 roa | -.5578723 .3217238 -1.73 0.083 -1.188439 .0726948 growth | .0879266 .0303994 2.89 0.004 .0283449 .1475083 toptenrate | .0032304 .0013442 2.40 0.016 .0005958 .0058649 dual | -.079793 .050784 -1.57 0.116 -.1793278 .0197418 independent | -1.019378 .3750888 -2.72 0.007 -1.754539 -.2842178 magpay | 2.68412 3.683406 0.73 0.466 -4.535223 9.903462 boardshare | -1.018715 .1736654 -5.87 0.000 -1.359093 -.6783367 invrec | -.178663 .132241 -1.35 0.177 -.4378506 .0805247 analyst | -.0033943 .0025893 -1.31 0.190 -.0084692 .0016805 opin | .1167083 .0995624 1.17 0.241 -.0784304 .3118469 aud | .0260504 .040476 0.64 0.520 -.053281 .1053818 | year | 2009 | .0515999 .1140643 0.45 0.651 -.171962 .2751618 2010 | .2019266 .109724 1.84 0.066 -.0131284 .4169816 2011 | -.0315077 .1130668 -0.28 0.781 -.2531146 .1900991 2012 | .0078514 .1080623 0.07 0.942 -.2039468 .2196496 2013 | -.0658881 .108644 -0.61 0.544 -.2788264 .1470502 2014 | -.0454649 .1080572 -0.42 0.674 -.2572531 .1663232 2015 | -.0564083 .1072322 -0.53 0.599 -.2665795 .1537629 2016 | -.092317 .1085734 -0.85 0.395 -.3051169 .1204829 2017 | -.2501913 .1096922 -2.28 0.023 -.4651841 -.0351984 2018 | -.0990409 .1078524 -0.92 0.358 -.3104276 .1123459 2019 | .0223512 .1059781 0.21 0.833 -.1853622 .2300645 | industry | 2 | -.2381904 .204628 -1.16 0.244 -.6392539 .1628731 3 | -.3332723 .2179531 -1.53 0.126 -.7604526 .093908 4 | -.4154694 .2532426 -1.64 0.101 -.9118158 .080877 5 | .0617126 .1984843 0.31 0.756 -.3273095 .4507348 6 | .0875193 .2074342 0.42 0.673 -.3190443 .4940829 | _cons | -4.344332 .5487066 -7.92 0.000 -5.419777 -3.268887 -------------+---------------------------------------------------------------- hazard | lambda | .0035703 .001916 1.86 0.062 -.0001849 .0073256 -------------+---------------------------------------------------------------- rho | 0.04143 sigma | .08617484 ------------------------------------------------------------------------------ */ *- （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 | -.0208381 .0114122 -1.83 0.068 -.0432055 .0015293 institution | -.0001099 .0002544 -0.43 0.666 -.0006085 .0003886 size | .0065137 .0025792 2.53 0.012 .0014586 .0115689 leverage | .0287133 .0115701 2.48 0.013 .0060363 .0513904 roa | -.1696233 .0264688 -6.41 0.000 -.2215012 -.1177455 growth | .0406549 .0028657 14.19 0.000 .0350382 .0462716 toptenrate | .0005991 .0001313 4.56 0.000 .0003418 .0008563 dual | .004107 .0047744 0.86 0.390 -.0052507 .0134647 independent | -.1094472 .0341261 -3.21 0.001 -.1763332 -.0425612 magpay | .7141059 .3357257 2.13 0.033 .0560956 1.372116 boardshare | -.1155465 .030857 -3.74 0.000 -.1760251 -.055068 invrec | -.0200946 .0138371 -1.45 0.146 -.0472147 .0070255 analyst | -3.14e-06 .0002186 -0.01 0.989 -.0004316 .0004253 opin | -.0348136 .0083493 -4.17 0.000 -.051178 -.0184492 aud | .0026224 .0036027 0.73 0.467 -.0044387 .0096835 | year | 2009 | -.0133637 .0087203 -1.53 0.125 -.0304551 .0037278 2010 | -.0023983 .0085558 -0.28 0.779 -.0191672 .0143707 2011 | .007618 .0084863 0.90 0.369 -.0090149 .0242508 2012 | -.018136 .0082114 -2.21 0.027 -.03423 -.002042 2013 | -.0343626 .0082714 -4.15 0.000 -.0505743 -.0181509 2014 | -.0268081 .0082289 -3.26 0.001 -.0429364 -.0106798 2015 | .0043649 .0082808 0.53 0.598 -.0118652 .020595 2016 | -.0049588 .0082864 -0.60 0.550 -.0212 .0112823 2017 | -.0219129 .0084003 -2.61 0.009 -.0383771 -.0054486 2018 | -.0298953 .008332 -3.59 0.000 -.0462257 -.0135649 2019 | -.028234 .0081615 -3.46 0.001 -.0442303 -.0122377 | industry | 2 | -.0652461 .0191914 -3.40 0.001 -.1028605 -.0276318 3 | -.0279979 .0207284 -1.35 0.177 -.0686249 .012629 4 | -.0100713 .0240308 -0.42 0.675 -.0571708 .0370283 5 | -.0641525 .0188057 -3.41 0.001 -.101011 -.027294 6 | -.0458935 .0192651 -2.38 0.017 -.0836525 -.0081345 | _cons | -.1057979 .077478 -1.37 0.172 -.257652 .0460562 ---------------+---------------------------------------------------------------- coz2 | Laginstitution | .000229 .0017509 0.13 0.896 -.0032026 .0036607 Lagsize | .1220394 .0129369 9.43 0.000 .0966835 .1473953 Lagleverage | .3218294 .0709585 4.54 0.000 .1827532 .4609055 Lagroa | -.3267603 .204901 -1.59 0.111 -.7283589 .0748383 Laggrowth | -.0533554 .0214025 -2.49 0.013 -.0953036 -.0114072 Lagtoptenrate | .0049017 .0008002 6.13 0.000 .0033334 .0064701 Lagdual | -.1151072 .030836 -3.73 0.000 -.1755447 -.0546697 Lagindependent | -1.470211 .2309171 -6.37 0.000 -1.9228 -1.017622 Lagmagpay | 1.025986 2.201411 0.47 0.641 -3.2887 5.340672 Lagboardshare | -1.582939 .1017037 -15.56 0.000 -1.782274 -1.383604 Laginvrec | -.6018562 .0724373 -8.31 0.000 -.7438308 -.4598817 Laganalyst | -.0040897 .0015752 -2.60 0.009 -.0071769 -.0010024 Lagopin | .0822428 .0619452 1.33 0.184 -.0391676 .2036531 Lagaud | .0487422 .0246017 1.98 0.048 .0005237 .0969607 _cons | -3.507714 .2810605 -12.48 0.000 -4.058582 -2.956845 ---------------+---------------------------------------------------------------- /mills | lambda | .0948201 .0173606 5.46 0.000 .0607939 .1288462 ---------------+---------------------------------------------------------------- rho | 0.84225 sigma | .11257936 -------------------------------------------------------------------------------- */ *- （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 */c o z 3 coz3 coz3的基准回归及稳健性检验的第三段代码见公众号下载附件（公众号后台回复
heckman即可获取下载链接）。
-
-
新药临床试验中重复测量资料的混合效应模型-陈峰.pdf
2020-02-22 15:15:26用平均差值推论处理的作用。 随机区组设计：处理只能在区组内随机分配，每个实验单位接受的处理是不相同的；要求每个区组内实验单位彼此独立。 重复测量设计：“处理”是在区组（受试者）间随机分配，区组内的各... -
序列平均模型提高GPS控制测量基线解算精度的探讨
2020-07-12 21:41:18文中利用阳山金矿近几年采集的GPS控制测量野外数据,使用序列平均模型,采取4种方案进行基线解算数据比较,同时充分考虑卫星残差带来的影响,对序列平均模型提高精度的效果进行评估,发现序列平均模型是一个能提高GPS基线... -
【更新通知】手把手教你Stata软件操作与案例分析更新，速来！
2020-12-19 08:55:42继3大政策效应评价方法、面板微观计量模型、空间计量模型、应用面板数据模型四大主题套餐后，手把手教你Stata系列课程推出多期DID、平行趋势检验系列专题。该专题包含多期DID及平行趋势检验；双重差分(倍分法、倍差... -
拓端tecdat|R语言 线性混合效应模型实战案例
2019-06-13 08:46:48例如，多级模型本身可以称为分级线性模型，随机效应模型，多级模型，随机截距模型，随机斜率模型或汇集模型。根据学科，使用的软件和学术文献，许多这些术语可能指的是相同的一般建模策略。 读入数据 多级模型... -
单因素试验固定效应模型方差分析
2020-07-01 20:58:11单因素试验固定效应模型方差分析 观测值的线性模型 平方和与自由度分解 例题与SPSS求解 非平衡单因素试验SPSS求解 一、观测值的线性模型 单因素试验线性可加模型为： Yij为第i个处理的第j个观测值；U为所有观测值... -
【视频】线性混合效应模型(LMM,Linear Mixed Models)和R语言实现案例
2022-03-03 22:21:31视频：线性混合效应模型(LMM,Linear Mixed Models)和R语言实现 线性混合效应模型(LMM,Linear Mixed Models)和R语言实现案例 时长12:13 什么是混合效应建模，为什么要使用？ 统计分析中许多问题的传统方法. -
CATE符合ML-有条件的平均处理效果和机器学习-研究论文
2021-05-20 15:19:55结合计量经济学理论，他们不仅可以估计平均值，而且可以估计个性化的治疗效果-条件平均治疗效果（CATE）。 在本教程中，我们概述了新颖的方法，对其进行了详细说明，并通过Quantlets将其应用于实际数据应用程序中。... -
拓端tecdat|基于R语言混合效应模型（mixed model）案例研究
2019-06-17 17:16:25混合模型的输出将给出一个解释值列表，其效应值的估计值和置信区间，每个效应的p值以及模型拟合程度的至少一个度量。如果您有一个变量将您的数据样本描述为您可能收集的数据的子集，则应该使用混合模型而不是简单的... -
论文研究 - 基于自回归综合移动平均模型的改进广义自回归条件异方体模型在数据分析中的应用
2020-05-18 01:16:41其次，为了描述财务时间序列的高峰和肥尾以及杠杆效应，本文使用基于自回归综合移动平均模型的偏斜非对称幂自回归条件异方差模型来分析销售数据。 实证分析表明，考虑偏态分布的模型是有效的。 -
数据分析36计(24)：因果推断结合机器学习估计个体处理效应
2021-02-20 09:00:00个体异质性为何重要传统的因果推断分析，主要关注焦点是平均处理效应（Average Treatment Effect）。许多科学和工程都会面临这样的挑战，从个性化的医疗救治方案，到定制型的营... -
R语言线性混合效应模型实战案例
2019-04-08 17:23:00例如，多级模型本身可以称为分级线性模型，随机效应模型，多级模型，随机截距模型，随机斜率模型或汇集模型。根据学科，使用的软件和学术文献，许多这些术语可能指的是相同的一般建模策略。 读入数据 多级模型... -
R语言︱线性混合模型理论与案例探究（固定效应&随机效应）
2016-06-11 13:31:04线性混合模型与普通的线性模型不同的地方是除了有固定效应外还有随机效应。 ___________________________________________________________________________________ 一、线性混合模型理论 由两个部分来决定... -
拓端tecdat|R语言用lme4多层次（混合效应）广义线性模型（GLM），逻辑回归分析教育留级调查数据
2021-06-18 15:25:55本教程为读者提供了使用频率学派的广义线性模型（GLM）的基本介绍。具体来说，本教程重点介绍逻辑回归在二元结果和计数/比例结果情况下的使用，以及模型评估的方法。本教程使用教育数据例子进行模型的应用。此外，本... -
手把手教你Stata软件操作与案例分析
2020-12-19 08:55:461.5多重共线性 1.6异方差检验 1.7 WLS加权最小二乘法权重设置 1.8自相关检验 1.9一般化处理自相关或异方差 1.10离散选择模型（含Logit模型、Probit模型） 1.11一招搞定Logit、Probit离散选择模型边际效应 ... -
拓端tecdat|R语言混合效应逻辑回归（mixed effects logistic）模型分析肺癌数据
2021-04-22 18:45:08混合效应逻辑回归用于建立二元结果变量的模型，其中，当数据被分组或同时存在固定和随机效应时，结果的对数几率被建模为预测变量的线性组合。 混合效应逻辑回归的例子 例1：一个研究人员对40所不同大学的申请... -
利用序列平均的高精度GPS基线解算模型 (2012年)
2021-05-09 10:35:35系统分析了改进随机模型和改进函数模型两类GPS基线解算模型的优缺点，在此基础上，提出了一种基于序列平均的高精度GPS基线解算模型，即采用动态单历元技术进行静态基线解算，充分利用多路径效应的低频特性，采用小波... -
面板数据固定效应 vs. 随机效应
2021-01-12 20:25:26原标题：面板数据固定效应 vs. 随机效应来源：不止点滴 一般来说，经济数据有三种类型：横截面数据(包括混合横截面数据)、面板数据和时间序列数据。对于应用微观研究而言，主要还是采用前两种数据类型，时间序列数据... -
拓端tecdat|R语言混合线性模型、多层次模型、回归模型分析学生平均成绩GPA和可视化
2021-07-22 17:49:21原文出处：拓端数据部落公众号 混合模型在统计学领域已经存在了很长时间。例如，标准的方差分析方法可以被看作是混合模型的特殊情况。... 随机效应 随机系数 变化的系数 截距和/或斜率作为结果 分层线性模型 . -
线性混合效应模型Linear Mixed-Effects Models的部分折叠Gibbs采样
2018-08-15 17:55:00本文介绍了线性混合效应模型的新型贝叶斯分析。该分析基于部分折叠的方法，该方法允许某些组件从模型中部分折叠。得到的部分折叠的Gibbs（PCG）采样器被构造成适合线性...已经开发出混合效应模型来处理相关响应数... -
拓端tecdat|R语言经济学：动态模型平均(DMA)、动态模型选择(DMS)预测原油价格时间序列
2021-05-11 17:21:32简要地提供了在经济学中使用模型平均和贝叶斯方法的论据，使用了动态模型平均法（DMA），并与ARIMA、TVP等方法进行比较。希望对经济和金融领域的从业人员和研究人员有用。 动机 事实上，DMA将计量经济学建模的几个... -
具有相关关系的数据处理：线性混合模型与广义线性混合模型
2019-11-03 16:09:05进行数据分析时，会发现有时候一个模型中的变量之间可能具有相关性（correlation），比如面积和长度就具有高度的相关性，如果同时对这些参数建模，就存在共线性问题，所以一般是只针对其中一个参数建模。而这种... -
强可忽略处理分配下因果推断的结构回归模型关 (2000年)
2021-05-07 14:39:10建立强可忽略处理分配条件下因果推断的结构回归模型,估计平均处理效应．在正态分布假设下,总体参数的极大似然估计是渐近正态无偏估计．提出的方法可推广到具有测量误差的一般情形． -
拓端tecdat|R语言用线性混合效应（多水平/层次/嵌套）模型分析声调高低与礼貌态度的关系
2021-09-04 00:33:45线性混合效应模型与我们已经知道的线性模型有什么不同？ 线性混合模型（有时被称为 "多层次模型 "或 "层次模型"，取决于上下文）是一种回归模型，它同时考虑了（1）被感兴趣的自变量（如lm()）所解释的变化--固定... -
[线性模型总结] 线性回归+方差分析+协方差分析+混合效应+面板数据模型
2020-05-01 22:54:26线性模型是一类统计模型的总称，它包括了线性回归模型、方差分析模型、协方差分析模型和线性混合效应模型（或称方差分量模型）等。许多生物、医学、经济、管理、地质、气象、农业、工业、工程技术等领域的现象都可以...
