👇 连享会 · 推文导航 |
www.lianxh.cn
🍓 课程推荐:
连享会:2025 寒假班
嘉宾:连玉君(初级|高级);杨海生(前沿)
时间:2025 年 1 月 13-24 日
咨询:王老师 18903405450(微信)
温馨提示:
文中链接在微信中无法生效。请点击底部
「阅读原文」
。或直接长按/扫描如下二维码,直达原文:
译者
:陈佳慧 (北京理工大学)
邮箱
:[email protected]
编者按:
本文主要摘译自以下文章,特此感谢!
Fernando Rios-Avila. Blog. Story of marginal effects: The IV-Probit:
margins
and
ml
. -Link-
Rios-Avila, F., & Canavire-Bacarreza, G. (2018). Standard-error Correction in Two-stage Optimization Models: A Quasi–maximum Likelihood Estimation Approach. The Stata Journal, 18(1), 206–222. -Link- -PDF- -Google-
1. 简介
在之前的文章中,作者已经展示了如何在
ml
之后使用
margins
命令进行线性回归模型 (假设正态) 和 probit 模型的估计。今天,作者将提供一个更复杂的例子:双方程模型的估计,特别是 "IVprobit" (工具变量 probit)。
这个模型因其不一致的边际效应估计而声名狼藉,在 Statalist 上引发了多个讨论。其问题源于不同版本的 Stata 中
margins
命令产生的不同版本的“边际效应”。从计算角度来看,每个版本都是正确的,但并不总是符合直观理解。
作者的观点如下:
Stata 13 估计了 IVprobit MLE 的正确边际效应,但未能为两步法提供正确的边际效应。
Stata 14 和 Stata 15 估计的是全信息边际效应,尽管在技术上是正确的,但与常识相悖。Wooldridge 教授对此问题进行了广泛讨论,并建议采用 MLE 或两步边际效应法。
最近,Stata 16 发布了一个更新,现在能够生成等同于 MLE 估计的两步边际效应,这很可能是近期讨论的结果。
接下来将介绍作者版本的 IVprobit 模型估计方法以及预测程序的设置方法,并解释为何使用多种类型的边际效应。
2. 设置
首先,请确保内存中包含以下程序。除非你已经将该程序保存到了可访问的 ado 文件中 (最有可能的是 "ado/personal" 文件夹)。该程序将允许你添加或修改
e()
中的信息,所有估算命令都会在
e()
中存储信息。
clear all program adde, eclass ereturn `0' end
IVprobit 模型是一种用于二元因变量 (0-1) 的非线性模型,特别适用于存在一个或多个控制因素内生性的情况。在这种情况下,模型用于估计在一系列特征条件下的“成功概率 (y=1)”,但我们关注的特征是一个连续但内生的解释变量。
当一个变量是内生变量时,我们无法将模型的估计结果直接解释为因果效应。这是因为内生变量的变化可能与未观测变量的变化同时发生。因此,如果结果发生变化,我们无法确定是由于内生变量的变化还是未观测变量的变化所引起的 (毕竟,它们是相互关联的)。
要处理未观测到的混杂因素,有两种选择:
可以使用工具分离相关变量的外生变化 (例如使用 2SLS 方法)。
使用工具来获得可以直接控制的内生成分的近似值 (控制函数法)。
事实上,IV-probit 就是后者的应用:控制函数法。从形式上看,IV-probit 模型可以写成如下形式:
其中,误差
,
服从二元正态分布:
在此模型中,
和
是外生变量,其中
是一组工具变量,
则是模型中连续但内生的变量。我们无法直接观察到潜在变量
,而只能观察到
,它仅取 0 或 1 的值。那么我们如何估计这个模型呢?让我们一步步来进行分析。
方程 (1) 可以直接估计,因为它仅仅是外生变量的函数。因此,我们可以使用标准的 OLS 方法 (如 ivprobit 两步法所采用的) 或在假设误差项正态性的前提下,通过 MLE 方法来估计该方程。
更值得注意的是方程 (2)。由于
不等于 0,
因此具有内生性。不过,
可以分解为两个部分:一个包含内生成分,另一个是外生的,与所有其他变量无关。
要实现这一点,如果
和
遵循二元正态分布,那么在给定
的条件下,
将具有以下分布:
这意味着,
写成以下形式:
在这种情况下,通过构造,
与
或
不相关。将式 (4) 代入式 (1),可得:
假设我们可以观察到
,则这个方程可以直接估计。然而,为了使用 probit 模型进行估计,还需要对方程进行重新标度,以确保重新标度后的方程中误差项
的方差为 1。
根据方程 (3),我们知道
的方差为
。因此,需将方程 (5) 中的所有项除以
,然后使用标准 probit 模型估计以下模型 (或其简化模型)。
只要知道如何从系统中估计标准误差,使用这两种方程没有区别。
3. 实际估算
用于估计 ivprobit 模型的方法至少有三种。
第一种方法
是两步法。该方法使用 OLS 对方程 (1) 进行估计,从而得到预测残差,然后将其代入方程 (7)。这是通过一个简单的概率模型来估计的。
这种方法存在两个问题:
它仅能提供“重新调整”系数的估计值,而非结构系数。
它无法提供正确的标准误差估计,因为没有考虑到残差带来的第一步误差。
一些教科书建议,解决这一问题可以简单应用 delta 方法或使用 bootstrap。然而,实际上需要考虑到方程 (1) 的残差是估计的,而不是真实的残差
。
第二种方法
是使用全信息最大似然法,同时估计方程 (1) 和方程 (5)。该方法假设误差服从二元正态分布,因此不仅可以获得结构参数的估计值,还可以获得“链接”参数
和
,并提供正确的标准误差。在这种策略下,单个观测值对似然函数的贡献变为:
请注意,这里并未在 probit 方程中插入
,而是明确地添加了
。这样做是为了解释第一阶段中的测量误差。
还有
第三种选择
,作者称之为
Two-Step-MLE
。之所以这样命名,是因为 IVprobit 将使用公式 (1) 和公式 (7) 进行估计。称其为 MLE,是因为两个方程同时使用最大似然估计 (MLE) 进行估计:
与标准 FIML 的不同之处在于,该方法仅估计重新标度的系数,并且两个方程之间的联系为
而不是
。尽管如此,在这个简化的例子中,两个方程识别的模型是完全相同的。如果您对 -ml- 的这种用法感兴趣,请参阅作者的论文 Rios-Avila 和 Canavire-Bacarreza (2018)。
然而,与通常的两步方法相比,由于模型是同时估计的,因此所有系数的标准误差都是正确估计的,而无需进一步计算 (没有 delta 方法也没有 bootstrap)。
在该模型中还有最后一点需要注意。
θ
和
之间存在密切关系,这将影响到重新标定的参数:
4. 对数似然函数
虽然已经证明使用 FIML 和两步 ML 方法会提供相同的结果,但作者坚持使用两步法,因为它允许在不同的 Stata 版本中推导边际效应,从而更清楚地说明边际变化的情况。
下面的程序使用两步方法 (公式 8) 定义了 IVprobit 的对数似然函数,并使用了以下规格说明:
program myivprobit_2sls args lnf xb theta zb lnsigma qui { local y1 $ML_y1 local y2 $ML_y2 local u2 (`y2'-`zb') tempvar xb_zb p1 p0 gen double `xb_zb'= `xb'+`theta'*((`u2')/exp(`lnsigma')) gen double `p1' = normal( `xb_zb') gen double `p0' = normal(-`xb_zb') tempvar lnf1 lnf2 gen double `lnf1' = log(normalden(`y2', `zb', exp(`lnsigma'))) gen double `lnf2' = log(`p1') if `y1'==1 replace `lnf2' = log(`p0') if `y1'==0 replace `lnf' = `lnf1' + `lnf2' } end
5. 预测程序
最后,我们来讨论一个颇具争议的部分:成功概率的预测!
之所以存在争议,是因为有两个候选方案可以确定这个表达式。
第一个候选方案与结构方程 2 相关。基本上,如果我们能够估算出无标度系数,那么预测结果就可以通过以下方式确定:
或者,如果更喜欢基于重新标定系数的版本:
因此,边际效应可以通过单独分析其中一个方程来获得。只有当我们使用 FIML 或使用重新标定的系数估计结构方程时,才能直接确定该表达式的标准误差,并确保在计算标准误差时考虑到第一阶段的估计误差。
第二种方法是使用方程 (7) 估计边际效应:
但是,因为
从未观察到,通常建议平均 (但不忽略)
在等式上的影响:
底线:如果使用两步法,边际效应的估计可以将
视为模型中的另一个外生变量,关键在于如何正确估计标准误差。
因此,建议将这两个选项编入“预测”程序中。
program myivprobit_p syntax newvarname [if] [in] , [ pr1 pr2 *] if "`pr1'`pr2'" =="" { ml_p `0' } tokenize `e(depvar)' local y1 `1' local y2 `2' marksample touse, novarlist if "`pr1'" !="" { tempvar xb zb theta lnsigma _predict double `xb' , eq(#1) _predict double `theta', eq(#2) _predict double `zb' , eq(#3) _predict double `lnsigma', eq(#4) gen `typlist' `varlist' = /// normal(`xb'+`theta'*(`y2'-`zb')/exp(`lnsigma')) if `touse' label var `varlist' "P(y=1|X) two-step" } else if "`pr2'"!="" { tempvar xb zb theta lnsigma _predict double `xb' , eq(#1) _predict double `theta' , eq(#2) gen `typlist' `varlist' = /// normal(`xb'/sqrt(1+`theta'^2)) if `touse' label var `varlist' "P(y=1|X) FIML" } end
第一个选项
pr1
将估计预测概率,就像使用两步法估计模型一样,而第二个选项将根据结构方程估计预测概率。
下面为估计模型并将结果与内置的
ivprobit
命令进行比较:
. webuse laborsup, clear . global y1 fem_work . global z1 fem_educ kids . global y2 other_inc . global z2 male_educ . *Built in command: . ivprobit $y1 $z1 ($y2 = $z2), two Two-step probit with endogenous regressors Number of obs = 500 Wald chi2(3) = 93.97 Prob > chi2 = 0.0000 ------------------------------------------------------------------------------ | Coefficient Std. err. z P>|z| [95% conf. interval] -------------+---------------------------------------------------------------- other_inc | -0.058 0.009 -6.26 0.000 -0.077 -0.040 fem_educ | 0.227 0.028 8.08 0.000 0.172 0.283 kids | -0.196 0.050 -3.95 0.000 -0.293 -0.099 _cons | 0.396 0.498 0.79 0.427 -0.581 1.372 ------------------------------------------------------------------------------ Wald test of exogeneity: chi2(1) = 6.50 Prob > chi2 = 0.0108 Endogenous: other_inc Exogenous: fem_educ kids male_educ . *my ivprobit two-step . ml model lf myivprobit_2sls ($y1 = $z1 $y2 ) /// > (theta:) ($y2 = $z1 $z2 ) (lnsigma:) , /// > technique(nr bhhh) init(lnsigma:_cons = 2.81 ) maximize nolog . ml display Number of obs = 500 Wald chi2(3) = 94.01 Log likelihood = -2368.2062 Prob > chi2 = 0.0000 ------------------------------------------------------------------------------ | Coefficient Std. err. z P>|z| [95% conf. interval] -------------+---------------------------------------------------------------- eq1 | fem_educ | 0.227 0.028 8.08 0.000 0.172 0.283 kids | -0.196 0.050 -3.95 0.000 -0.293 -0.099 other_inc | -0.058 0.009 -6.26 0.000 -0.077 -0.040 _cons | 0.396 0.498 0.79 0.427 -0.581 1.372 -------------+---------------------------------------------------------------- theta | _cons | 0.401 0.163 2.46 0.014 0.082 0.720 -------------+---------------------------------------------------------------- eq3 | fem_educ | 0.335 0.283 1.19 0.236 -0.219 0.889 kids | 0.833 0.548 1.52 0.128 -0.240 1.906 male_educ | 2.845 0.283 10.06 0.000 2.291 3.399 _cons | 9.873 5.029 1.96 0.050 0.016 19.730 -------------+---------------------------------------------------------------- lnsigma | _cons | 2.813 0.032 88.97 0.000 2.751 2.875 ------------------------------------------------------------------------------
可以看到,除了因四舍五入误差和自由度造成的差异外,其他结果几乎相同。同样令人欣慰的是,当比较 ivprobit-mle 和重标系数时,结果也是一样的:
. *FIML . ivprobit $y1 $z1 ($y2 = $z2), ml nolog Probit model with endogenous regressors Number of obs = 500 Wald chi2(3) = 163.88 Log likelihood = -2368.2062 Prob > chi2 = 0.0000 ---------------------------------------------------------------------------------------------- | Coefficient Std. err. z P>|z| [95% conf. interval] -----------------------------+---------------------------------------------------------------- other_inc | -0.054 0.006 -8.92 0.000 -0.066 -0.042 fem_educ | 0.211 0.027 7.86 0.000 0.158 0.264 kids | -0.182 0.048 -3.81 0.000 -0.276 -0.088 _cons | 0.367 0.448 0.82 0.412 -0.511 1.245 -----------------------------+---------------------------------------------------------------- corr(e.other_inc,e.fem_work)| 0.372 0.130 0.095 0.596 sd(e.other_inc)| 16.666 0.527 15.665 17.732 ---------------------------------------------------------------------------------------------- Wald test of exogeneity (corr = 0): chi2(1) = 6.70 Prob > chi2 = 0.0096 Endogenous: other_inc Exogenous: fem_educ kids male_educ
. *my ivprobit two-step . ml model lf myivprobit_2sls ($y1 = $z1 $y2 ) (theta:) ($y2 = $z1 $z2 ) (lnsigma:), /// > technique(nr bhhh) init(lnsigma:_cons = 2.81 ) maximize nolog . adde local predict myivprobit_p . est store myivp . *with rescaled coefficients: . nlcom (other_inc: _b[other_inc]/sqrt(1+_b[theta:_cons]^2)) /// > (fem_educ: _b[fem_educ]/sqrt(1+_b[theta:_cons]^2)) /// > (kids: _b[kids]/sqrt(1+_b[theta:_cons]^2)) /// > (cons: _b[_cons]/sqrt(1+_b[theta:_cons]^2)) other_inc: _b[other_inc]/sqrt(1+_b[theta:_cons]^2) fem_educ: _b[fem_educ]/sqrt(1+_b[theta:_cons]^2) kids: _b[kids]/sqrt(1+_b[theta:_cons]^2) cons: _b[_cons]/sqrt(1+_b[theta:_cons]^2) ------------------------------------------------------------------------------ | Coefficient Std. err. z P>|z| [95% conf. interval] -------------+---------------------------------------------------------------- other_inc | -0.054 0.006 -8.92 0.000 -0.066 -0.042 fem_educ | 0.211 0.027 7.86 0.000 0.158 0.264 kids | -0.182 0.048 -3.81 0.000 -0.276 -0.088 cons | 0.367 0.448 0.82 0.412 -0.511 1.245 ------------------------------------------------------------------------------
再次,显示完全相同的结果。
6. 边际效应的故事
现在将带大家了解 IVprobit 边际效应的故事。
6.1 Stata 13
在 Stata 13 中,IV 概率的边际效应是通过使用结构方程系数进行估计的:
因此,边际效应被定义为:
注意:这里作者是使用 Quarto 重写他的一些旧文章。因此,他只能动态使用 Stata 17。因此,除非你拥有相同版本的 Stata,否则你在下面看到的代码将无法重现。
如果您有权访问 Stata 13,您将能够重现以下输出:
. margins, dydx(*) predict(pr) Average marginal effects Number of obs = 500 Model VCE : OIM Expression : Probability of positive outcome, predict(pr) dy/dx w.r.t. : other_inc fem_educ kids male_educ ------------------------------------------------------------------------------ | Delta-method | dy/dx Std. Err. z P>|z| [95% Conf. Interval] -------------+---------------------------------------------------------------- other_inc | -.014015 .0009836 -14.25 0.000 -.0159428 -.0120872 fem_educ | .0545129 .0066007 8.26 0.000 .0415758 .06745 kids | -.0470199 .0123397 -3.81 0.000 -.0712052 -.0228346 male_educ | 0 (omitted) ------------------------------------------------------------------------------
这种边际效应在 myivprobit 之后使用
pr2
模拟:
. est restore myivp . margins, dydx(*) predict(pr2) force Average marginal effects Number of obs = 500 Model VCE: OIM Expression: P(y=1|X) FIML, predict(pr2) dy/dx wrt: fem_educ kids other_inc male_educ ------------------------------------------------------------------------------ | Delta-method | dy/dx std. err. z P>|z| [95% conf. interval] -------------+---------------------------------------------------------------- fem_educ | 0.055 0.007 8.26 0.000 0.042 0.067 kids | -0.047 0.012 -3.81 0.000 -0.071 -0.023 other_inc | -0.014 0.001 -14.25 0.000 -0.016 -0.012 male_educ | 0.000 (omitted) ------------------------------------------------------------------------------
可以看到,
margins
的输出将 male_educ 列入了外生变量列表。这是因为它在模型中至少是一个方程 (第一个方程) 的解释变量。然而,由于在第二个方程中没有包含这个变量,因此它对第二个方程没有影响。
6.2 Stata 14.1
当使用 Stata 14.1 时,IVprobit 后的概率计算方法发生了变化。正如 "whatsnew" 材料中所说,新的计算方法将考虑内生性。
具体来说,他们使用的是作者所说的 2sls 预测概率,即公式 (5),但需要注意的是
被替换为
:
上述两个等式基本相同,但在用软件估算边际效应时,它们却有重要区别。具体来说,成功的概率现在是
的函数。
那么 Stata 14.1 做了什么,为了估计边际效应,偏导数基于公式 (11)。
从技术角度来看,这些偏导数是正确的,因为它们同时捕捉到了所有变量对成功概率的直接和间接影响。类似于总导数,而不是偏导数。
然而,问题在于,这假定我们可以实际观察到当其他变量发生变化时,未观察到的部分是如何变化的。然而,标准的回归分析认为,这些未观察到的因素应该被视为固定的,而我们应该对未观察到的因素的边际效应进行平均估算。因此,上述每个导数的第二项都应为零。
尽管如此,如果用 Stata 14.2 估算边际效应,会得到以下结果:
. margins, dydx(*) predict(pr) Average marginal effects Number of obs = 500 Model VCE : OIM Expression : Probability of positive outcome, predict(pr) dy/dx w.r.t. : other_inc fem_educ kids male_educ ------------------------------------------------------------------------------ | Delta-method | dy/dx Std. Err. z P>|z| [95% Conf. Interval] -------------+---------------------------------------------------------------- other_inc | -.0097802 .0014994 -6.52 0.000 -.012719 -.0068414 fem_educ | .0623273 .007099 8.78 0.000 .0484135 .076241 kids | -.0614265 .0139446 -4.41 0.000 -.0887574 -.0340956 male_educ | -.0194406 .0022103 -8.80 0.000 -.0237728 -.0151084 ------------------------------------------------------------------------------
为了使用
myivprobit
复制这一点,作者使用选项 pr1 来估算边际效应,并要求在不使用链式规则 (nochain) 的情况下估计导数。这样可以确保考虑到以下所有变量变化的影响:预测结果
中的
,
和
。
. est restore myivp . margins, dydx(*) predict(pr1) force nochain Average marginal effects Number of obs = 500 Model VCE: OIM Expression: P(y=1|X) two-step, predict(pr1) dy/dx wrt: fem_educ kids other_inc male_educ ------------------------------------------------------------------------------ | Delta-method | dy/dx std. err. z P>|z| [95% conf. interval] -------------+---------------------------------------------------------------- fem_educ | 0.062 0.006 10.33 0.000 0.051 0.074 kids | -0.061 0.013 -4.78 0.000 -0.087 -0.036 other_inc | -0.010 0.001 -9.95 0.000 -0.012 -0.008 male_educ | -0.019 0.007 -2.60 0.009 -0.034 -0.005 ------------------------------------------------------------------------------
6.3 Stata 16
早期版本的 Stata 16 存在与上述问题非常相似的情况。不过,由于本文早期版本的某些原因,Stata 于 2020 年 11 月对
ivprobit
及其他相关命令的边际效应估计方法进行了修正 (参见 2020 年 11 月 19 日的更新)。