精华内容
下载资源
问答
  • 2021-08-07 19:24:13

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

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

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

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

    • 样本选择偏差。样本选择偏差的非随机选择机制在于对样本的选择不随机。在样本数据的采集过程中,只对某部分群体进行调查,但这部分群体与其他群体在某些方面的特征差异较大,因此根据这样的样本做回归得到的普适性结论并不可信。体现在具体的数据集中就是,数据集中只有特定群体的样本,或者,虽然有全部群体的所有解释变量数据,但除特定群体之外的其他群体的被解释变量数据缺失,在这两种情况下进行的回归,都将直接忽视其他群体的样本信息(y缺失的样本在参与回归时将被drop掉)。实质上,样本选择偏差说的就是参与回归的样本不能代表总体从而产生估计偏误的问题

    • 自选择偏差。自选择偏差的非随机选择机制在于对自变量的选择不随机。在使用DID方法评估政策效应时,一个明显的事实就是,相对于未实施政策的地区(控制组),实施政策的地区(处理组)通常情况下经济发展都较为发达、各类基础设施建设都较为完善,而所谓的“政策效果评估”也即考察政策的经济效应,因此地区是否参与政策这一行为是内生的。体现在回归方程中就是,经济指标(如,GDP、人均GDPGDP增长率等)作为被解释变量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只能解决依可测变量选择的内生问题,而将PSMDID结合(即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 估计思路的对比

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

    相同点在于:

    1. 都是两步估计法。Heckman于1979年提出的两步估计法最开始是用于解决样本选择偏差的,即最初的Heckman两步法指的就是样本选择模型,后来有学者借鉴这种两步估计法的思想,应用于解决自选择偏差的处理效应模型。这两个模型在估计思路上是一脉相承的,而正是因为这种相似性,所以才导致各个学者对这两个模型的错误理解与错误应用,这种错误在现阶段的文献中较为常见。

    2. 都可以使用MLE进行模型的整体估计。两步估计法(如2SLS、PSM - DID以及这里的样本选择模型和处理效应模型等)一个明显的缺陷是,第一步估计的误差将被带入第二步,导致效率损失。而使用MLE从整体上进行参数估计可以避免这种问题,但如果样本量过大,MLE估计耗时较长,且MLE对分布的假设较为严格,因此需要在估计的精准性、操作的简便性等方面进行权衡。

    3. 第一阶段回归都需要引入外生变量,同时应包括第二阶段的所有外生解释变量。引入的外生变量需满足相关性和外生性的要求,即与选择方程中的被解释变量在理论上和统计上均具有相关性,而与第二步回归的被解释变量不具有直接的相关关系。引入外生变量的目的是确保第一步计算得到的IMR在引入原回归方程后不与干扰项相关。该外生变量在处理效应模型中可以直接称作工具变量。此外,如果核心解释变量D是DID模型的did项,那么为了防止出现多重共线性,应该尽量避免使用D的滞后项D_lag作为工具变量。事实上,如果找到了一个良好的工具变量,也完全能够使用2SLS解决内生性问题。此外,两个模型除了都需要在第一阶段引入至少一个外生变量,第一阶段回归中的其余控制变量也应该是第二阶段回归中所有的控制变量,即应该包括所有的外生解释变量,原因在于保证两阶段估计的一致性,详情请看陈强教授的推文『工具变量法(五): 为何第一阶段回归应包括所有外生解释变量』。然而,部分文献在第一阶段并未包括第二阶段所有的外生解释变量,少部分文献甚至根本就不引入第二阶段的外生解释变量(如,考虑滞后效应,直接引入第二阶段外生解释变量的滞后项),并且在Stata处理效应模型的官方命令etregresshelp文件的演示案例中,第一阶段回归也并未包括所有的外生解释变量,原因可能在于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。在处理效应模型第一阶段回归方程中,被解释变量是原方程的核心解释变量DD取值为0或1,且不存在缺失值,该变量还同时参与了第二阶段回归,此外第一阶段引入的外生变量(或称工具变量)直接影响的是D

    3. 各阶段样本参与回归的数目不同。假设除关键变量,其余变量都不存在缺失值,那么对于样本选择模型来说,第一阶段回归的解释变量均不存在缺失值,被解释变量y_dummy取值为0或1,也不存在缺失值,因此选择方程中参与回归的样本是全样本,第二阶段由于被解释变量y本身就存在缺失值,因此参与第二阶段回归的样本不是全样本,从而第一阶段的样本多于第二阶段。对于处理效应模型来说,所有变量均不存在缺失值,因此两阶段参与回归的样本是相同的,虽然在第一阶段引入滞后项D_lag作为工具变量的情况下会损失一部分样本,但由于计算出来的IMR同样也存在缺失值,从而第二阶段参与回归的样本也将与第一阶段相同。

    4. 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以及控制变量x1x2,不需要写入IMR。选择项则包括对选择方程的具体设置,以及其他的细节。

    • select()表示写入选择方程,括号内的是选择方程的具体变量,有两种写法。
      • 第一种写法:直接写入控制变量(x1x2)和外生变量z1。此时严格要求原回归的被解释变量y存在缺失值,且缺失值不能以0作为标记。
      • 第二种写法:首先写入y是否被观测到的虚拟变量y_dummyy存在缺失值的样本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,解释变量是educage;选择方程中额外引入了两个外生解释变量marriedchildren

    由于被解释变量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,但显著性未知,该值等于rhosigma的乘积,其中:
      • 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()是标准正态的概率密度函数pdfnormal()是标准正态的累积分布函数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均有取值。
    • imrimr2保持高度一致(imr四舍五入了),因为两者都是用两步估计法计算出来的。
    • imr1imr2存在差距,但差别不大。
    *- 对比imr、imr1和imr2
    
    order imr1 imr2 imr wage_dummy
    

    最后是五个模型回归结果的对比,可以发现:

    • 两个解释变量(educationage)的回归系数的符号、大小、显著性在模型 ( 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中的官方命令是treatregetregress,其中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是核心解释变量,并且我们其怀疑存在自选择问题。agegradesmsablacktenure是基准回归方程的控制变量;blacktenure是选择方程的控制变量,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的三倍,这说明在未考虑自选择偏差的情况下,低估了unionwage的影响。
    • 其余解释变量估计系数的符号和显著性基本没变,但估计系数的数值改变较大,这说明自选择偏差带来的估计偏误不能忽视。
    • IMRlambda)的估计系数为-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显著,这说明原模型中的自选择偏差不可忽视,在未考虑自选择偏差的情况下,将严重低估unionwage的影响。
    *-(四)两步估计法
    
    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的取值不随机,也即本质上是自选择偏差;或者,样本选择偏差与自选择偏差这两种说法混着用。这样的做法也确实给读者造成了困扰,笔者之前阅读过多篇论文,确实感觉不同论文这对方面的说明大相径庭。而正是由于对问题的界定不清晰,所以造成了对模型使用的偏误。
    • 第二,混淆样本选择模型和处理效应模型。现阶段文献对样本选择模型和处理效应模型的实际操作主要有三种:heckmanetregress和手工两步法,其中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)
    个体、村庄
    控制变量
    YESYES
    IMR0.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=α1missit1+1JβjXj,it1+δ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=α2missit1+1JβjZj,it1+ε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} missit1是论文的核心解释变量企业资源错配,分别使用企业资源错配指标 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,it1是影响企业出口决策的第 j j j个控制变量,包括一个外生解释变量 e x d u m i t − 1 exdum_{it-1} exdumit1,即企业前一期是否出口的虚拟变量,作者认为该变量满足相关性和外生性的要求,以及其他控制变量。

    公式 ( 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,it1表示影响企业出口强度的控制变量,这些控制变量包括第一阶段的所有控制变量(除 e x d u m i t − 1 exdum_{it-1} exdumit1),以及一个根据第一阶段回归拟合值计算的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 exshare
    m i s s 1 t − 1 miss1_{t-1} miss1t10.012***
    (0.0030)
    0.008***
    (0.0009)
    e x d u m t − 1 exdum_{t-1} exdumt12.791***
    (0.0099)
    控制变量YESYES
    行业、年份YESYES
    W a l d   c h i 2 Wald~chi^2 Wald chi211,949.25***(右同)
    ρ \rho ρ-0.3858***(右同)
    N N N177,386177,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即可获取下载链接)。

    更多相关内容
  • 该方法可以克服逐步检验法和 Sobel 检验法在处理小样本量,小中介效应值,或者中介效应值不呈正态分布的情况下统计功效不高的缺点,并且能有效解决变量的测量误差以及多重中介模型的问题。在对该方法的原理进行介绍...
  • 基于结构方程模型多重中介效应分析
  • stata命令:中介效应分析
  • 固定效应模型

    万次阅读 多人点赞 2021-09-16 15:56:48
    (这里的主要原因是:若个体固定效应模型是采用Within回归(xtreg , fe),它会将不随时点变化的量都减去了,所以,如果模型中不随时点变化的虚拟变量(包括个体固定效应项)的属个数如果大于N(无截距项情形;...

    一、面板数据优点

    1. 可以解决遗漏变量的问题:遗漏变量由于不可观测的个体差异或“异质性”造成的,如果这种个体差异“不随时间而改变”,则面板数据提供了解决遗漏变量问题的又一利器。

    2. 提供更多个体动态行为的信息:由于面板数据同时有横截面与时间两个维度,优势它可以解决单独的截面数据或时间序列数据所不能解决的问题。

    3. 样本容量较大:由于同时有截面维度与时间维度,通常面板数据的样本容量更大,从而可以提高估计的精确度。

    估计面板数据长假定个体的回归方程拥有相同的斜率,但可以有不同的截距,以此来捕捉异质性

    扰动项由u_{i}+\varepsilon _{i,t}两部分构成,不可观测随机变量u_{i}代表个体异质性的截距项,\varepsilon _{i,t}为随个体与时间而改变的扰动项。假设\varepsilon _{i,t}为独立同分布且与u_{i}不相关。

    如果u_{i}与某个解释变量相关,则进一步称为固定效应模型;如果与所有解释变量均不相关,则称为随机变量模型 。

    二、混合回归

     如果个体都拥有完全一样的回归方程,则方程可以写成

     

    不包括常数项,就可以吧所有的数据放在一起,像对待横截面数据那样进行OLS回归。对于面板数据虽然可以假设所有的扰动项之间相互独立,但同一个体在不同时期的扰动项之间往往存在自相关。对标准误的估计应该使用聚类稳健的标准误,而所谓的聚类,就是每个个体不同时期的所有观测值所组成。同一聚类的观测值允许存在相关性,不同聚类的观测值则不相关。

    混合回归的基本假设是不存在个体效应。混合回归也被成为“总体平均估计量”可以理解为把个体效应都给平均掉了。

    三、个体固定效应模型

    组内离差——\beta _{FE}被称为“组内估计量”

    (1)即使个体特征u_{i}与解释变量x_{it}相关,只要使用组内估计量,就可以得到一致估计。

    (2)\beta _{FE}无法估计不随时间变化的变量的影响。

    (3)扰动项与各期解释变量均不相关,是比较强的假定。

    LSDV(Least Square Dummy Variable Model):如果在原方程中引入(n-1)个虚拟变量(如果没有截距项则引入n个虚拟变量)来代表不同的个体,可以得到相同的结果。但如果做完LSDV后删除不显著的个体虚拟变量则结果不一致。

    LSDV的优点是能够得到对个体异质性u_{i}的估计;缺点是如果n很大,则需要在方程中引入多个虚拟变量。

     四、时间固定效应

    个体固定效应模型解决了不随时间而变但随个体而异的遗漏变量问题,引入时间固定效应,则可以解决不随个体而变但随时间而变的遗漏变量问题。

     

    五、 行业固定效应

    在实证论文中我们通常可以看到文章没有做个体的固定效应,而是做了行业的固定效应。因为控制个体固定效应相当于为每一个个体(公司)都设置不同的截距项,而行业固定效应只需要为每一个行业设置不同的截距项,自然模型的要求低了很多。

    对于同时控制个体和行业

    控制个体时是随个体变化,不随时间变化的,所以,当你考虑其他不随时间变量的因素(行业、省份、区域、企业性质、银行所有制性质等)时,其实他们的信息都在个体中反应出来了,所以再设置不随时间变化的变量时,就是多余的了。(这里的主要原因是:若个体固定效应模型是采用Within回归(xtreg    , fe),它会将不随时点变化的量都减去了,所以,如果模型中不随时点变化的虚拟变量(包括个体固定效应项)的属个数如果大于N(无截距项情形;有截距项就是N-1个),它只能估计出前N个,其他的都不在模型中;若是采用LSDV法估计个体固定效应模型(reg     i.number),是设置了N-1个虚拟变量实现的,如果再往模型里加不随时点变化的虚拟变量(如行业、区域等),模型是会将它们排除在模型里面的。)所以,一些文献关于,在有个体固定效应的基础上,考虑控制(行业、省份、区域、企业性质、银行所有制性质等)这类不随时间变化的因素的影响时,不知道他们是如何控制的。

    除时间FE,其他非时变的FE均可由个体FE线性表出,这就意味着如果模型中控制多个非时变的FE,其他FE总能被个体FE表出,即存在多重共线性的问题,这样的FE将被omitted。因此,许多论文不会在模型中同时控制个体FE和行业FE。

    然而,这并不是说同时控制个体FE和行业FE是不可行的。一种特殊情况是,如果企业所属行业发生变更(如环境规制政策实施前后,部分制造业企业选择变更行业以规避政策的不利影响或套取政策红利,虽然后一种情况比较少),在这种情况下行业FE将不再是非时变的了,因此行业FE就不会再被个体FE线性表出。况且,就算不存在企业跨行转移的情况,也可以通过附上时变因素来规避共线性的问题。

    六、交互固定效应

    年份FE的同质性就是假定在同一年份某一不可观测因素(如政策冲击、经济周期等)对所有企业的结果变量y的作用方向、作用大小是一样的。但是,现实的经济冲击并不会对所有企业产生一致的同质性影响,不同企业因自身实力、价值链地位、所有者性质等的不同在面对同一经济冲击时做出的战略性反应不同,从而导致最终的结果不同。

    行业FE表征企业所属行业的不可观测的典型特征对企业的同质性影响,换言之,如果怀疑行业的某些特征对行业内所有企业的y均存在影响(如金融业企业一般都比较“赚钱”),并且对行业内的不同企业的作用大小不存在明显差异,那么行业FE就可以代表这样的行业特征。行业FE假定同一行业中的样本行业特征是近似一致的。(这一假定从数据集中也可以看出来,即同一行业样本的indfe#均赋值为1(属于行业#),或者均赋值为0(不属于行业#)。其他FE同理。)

    总结来说就是,控制时间FE仅仅考虑到了时间维度上的同质性经济冲击,但现实中的经济冲击将对不同类型企业产生异质性影响,为将这些不可观测的异质性冲击因素控制住,回归方程需要引入交互FE。

    交互FE比单独的FE更严格,交互FE本质上包含了单个FE,引入过多的虚拟变量可能导致核心解释变量统计上不显著,甚至造成符号与预期相反,这种情况下就需要仔细斟酌一下,到底是经济系统本身就是这种运行规律?还是说过多的虚拟变量导致某些控制变量被omitted,从而影响了估计结果?切不能简单地“见Star行事”,因为某些情况下基于这样的交互FE得出的结果更能反映经济系统本身的运行规律,且不显著的回归结果某种程度上可以讨论出影响机制,增强论文的故事性,比如分样本回归。

    但是,有一种情况建议使用交互FE。假设基于某个政策做一个DID,“两高一剩”行业企业treated赋值为1,其他企业赋值为0;2012年及以后post赋值为1,以前赋值为0;被解释变量是企业TFP。观察这一模型的数据结构可以发现,被解释变量是企业级别,核心解释变量是行业 - 年份级别。那么,为了控制企业级别的不可观测因素对企业TFP的影响,同时为了控制样本期间其他所有行业级别的环境规制政策对企业TFP的影响,模型就需要引入企业FE和行业 - 年份FE,至于行业代码具体细化到什么程度,这就是另外的故事了。

    七、面板数据的stata命令

    1. 面板数据设定

    xtset panelvar timevar
    xtdes    //显示面板数据的结构,是否为平衡面板
    xtsum    //显示组内、组间与整体的统计指标
    xttab varname     //显示组内、组间与整体的分布频率,tab指的是tabulate
    xtline varname    //对每个个体分别显示该变量的时间序列图;如果希望所有个体的时间序列图叠放在一起,可加上选择项overlay

    2.  混合回归

    reg y x1 x2 x3,vce(cluster id)         //聚类标准误估计准确

    3. 固定效应

    xtreg y x1 x2 x3,fe r
    *-LSDV法的stata命令:
    reg y x1 x2 x3 i.id, r    
    //其中r表示聚类稳健标准误,使用vce(cluster id)能达到一样的效果
     i.id则表示根据变量id而生成的虚拟变量

    (1)如果通过LSDV法得到大多数个体虚拟变量都很显著,即认为存在个体效应,不应该使用混合回归。

    (2)双向固定效应模型

    xtreg y x1 x2 x3 i.year,fe r

    (3)交互固定效应

    regife ln_w tenure, f(id year, 1)                   //考虑一维交互固定效应
    regife ln_w tenure, a(id) f(id year, 1)             //加入个体固定效应
    regife ln_w tenure, a(id year) f(id year, 1)        // 加入个体和时间固定效应
    regife ln_w tenure, f(fid = id fyear = year, 1)     //生成因子载荷和共同因子并保存在新变量fid、fyear中
    regife ln_w tenure, f(id year, 1) residuals(newvar) //保存残差项

    eg.加入省份*年固定效应

    encode province,gen(province2) //首先对字符变量解码成数值变量
    regife ln_Cash_ratio1 lnnumber_in_5km $control, a(stkcd year) ife(province2 year, 1)

    rehdfe也可以做交互效应

    reghdfe ln_w grade age ttl_exp tenure not_smsa , absorb(idcode#occ)
    reghdfe lny x1 x2,absorb(i.region#i.industry i.region#i.year i.industry#i.year)
    
    gen plus=province2*year
    reghdfe ln_Cash_ratio1 lnnumber_in_5km $control, a(stkcd year plus) 
    

    (1)交互的时候要注意考虑维度限制,即自由度够不够。如果是城市面板,那用省和年交互,可能自由度才够。但是如果用城市id和year都比较大,那么会产生非常大的矩阵,自由度也不够,使得Stata无法运算

    (2)加了 id#year,没加单独的;还是加了单独,都是控制year和id的特征,核心变量系数一般不会这两种选择而不同,对结果影响不大。但是如果你的 id 和 year 都比较大,那么运算速度是一个问题。

    (3)regife 与 reghdfe 原理是不一样的

    4. 豪斯曼检验

    xtreg y x1 x2 x3,fe       //固定效应估计
    estimates store FE
    xtreg y x1 x2 x3,re        //随机效应估计
    estimates store RE
    hausman FE RE,constant sigmamore          //豪斯曼检验

    八、聚类(解决自相关问题)

    加入个体层面的固定效应将使得残差项里不再包括那些不随时间变化的遗漏变量,然而,残差项的方差是否一致,是否自相关的问题都未得到解决。而聚类稳健标准误(比如你对公司聚类)则使得公司层面的异方差问题得到解决。

    即固定效应的目的是解决内生性问题,为了使得系数估计无偏且一致,使得因果识别更干净。可以理解为一组控制变量。稳健标准误的目的是为了系数估计量估计的准确性,系数估计量的准确性会取决于扰动项的相关性,因此不同的标准误的形式对应着关于不同的个体之间扰动项相关性的假设

    reg y1 x1 x2 x3 x4,vce(cluster state)      //"state"为聚类变量的聚类稳健标准误
    reg y1 x1 x2 x3 x4                        //普通标准误

    个体变量在不同时期之间的扰动项一般会存在自相关,磨人的普通标准无计算方法假设扰动项为独立同分布的,故普通标准误的估计并不准确。

    PS(只要使用了稳健标准误,就可以和异方差和平共处了)

    reg y x,robust

    展开全文
  • 根据多个中介变量之间是否存在相互影响,多重中介模型可以分为单步多重中介模型(single-step mul-tiple mediator model)和多步多重中介模型(multiple-step multiple mediator model)。单步多重中介模型是指多个...

    1、数据来源:自主整理

    2、时间跨度:无

    3、区域范围:无

    4、指标说明:

    多重中介模型即存在多个中介变量的模型。根据多个中介变量之间是否存在相互影响,多重中介模型可以分为单步多重中介模型(single - step mul-tiple mediator model)和多步多重中介模型(multiple- step multiple mediator model)。单步多重中介模型是指多个中介变量之间不存在相互影响(图2 中去掉路径便是一个单步多重中介模型),又称为并行多重中介模型。多步多重中介模型是指多个中介变量之间存在相互影响,多个中介变量表现出顺序性特征,形成中介链(如图 2 中的 路径),又称为链式多重中介模型。

    以图 2 所示的含有两个中介变量 和  的多重中介模型为例(其中的圆圈表示变量可以是显变量或潜变量),此时的多重中介效应分析可以从三个角度进行,(1)特定路径的中介效应(specific mediationeffect),如  、 和  ,(2)总的中介效应(total mediation effect),即  ,(3)对比中介效应,如   和  。

    多重中介模型相对于简单中介模型具有三大优势。首先,可以得到总的中介效应。其次,可以在控制其他中介变量(如控制 )的前提下,研究每个中介变量(如 )的特定中介效应。这种做法可以减少简单中介模型因为忽略其他中介变量而导致的参数估计偏差。第三,可以得到对比中介效应,使得研究者能判断多个中介变量的效应(如  和  )中,哪一个效应更大,即判断哪一个中介变量的作用更强。这样,对比中介效应能使研究者判断多个中介变量理论(如 和 )中,哪个中介变量理论更有意义。因此研究多重中介模型更具理论和实践意义。

    部分代码如下:

    相关研究:

    [1]张正堂. 人力资源管理活动与企业绩效的关系:人力资源管理效能中介效应的实证研究[J]. 经济科学, 2006(2):11.

    [2]解学梅, 左蕾蕾. 企业协同创新网络特征与创新绩效:基于知识吸收能力的中介效应研究[J]. 南开管理评论, 2013, 16(3):10.

    [3]刘云, 石金涛. 组织创新气氛对员工创新行为的影响过程研究——基于心理授权的中介效应分析[J]. 中国软科学, 2010(3):12.

    [4]郭国庆, 张中科, 陈凯,等. 口碑传播对消费者品牌转换意愿的影响:主观规范的中介效应研究[J]. 管理评论, 2010, 22(12):8.

    下载链接多重中介STATA代码.zip-数据集文档类资源-CSDN下载

    展开全文
  • 建立了空气孔径呈现横向线性变化的三角形网格梯度折射率光子晶体平板透镜二维模型,应用多重散射方法(MSM),对该模型在波长平面光入射的情况下, TE模式的电磁场分布进行了数值模拟,从而验证了该光子晶体透镜的聚焦...
  • 面板模型混合效应模型This article shows how tree-boosting (sometimes also referred to as “gradient tree-boosting”) can be combined with mixed effects models using the GPBoost algorithm. Background is...

    面板模型混合效应模型

    This article shows how tree-boosting (sometimes also referred to as “gradient tree-boosting”) can be combined with mixed effects models using the GPBoost algorithm. Background is provided on both the methodology as well as on how to apply the GPBoost library using Python. We show how (i) models are trained, (ii) parameters tuned, (iii) model are interpreted, and (iv) predictions are made. Further, we do a comparison of several alternative approaches.

    本文展示了如何使用GPBoost算法将树增强(有时也称为“梯度树增强”)与混合效果模型结合使用。 提供了方法论以及如何使用Python应用GPBoost库背景知识 。 我们展示了如何(i)训练模型,(ii)调整参数,(iii)解释模型,以及(iv)进行预测。 此外,我们对几种替代方法进行了比较。

    介绍 (Introduction)

    Tree-boosting with its well-known implementations such as XGBoost, LightGBM, and CatBoost, is widely used in applied data science. Besides state-of-the-art predictive accuracy, tree-boosting has the following advantages:

    借助XGBoost,LightGBM和CatBoost等著名的实现来增强树性能在应用数据科学中得到了广泛的应用。 除具有最新的预测准确性外,增强树功能还具有以下优点:

    • Automatic modeling of non-linearities, discontinuities, and complex high-order interactions

      自动建模非线性,不连续和复杂的高阶相互作用
    • Robust to outliers in and multicollinearity among predictor variables

      稳健的预测变量中的异常值和多重共线性
    • Scale-invariance to monotone transformations of the predictor variables

      尺度不变性到预测变量的单调变换
    • Automatic handling of missing values in predictor variables

      自动处理预测变量中的缺失值

    Mixed effects models are a modeling approach for clustered, grouped, longitudinal, or panel data. Among other things, they have the advantage that they allow for more efficient learning of the chosen model for the regression function (e.g. a linear model or a tree ensemble).

    混合效果模型是针对聚类,分组,纵向或面板数据的建模方法。 除其他优点外,它们还具有以下优点:允许更有效地学习所选的回归函数模型(例如线性模型或树集合)。

    As outlined in Sigrist (2020), combined gradient tree-boosting and mixed effects models often performs better than (i) plain vanilla gradient boosting, (ii) standard linear mixed effects models, and (iii) alternative approaches for combing machine learning or statistical models with mixed effects models.

    Sigrist(2020)所述, 结合梯度树增强和混合效应模型的性能通常比(i)普通香草梯度增强,(ii)标准线性混合效应模型和(iii)结合机器学习或统计的替代方法要好。具有混合效果模型的模型。

    建模分组数据 (Modeling grouped data)

    Grouped data (aka clustered data, longitudinal data, panel data) occurs naturally in many applications when there are multiple measurements for different units of a variable of interest. Examples include:

    当对感兴趣变量的不同单位进行多次测量时,在许多应用程序中自然会出现分组数据(又名聚类数据,纵向数据,面板数据) 。 示例包括:

    • One wants to investigate the impact of some factors (e.g. learning technique, nutrition, sleep, etc.) on students’ test scores and every student does several tests. In this case, the units, i.e. the grouping variable, are the students and the variable of interest is the test score.

      一个人想调查某些因素(例如学习技术,营养,睡眠等)对学生考试成绩的影响,而每个学生都进行几次考试。 在这种情况下,单位,即分组变量,是学生,而感兴趣的变量是测试分数。
    • A company gathers transaction data about its customers. For every customer, there are several transactions. The units are then the customers and the variable of interest can be any attribute of the transactions such as prices.

      公司收集有关其客户的交易数据。 对于每个客户,都有几笔交易。 单位就是客户,兴趣变量可以是交易的任何属性,例如价格。

    Basically, such grouped data can be modeled using four different approaches:

    基本上,可以使用四种不同的方法对此类分组数据进行建模:

    1. Ignore the grouping structure. This is rarely a good idea since important information is neglected.

      忽略分组结构 。 因为忽略了重要信息,所以这很少是一个好主意。

    2. Model each group (i.e. each student or each customer) separately. This is also rarely a good idea as the number of measurements per group is often small relative to the number of different groups.

      分别对每个小组(即每个学生或每个客户)建模 。 这也不是一个好主意,因为每组的测量数量相对于不同组的数量通常很小。

    3. Include the grouping variable (e.g. student or customer ID) in your model of choice and treat it as a categorical variable. While this is a viable approach, it has the following disadvantages. Often, the number of measurements per group (e.g. number of tests per student, number of transactions per customer) is relatively small and the number of different groups is large (e.g. number of students, customers, etc.). In this case, the model needs to learn many parameters (one for every group) based on relatively little data which can make the learning inefficient. Further, for trees, high cardinality categorical variables can be problematic.

      在您选择的模型中包括分组变量(例如,学生或客户ID),并将其视为分类变量。 尽管这是一种可行的方法,但它具有以下缺点。 通常,每组的测量数量(例如,每个学生的测试数量,每个客户的交易数量)相对较小,而不同组的数量却很大(例如,学生,客户数量等)。 在这种情况下,模型需要基于相对较少的数据来学习许多参数(每组一个),这会使学习效率低下。 此外,对于树木,高基数类别变量可能会出现问题。

    4. Model the grouping variable using so-called random effects in a mixed effects model. This is often a sensible compromise between the approaches 2. and 3. above. In particular, as illustrated below and in Sigrist (2020), this is beneficial compared to the other approaches in the case of tree-boosting.

      在混合效应模型中使用所谓的随机效应对分组变量进行建模。 这通常是上述方法2和方法3之间的明智折衷。 尤其是,如下面和Sigrist(2020)中所示,与在树增强情况下的其他方法相比这是有益的。

    方法论背景 (Methodological background)

    For the GPBoost algorithm, it is assumed that the response variable y is the sum of a non-linear mean function F(X) and so-called random effects Zb:

    对于GPBoost算法,假定响应变量y是非线性均值函数F(X)与所谓的随机效应Zb的和

    y = F(X) + Zb + e

    y = F(X)+ Zb + e

    where

    哪里

    • y the response variable (aka label)

      y响应变量(也称为标签)
    • X contains the predictor variables (aka features) and F() is a potentially non-linear function. In linear mixed effects models, this is simply a linear function. In the GPBoost algorithm, this is an ensemble of trees.

      X包含预测变量(又称特征),F()是潜在的非线性函数。 在线性混合效果模型中,这只是线性函数。 在GPBoost算法中,这是一棵树木。
    • Zb are the random effects which are assumed to follow a multivariate normal distribution

      Zb是假定遵循多元正态分布的随机效应
    • e is an error term

      e是错误项

    The model is trained using the GPBoost algorithm, where trainings means learning the (co-)variance parameters (aka hyper-parameters) of the random effects and the regression function F(X) using a tree ensemble. The random effects Zb can be estimated (or predicted, as it is often called) after the model has been learned. In brief, the GPBoost algorithm is a boosting algorithm that iteratively learns the (co-)variance parameters and adds a tree to the ensemble of trees using a gradient and/or a Newton boosting step. The main difference to existing boosting algorithms is that, first, it accounts for dependency among the data due to clustering and, second, it learns the (co-)variance of the random effects. See Sigrist (2020) for more details on the methodology. In the GPBoost library, (co-)variance parameters can be learned using (accelerated) gradient descent or Fisher scoring, and trees are learned using the LightGBM library. In particular, this means that the full functionality of LightGBM is available.

    使用GPBoost算法对模型进行训练,其中训练意味着 使用树集合 学习 随机效应 的(共)方差参数(aka超参数) 和回归函数F(X) 。 在学习模型之后,可以估计(或预测,通常称为)随机效应Zb。 简而言之,GPBoost算法是一种增强算法,它迭代地学习(协)方差参数,并使用梯度和/或牛顿增强步骤将一棵树添加到树的集合中。 与现有增强算法的主要区别在于,首先,它考虑了由于聚类导致的数据之间的依赖性,其次,它学习了随机效应的(协)方差。 有关该方法的更多详细信息,请参见Sigrist(2020) 。 在GPBoost库中,可以使用(加速)梯度下降或Fisher评分来学习(协)方差参数,而可以使用LightGBM库来学习树。 特别是,这意味着可以使用LightGBM的全部功能。

    如何在Python中使用GPBoost库 (How to use the GPBoost library in Python)

    In the following, we show how combined tree-boosting and mixed effects models can be applied using the GPBoost library from Python. Note that there is also an equivalent R package. More information on this can be found here.

    在下面的内容中,我们展示了如何使用Python的GPBoost库应用组合的树加速和混合效果模型。 请注意,还有一个等效的R包。 有关此的更多信息,请参见此处

    安装 (Installation)

    pip install gpboost -U

    模拟数据 (Simulate data)

    We use simulated data here. We adopt a well known non-linear function F(X). For simplicity, we use one grouping variable. But one could equally well use several random effects including hierarchically nested ones, crossed ones, or random slopes. The number of samples is 5'000 and the number of different groups or clusters is 500. We also generate test data for evaluating the predictive accuracy. For the test data, we include both known, observed groups as well as novel, unobserved groups.

    我们在这里使用模拟数据。 我们采用了众所周知的非线性函数F(X) 。 为简单起见,我们使用一个分组变量。 但是同样可以很好地使用几种随机效果,包括层次嵌套的效果,交叉效果或随机斜率。 样本数量为5,000,不同组或群集的数量为500。我们还生成测试数据以评估预测准确性。 对于测试数据,我们既包括已知的观察组,也包括新颖的未观察组。

    import gpboost as gpb
    import numpy as np
    import sklearn.datasets as datasets
    import time
    import pandas as pd# Simulate data
    ntrain = 5000 # number of samples for training
    n = 2 * ntrain # combined number of training and test data
    m = 500 # number of categories / levels for grouping variable
    sigma2_1 = 1 # random effect variance
    sigma2 = 1 ** 2 # error variance
    # Simulate non-linear mean function
    np.random.seed(1)
    X, F = datasets.make_friedman3(n_samples=n)
    X = pd.DataFrame(X,columns=['variable_1','variable_2','variable_3','variable_4'])
    F = F * 10**0.5 # with this choice, the fixed-effects regression function has the same variance as the random effects
    # Simulate random effects
    group_train = np.arange(ntrain) # grouping variable
    for i in range(m):
    group_train[int(i * ntrain / m):int((i + 1) * ntrain / m)] = i
    group_test = np.arange(ntrain) # grouping variable for test data. Some existing and some new groups
    m_test = 2 * m
    for i in range(m_test):
    group_test[int(i * ntrain / m_test):int((i + 1) * ntrain / m_test)] = i
    group = np.concatenate((group_train,group_test))
    b = np.sqrt(sigma2_1) * np.random.normal(size=m_test) # simulate random effects
    Zb = b[group]
    # Put everything together
    xi = np.sqrt(sigma2) * np.random.normal(size=n) # simulate error term
    y = F + Zb + xi # observed data
    # split train and test data
    y_train = y[0:ntrain]
    y_test = y[ntrain:n]
    X_train = X.iloc[0:ntrain,]
    X_test = X.iloc[ntrain:n,]

    学习和做出预测 (Learning and making predictions)

    The following code shows how one trains a model and makes predictions. As can be seen below, the learned variance parameters are close to the true ones. Note that when making predictions, one can make separate predictions for the mean function F(X) and the random effects Zb.

    以下代码显示了如何训练模型并进行预测。 如下所示,学习的方差参数接近真实参数。 注意,进行预测时,可以对均值函数F(X)和随机效应Zb进行单独的预测。

    # Define and train GPModel
    gp_model = gpb.GPModel(group_data=group_train)
    # create dataset for gpb.train function
    data_train = gpb.Dataset(X_train, y_train)
    # specify tree-boosting parameters as a dict
    params = { 'objective': 'regression_l2', 'learning_rate': 0.1,
    'max_depth': 6, 'min_data_in_leaf': 5, 'verbose': 0 }
    # train model
    bst = gpb.train(params=params, train_set=data_train, gp_model=gp_model, num_boost_round=32)
    gp_model.summary() # estimated covariance parameters
    # Covariance parameters in the following order:
    # ['Error_term', 'Group_1']
    # [0.9183072 1.013057 ]
    # Make predictions
    pred = bst.predict(data=X_test, group_data_pred=group_test)
    y_pred = pred['fixed_effect'] + pred['random_effect_mean'] # sum predictions of fixed effect and random effect
    np.sqrt(np.mean((y_test - y_pred) ** 2)) # root mean square error (RMSE) on test data. Approx. = 1.25

    参数调整 (Parameter tuning)

    A careful choice of the tuning parameters is important for all boosting algorithms. Arguably the most important tuning parameter is the number of boosting iterations. A too large number will often result in over-fitting in regression problems and a too small value in “under-fitting”. In the following, we show how the number of boosting iterations can be chosen using cross-validation. Other important tuning parameters include the learning rate, the tree-depth, and the minimal number of samples per leaf. For simplicity, we do not tune them here but use some default values.

    仔细选择调整参数对于所有升压算法都很重要。 可以说,最重要的调整参数是加速迭代的次数。 数量太大通常会导致回归问题过度拟合,而“欠拟合”值太小。 在下面,我们展示了如何使用交叉验证来选择增强迭代的次数。 其他重要的调整参数包括学习率,树深度和每片叶子的最少样本数。 为简单起见,我们在这里不对其进行调整,而是使用一些默认值。

    # Parameter tuning using cross-validation (only number of boosting iterations)
    gp_model = gpb.GPModel(group_data=group_train)
    cvbst = gpb.cv(params=params, train_set=data_train,
    gp_model=gp_model, use_gp_model_for_validation=False,
    num_boost_round=100, early_stopping_rounds=5,
    nfold=4, verbose_eval=True, show_stdv=False, seed=1)
    best_iter = np.argmin(cvbst['l2-mean'])
    print("Best number of iterations: " + str(best_iter))
    # Best number of iterations: 32

    特征重要性和部分依赖图 (Feature importance and partial dependence plots)

    Feature importance plots and partial dependence plots are tools for interpreting machine learning models. These can be used as follows.

    特征重要性图和偏相关图是解释机器学习模型的工具。 这些可以如下使用。

    # Plotting feature importances
    gpb.plot_importance(bst)
    Image for post
    Feature importance plot
    特征重要性图

    Univariate partial dependence plots

    单变量偏相关图

    from pdpbox import pdp
    # Single variable plots
    pdp_dist = pdp.pdp_isolate(model=bst, dataset=X_train,
    model_features=X_train.columns,
    feature='variable_2',
    num_grid_points=100)
    pdp.pdp_plot(pdp_dist, 'variable_2', plot_lines=True)
    Image for post
    Partial dependence plot for variable 2
    变量2的偏相关图

    Multivariate partial dependence plots

    多元偏相关图

    # Two variable interaction plot
    inter_rf = pdp.pdp_interact(model=bst, dataset=X_train, model_features=X_train.columns,
    features=['variable_1','variable_2'])
    pdp.pdp_interact_plot(inter_rf, ['variable_1','variable_2'], x_quantile=True, plot_type='contour', plot_pdp=True)
    Image for post
    Two dimensional partial dependence plot for visualizing interactions
    二维局部依赖图,用于可视化交互

    SHAP值 (SHAP values)

    SHAP values and dependence plots are another important tool for model interpretation. These can be created as follows.

    SHAP值和依赖性图是模型解释的另一个重要工具。 这些可以如下创建。

    Edit: this is currently not yet fully supported by the shap Python package. It should be available soon (hopefully in the next days, see here for the current status). In the meantime, you have to copy-paste a few lines of code to your shap Python package. Just go to the location where your python packages are and add these green marked lines of code to the shap/tree_explainers/tree.py file.

    编辑:shap Python软件包目前尚未完全支持此功能。 它应该很快就可用(希望在接下来的几天中,请参阅 此处 了解当前状态)。 同时,您必须将几行代码复制粘贴到您的shap Python包中。 只需转到python包所在的位置,然后将 这些 带有 绿色标记的代码行 添加 到shap / tree_explainers / tree.py文件即可。

    import shap
    shap_values = shap.TreeExplainer(bst).shap_values(X_test)
    shap.summary_plot(shap_values, X_test)
    shap.dependence_plot("variable_2", shap_values, X_test)
    Image for post
    SHAP values
    SHAP值
    Image for post
    SHAP dependence plot for variable 2
    变量2的SHAP依赖图

    与替代方法的比较 (Comparison to alternative approaches)

    In the following, we compare the GPBoost algorithm to several existing approaches using the above simulated data. We consider the following alternative approaches:

    接下来,我们使用上述模拟数据将GPBoost算法与几种现有方法进行比较。 我们考虑以下替代方法:

    • A linear mixed effects model (‘Linear_ME’) where F(X) is a linear function

      线性混合效果模型('Linear_ME') ,其中F(X)是线性函数

    • Standard gradient tree-boosting ignoring the grouping structure (‘Boosting_Ign’)

      标准梯度树增强忽略分组结构('Boosting_Ign')

    • Standard gradient tree-boosting including the grouping variable as a categorical variables (‘Boosting_Cat’)

      标准梯度树增强功能,包括将分组变量作为分类变量('Boosting_Cat')

    • Mixed-effects random forest (‘MERF’) (see here and Hajjem et al. (2014) for more information)

      混合效应随机森林(“ MERF”) (有关更多信息,请参见此处Hajjem等(2014) )

    We compare the algorithms in terms of predictive accuracy measured using the root mean square error (RMSE) and computational time (clock time in seconds). The results are shown in the table below. The code for producing these results can be found below in the appendix.

    我们根据均方根误差(RMSE)和计算时间(以秒为单位的时钟时间)测得的预测准确性比较算法。 结果如下表所示。 产生这些结果的代码可以在下面的附录中找到。

    Image for post
    Comparison of GPBoost and alternative approaches.
    GPBoost与替代方法的比较。

    We see that GPBoost and MERF perform clearly best (and almost equally well) in terms of predictive accuracy. Further, the GPBoost algorithm is approximately 1000 times faster than the MERF algorithm. The linear mixed effects model (‘Linear_ME’) and tree-boosting ignoring the grouping variable (‘Boosting_Ign’) have clearly lower predictive accuracy. Tree-boosting with the grouping variable included as a categorical variable also shows lower predictive accuracy than GPBoost or MERF.

    我们看到,就预测准确性而言,GPBoost和MERF的表现明显最佳(并且几乎同样出色)。 此外,GPBoost算法比MERF算法快约1000倍。 线性混合效果模型('Linear_ME')和忽略分组变量的'boosting'(boosting_Ign')具有明显较低的预测准确性。 与GPBoost或MERF相比,将分组变量作为类别变量包括在内的树式提升也显示出较低的预测准确性。

    Note that, for simplicity, we do only one simulation run (see Sigrist (2020) for a much more detailed comparison). Except for MERF, all computations are done using the GPBoost library version 0.2.1 compiled with MSVC version 19.24.28315.0. Further, we use the MERF Python package version 0.3.

    请注意,为简单起见,我们仅进行一次模拟运行(有关 详细比较, 请参阅 Sigrist(2020) )。 除MERF外,所有计算均使用GPBoost库0.2.1版和MSVC 19.24.28315.0版进行编译。 此外,我们使用MERF Python软件包0.3版。

    结论 (Conclusions)

    GPBoost allows for combining mixed effects models and tree-boosting. If you apply linear mixed effects models, you should investigate whether the linearity assumption is indeed appropriate. The GPBoost model allows for relaxing this assumption. It may help you to find non-linearities and interactions and achieve higher predictive accuracy. If you are a frequent user of boosting algorithms such as XGBoost and LightGBM and you have categorical variables with potentially high-cardinality, GPBoost (which extends LightGBM) can make learning more efficient and result in higher predictive accuracy.

    GPBoost允许将混合效果模型和树加速结合在一起。 如果应用线性混合效应模型,则应调查线性假设是否确实合适。 GPBoost模型允许放宽此假设。 它可以帮助您发现非线性和相互作用,并获得更高的预测准确性。 如果您经常使用诸如XGBoost和LightGBM之类的增强算法,并且您的分类变量具有潜在的高基数,那么GPBoost(扩展了LightGBM)可以使学习 效率更高, 并获得更高的预测准确性。

    To the best of our knowledge, the GPBoost library is currently unmatched in terms of computational speed and predictive accuracy. Additional advantages are that GPBoost supports a range of model interpretation tools (variable importance values, partial dependence plots, SHAP values etc.). Further, it also supports other types of random effects such as Gaussian processes in addition to grouped or clustered random effects.

    据我们所知,GPBoost库目前在计算速度和预测准确性方面无与伦比。 GPBoost的其他优点是支持多种模型解释工具(可变重要性值,偏相关图,SHAP值等)。 此外,除了分组或聚类的随机效应之外,它还支持其他类型的随机效应,例如高斯过程。

    Hopefully, you have found this article useful. More information on GPBoost can be found in the companion article Sigrist (2020) and on github.

    希望您发现本文很有用。 有关GPBoost的更多信息,请参见配套文章Sigrist(2020)github

    翻译自: https://towardsdatascience.com/tree-boosted-mixed-effects-models-4df610b624cb

    面板模型混合效应模型

    展开全文
  • Stata面板数据做回归分析时: 用固定效应模型(xtreg,fe)时,五个解释变量有四个不显著 用多元回归,将年份做虚拟变量(reg i.year,r),五个解释变量只有两个不显著 请问这两个模型的区别是什么呢?哪个更适用呢
  • 多重中介的调节模型的检验及实证分析,彭月,吴海英,在效应分析模型中,有时存在多个中介变量,因此推广了一种新模型---有多重中介的调节模型,该模型分为有多个并列中介的调节模型
  • R| 混合效应模型,lme4

    万次阅读 2019-04-09 23:45:16
    ########################### 混合效应模型 ##################################### library(lme4) politeness=read.csv("/Users/Gavin/Downloads/politeness_data.csv") # head(politeness), tail(politeness),...
  • 单因素试验固定效应模型方差分析

    千次阅读 2020-07-01 20:58:11
    单因素试验固定效应模型方差分析 观测值的线性模型 平方和与自由度分解 例题与SPSS求解 非平衡单因素试验SPSS求解 一、观测值的线性模型 单因素试验线性可加模型为: Yij为第i个处理的第j个观测值;U为所有观测值...
  • 文本主要阐述了R语言中如何用plm包做面板数据分析,几种模型的具体实现。 包括pool模型,LSDV模型,组内模型(within),时间固定效应Panel模型 和个体和时间双维固定效应模型。 提供了一个数据集供读者以复现
  • 多重中介模型又可按照中介变量之间是否存在顺序关系,分为并行多重中介模型(parallel multiple mediation,图2)以及链式多重中介模型(serial multiple mediation)。 并行多重中介分析 并行多重中介分析表示X→Y
  • 这期推送简单谈一下我本人对固定效应与交互固定效应一些或许不太成熟的理解。 Note:Note:Note: 1、该文首发于微信公众号DMETP,欢迎关注;2、需要本次推送所使用的数据和代码的朋友,可以在公众号后台对话框内回复...
  • 今天给大家写广义混合效应模型Generalised Linear Random Intercept Model的第一部分 ,混合效应logistics回归模型,这个和线性混合效应模型一样也有好几个叫法: Mixed Effects Logistic Regression is sometimes...
  • 结构方程模型-中介效应检验(Amos)

    万次阅读 多人点赞 2020-04-23 20:55:29
    一、中介效应含义 考虑自变量X对因变量Y的影响,如果X通过影响变量M、N等其它变量而对Y产生影响,则称M、N等为中介变量。 下图展示了X通过M最终到Y的过程,a表示X到M的系数,b表示M到Y的系数,c表示X到Y的总效果...
  • 【综合应用】中介效应PLS-SEM模型STATA实战
  • 线性模型是一类统计模型的总称,它包括了线性回归模型、方差分析模型、协方差分析模型和线性混合效应模型(或称方差分量模型)等。许多生物、医学、经济、管理、地质、气象、农业、工业、工程技术等领域的现象都可以...
  • 基于多重运移机制的页岩气藏有限导流压裂井试井模型,王海涛,郭晶晶,页岩储层中烃类运移受多重机制的作用:既有渗流,也有解吸和扩散;此外,应力敏感效应在页岩储层中也很明显,而各向异性也是储层
  • 利用不动点指数理论给出了一类带有加法Allee效应的捕食-食饵模型正解存在的充分条件,讨论了正解的惟一性和稳定性,运用扰动理论研究了参数b2充分大时正解的多重性。结果表明在一定条件下系统存在多重解和惟一解。
  •  为了解决因页岩气在基质-缝网系统中的多重输运机制和多尺度流动效应导致的数学描述复杂、耦合求解困难等问题,综合考虑页岩气解吸附、扩散、渗流等特征因素,基于系列实验创建的页岩气高压等温吸附模型、广义渗透率...
  • 在方差分量模型中,把既含有固定效应,又含有随机效应模型,称为混合线性模型【信息来源:百度】一般线性模型中仅包含固定效应和噪声两项影响因素,也就是不会考虑随机因素对模型的影响。数据集:R的MASS包中oast...
  • 多重线性回归模型的最终建立不仅仅拟合个方程就完事了,还需要进行适用条件的考察、模型的诊断以及改进模型的再评估等。 (1)拟合多重线性回归模型; (2)适用条件考察:线性、独立性、正态性、同方差性; (3...
  • 这期推送将比较时间固定效应和时间趋势项的区别,并使用两种方法对模型中可能存在的trend进行识别。 Note:Note:Note: 1、该文首发于微信公众号DMETP,欢迎关注;2、需要本次推送所使用的数据和代码的朋友,可以在...
  • 一般线性模型和混合线性模型 生命科学的数学统计和机器学习 (Mathematical Statistics and Machine Learning for Life Sciences) This is the seventeenth article from my column Mathematical Statistics and ...
  • 具体而言,存在自选择偏差的回归方程中被遗漏的IMR计算公式为: 明显可以看到,在处理效应模型中,D取值为1的样本与D取值为0的样本的 IMR 计算公式不同,而且由于处理效应模型第二阶段回归中所有样本均参与了回归,...
  • 为了研究围压对煤岩强度和破坏特征的影响,对围压效应进行了理论分析,并利用FLAC2D模拟不同围压下煤岩的破坏情况。研究表明:煤岩极限破坏强度随围压的增大而增加,二者之间呈非线性关系;煤岩的轴向应力也随围压的增大...
  • 请问实证论文DID模型可以不控制年度和行业固定效应吗?即y=treat+post+treat*post+控制变量,但看大多数文章都控制了,但不是有多重共线性问题嘛
  • 您好,请问多重中介效应模型怎么看每个变量对y的影响呢?想知道哪一个变量在传导中发挥了更大的作用。

空空如也

空空如也

1 2 3 4 5 ... 20
收藏数 4,553
精华内容 1,821
关键字:

多重效应模型