Stata连享会 主页 || 视频 || 推文 || 知乎 || Bilibili 站
温馨提示: 定期 清理浏览器缓存,可以获得最佳浏览体验。
New!
lianxh
命令发布了:
随时搜索推文、Stata 资源。安装:
. ssc install lianxh
详情参见帮助文件 (有惊喜):
. help lianxh
连享会新命令:cnssc
,ihelp
,rdbalance
,gitee
,installpkg
⛳ Stata 系列推文:
作者:陈卓然 (中山大学)
邮箱:chenzhr25@mail2.sysu.edu.cn
目录
因果效应识别的最大挑战之一就是遗漏变量。在线性模型中,当一个控制变量与遗漏变量相关时,这个控制变量的系数将是有偏的。为看清这一点,不妨回忆一下最简单的二元线性回归模型。假定真实模型如下:
从而正确的估计结果为:
然而由于某种原因,我们实际估计的模型如下:
从而错误的估计结果为:
我们将式 (1) 和式 (2) 写成离差的形式:
将式 (3) 左右两边同时乘以
简单整理之后可得:
因此,
其中,
因此遗漏变量偏误的表达式为:
等式 (4) 的推导利用到
其中,
这其中有几点值得注意,如果
当然对于上述遗漏变量偏误的推导也有一些别的方法,比如:
不过上述的推导过程尽管看似复杂,但是很容易推广到更多元的情形之下,因此我们在这里将其较为详细地列示出来,当然如果读者掌握更高级的矩阵知识的话,也可以有更为简便的办法,在这里就不再赘述了。
由于在实证分析中,我们很难将全部的解释变量都纳入进模型当中,因此将不可避免地受到遗漏变量的影响。尽管如此,我们仍可以分析结果对于遗漏变量偏误的敏感性。
一些经典的研究遗漏变量对于实证结果的影响的文献 (Altoji,2005;Oster,2019) 均要求外生控制变量假设,即遗漏变量必须与控制变量不相关,而这一假设通常是一个很强而且不太合理的假设。本文主要探讨在 Diegert 等 (2022) 提出的敏感性分析方法,该方法允许内生控制变量的存在。在这一框架下,有两个核心问题:
接下来,我们将使用命令 regsensitivity
对上述两种敏感性进行分析。
这里采用 Bazzi 等 (2020) 的数据进行演示,我们首先来复现 Bazzi 等 (2020) 中表 3 第 7 列。
. lxhget regsensitivity.pkg, install
. lxhuse bfg2020.dta, clear
. local y avgrep2000to2016
. local x tye_tfe890_500kNI_100_l6
. local w1 log_area_2010 lat lon temp_mean rain_mean elev_mean d_coa d_riv d_lak ave_gyi
. local w0 i.statea
. local w `w1' `w0'
. local SE cluster(km_grid_cel_code)
. reg `y' `x' `w', `SE'
Linear regression Number of obs = 2,036
F(39, 379) = .
Prob > F = .
R-squared = 0.3321
Root MSE = 9.6368
(Std. err. adjusted for 380 clusters in km_grid_cel_code)
------------------------------------------------------------------------------------------
| Robust
avgrep2000to2016 | Coefficient std. err. t P>|t| [95% conf. interval]
-------------------------+----------------------------------------------------------------
tye_tfe890_500kNI_100_l6 | 2.055 0.349 5.88 0.000 1.368 2.741
log_area_2010 | 0.276 0.980 0.28 0.778 -1.651 2.203
lat | 2.265 1.101 2.06 0.040 0.100 4.430
lon | 0.011 0.291 0.04 0.970 -0.562 0.584
temp_mean | 1.627 1.068 1.52 0.128 -0.473 3.728
rain_mean | 0.016 0.005 3.58 0.000 0.007 0.026
elev_mean | 0.015 0.004 4.10 0.000 0.008 0.023
d_coa | 0.000 0.000 2.62 0.009 0.000 0.000
d_riv | 0.000 0.000 3.10 0.002 0.000 0.000
d_lak | 0.000 0.000 0.07 0.945 -0.000 0.000
ave_gyi | -3.780 10.810 -0.35 0.727 -25.035 17.475
statea |
5 | -4.214 3.386 -1.24 0.214 -10.872 2.445
8 | -27.317 6.247 -4.37 0.000 -39.600 -15.034
12 | 4.628 3.355 1.38 0.169 -1.968 11.224
13 | 0.540 2.644 0.20 0.838 -4.658 5.738
17 | -10.258 3.787 -2.71 0.007 -17.705 -2.811
18 | -5.924 3.497 -1.69 0.091 -12.801 0.952
19 | -18.027 4.514 -3.99 0.000 -26.903 -9.151
20 | 1.599 4.634 0.35 0.730 -7.512 10.710
21 | -0.504 3.186 -0.16 0.874 -6.768 5.760
22 | 0.894 3.277 0.27 0.785 -5.549 7.337
26 | -14.063 4.406 -3.19 0.002 -22.725 -5.401
27 | -18.103 4.821 -3.75 0.000 -27.583 -8.623
28 | -6.931 3.676 -1.89 0.060 -14.159 0.297
29 | -4.170 3.902 -1.07 0.286 -11.843 3.502
31 | -1.343 4.738 -0.28 0.777 -10.658 7.972
35 | -40.780 9.264 -4.40 0.000 -58.996 -22.564
36 | -9.822 4.689 -2.09 0.037 -19.041 -0.602
37 | -13.538 4.242 -3.19 0.002 -21.878 -5.197
38 | -11.982 5.512 -2.17 0.030 -22.821 -1.143
39 | -6.191 3.656 -1.69 0.091 -13.378 0.997
40 | 13.600 4.881 2.79 0.006 4.004 23.197
42 | -3.146 4.426 -0.71 0.478 -11.850 5.557
46 | -11.847 5.102 -2.32 0.021 -21.878 -1.816
47 | -3.541 2.794 -1.27 0.206 -9.035 1.953
48 | 12.826 4.174 3.07 0.002 4.619 21.033
51 | -0.812 4.048 -0.20 0.841 -8.771 7.147
54 | -3.244 3.570 -0.91 0.364 -10.264 3.776
55 | -18.929 4.504 -4.20 0.000 -27.785 -10.073
56 | -19.023 9.729 -1.96 0.051 -38.153 0.107
_cons | -73.535 57.847 -1.27 0.204 -187.277 40.206
------------------------------------------------------------------------------------------
接着,进行默认的敏感性分析。
. regsensitivity `y' `x' `w', compare(`w1')
Analysis : DMP (2022) Number of obs = 2,036
Beta(short) = 1.925
Treatment : tye_tfe890_500kNI_100_l6 Beta(medium) = 2.055
Outcome : avgrep2000to2016 R2(short) = 0.033
R2(medium) = 0.105
Var(Y) = 101.739
Var(X) = 0.901
Var(X_Residual) = 0.882
Hypothesis : Beta > 0 Breakdown point = 80.4%
Other Params : cbar = 1, rybar = +inf
--------------------------------------------------------------------------------
rxbar Beta
--------------------------------------------------------------------------------
0.000 [ 2.05, 2.05 ]
0.095 [ 1.91, 2.20 ]
0.196 [ 1.76, 2.35 ]
0.296 [ 1.59, 2.52 ]
0.397 [ 1.41, 2.70 ]
0.497 [ 1.20, 2.91 ]
0.592 [ 0.95, 3.16 ]
0.693 [ 0.61, 3.50 ]
0.793 [ 0.07, 4.04 ]
0.894 [ -1.05, 5.16 ]
0.989 [-58.90, 63.01 ]
0.989 [ -inf, +inf ]
--------------------------------------------------------------------------------
事实上,上述结果已经回答了我们在引言中提到的两个问题:
考虑如下模型:
其中,
符号 | 含义 |
---|---|
所谓可识别集合是指在现有假设下给定可观测数据的分布,regsensitivity bounds
的子命令完成。
regsensitivity bounds depvar indepvar controls, options...
其中,depvar
是被解释变量 indepvar controls
是解释变量 regsensitivity
中,解释变量的顺序至关重要。第一个解释变量
默认情况下,regsensitivity bounds
计算在 cbar
的选项,例如:
. regsensitivity bounds `y' `x' `w', compare(`w1') cbar(.1)
Analysis : DMP (2022) Number of obs = 2,036
Beta(short) = 1.925
Treatment : tye_tfe890_500kNI_100_l6 Beta(medium) = 2.055
Outcome : avgrep2000to2016 R2(short) = 0.033
R2(medium) = 0.105
Var(Y) = 101.739
Var(X) = 0.901
Var(X_Residual) = 0.882
Hypothesis : Beta > 0 Breakdown point = 119%
Other Params : cbar = .1, rybar = +inf
--------------------------------------------------------------------------------
rxbar Beta
--------------------------------------------------------------------------------
0.000 [ 2.05, 2.05 ]
0.394 [ 1.45, 2.66 ]
0.808 [ 0.74, 3.37 ]
1.202 [ -0.02, 4.13 ]
1.617 [ -0.93, 5.04 ]
2.032 [ -2.02, 6.13 ]
2.425 [ -3.32, 7.43 ]
2.840 [ -5.17, 9.28 ]
3.234 [ -7.86, 11.97 ]
3.648 [-13.63, 17.74 ]
4.042 [-75.24, 79.35 ]
4.063 [ -inf, +inf ]
--------------------------------------------------------------------------------
为将结果绘制出来,可以使用命令 regsensitivity plot
。
注意在上述命令中,我们加入了 compare(varlist)
的选项,这一选项是在指明哪些控制变量被放入
. regsensitivity bounds `y' `x' `w', cbar(.1)
Regression Sensitivity Analysis, Bounds
Analysis : DMP (2022) Number of obs = 2,036
Beta(short) = 1.708
Treatment : tye_tfe890_500kNI_100_l6 Beta(medium) = 2.055
Outcome : avgrep2000to2016 R2(short) = 0.027
R2(medium) = 0.332
Var(Y) = 136.320
Var(X) = 1.257
Var(X_Residual) = 0.882
Hypothesis : Beta > 0 Breakdown point = 29.7%
Other Params : cbar = .1, rybar = +inf
--------------------------------------------------------------------------------
rxbar Beta
--------------------------------------------------------------------------------
0.000 [ 2.05, 2.05 ]
0.128 [ 1.20, 2.91 ]
0.262 [ 0.25, 3.86 ]
0.396 [ -0.77, 4.88 ]
0.531 [ -1.91, 6.02 ]
0.665 [ -3.24, 7.35 ]
0.800 [ -4.88, 8.99 ]
0.934 [ -7.07, 11.18 ]
1.069 [ -10.45, 14.56 ]
1.203 [ -17.46, 21.57 ]
1.331 [-106.48, 110.59 ]
1.336 [ -inf, +inf ]
--------------------------------------------------------------------------------
. regsensitivity plot
不难发现,当全部的控制变量都被包含进 rxbar
值,使其在两种情形下是相同的。
由于 regsensitivity bounds
的回归结果被存储在 e(idset_table)
中,因此我们可以提取其中的 rxbar
的值,然后将其作为选项添加到从
. forvalues i=1/12{
1. local rxbar `rxbar' `= e(idset_table)[`i', 1]'
2. }
. regsensitivity bounds `y' `x' `w', compare(`w1') cbar(.1) rxbar(`rxbar')
Analysis : DMP (2022) Number of obs = 2,036
Beta(short) = 1.925
Treatment : tye_tfe890_500kNI_100_l6 Beta(medium) = 2.055
Outcome : avgrep2000to2016 R2(short) = 0.033
R2(medium) = 0.105
Var(Y) = 101.739
Var(X) = 0.901
Var(X_Residual) = 0.882
Hypothesis : Beta > 0 Breakdown point = 119%
Other Params : cbar = .1, rybar = +inf
--------------------------------------------------------------------------------
rxbar Beta
--------------------------------------------------------------------------------
0.000 [ 2.0548, 2.0548 ]
0.128 [ 1.8628, 2.2468 ]
0.262 [ 1.6550, 2.4546 ]
0.396 [ 1.4408, 2.6687 ]
0.531 [ 1.2198, 2.8898 ]
0.665 [ 0.9911, 3.1184 ]
0.800 [ 0.7540, 3.3555 ]
0.934 [ 0.5078, 3.6017 ]
1.069 [ 0.2513, 3.8583 ]
1.203 [-0.0166, 4.1262 ]
1.331 [-0.2829, 4.3924 ]
1.336 [-0.2937, 4.4032 ]
--------------------------------------------------------------------------------
不难看出,当全部的控制变量都被包括进 cbar
的选项当中。此时 regsensitivity bounds
命令将会绘制一张图,而不是以表格的形式呈现在命令窗口
. regsensitivity bounds `y' `x' `w',compare(`w1') cbar(0(.2)1)
在默认情况下,regsensitivity plot
将绘制到对于每一个
regsensitivity bounds
的结果显示了对于一个给定假设下,参数
regsensitivity
可以处理的假设形式为:对于任意的
regsensitivity bounds
显示的是去除掉州层面的固定效应的 breakdown
的子命令。例如:
. regsensitivity breakdown `y' `x' `w', compare(`w1') cbar(0(.1)1)
Analysis : DMP (2022) Number of obs = 2,036
Beta(short) = 1.925
Treatment : tye_tfe890_500kNI_100_l6 Beta(medium) = 2.055
Outome : avgrep2000to2016 R2(short) = 0.033
R2(medium) = 0.105
Var(Y) = 101.739
Var(X) = 0.901
Hypothesis : Beta > 0 Var(X_Residual) = 0.882
Other Params : rybar = +inf
--------------------------------------------------------------------------------
cbar rxbar(Breakdown)
--------------------------------------------------------------------------------
0.000 135.0%
0.100 119.5%
0.200 108.0%
0.300 99.3 %
0.400 92.7 %
0.500 87.6 %
0.600 83.9 %
0.700 81.4 %
0.800 80.4 %
0.900 80.4 %
1.000 80.4 %
--------------------------------------------------------------------------------
. regsensitivity plot
为了检验给定某个 beta(b,lb)
。其中,lb
是 lower bound 的意思。选项 beta
也可以接受一个数列以检验一系列的假设。例如,对于一系列的
. regsensitivity breakdown `y' `x' `w', compare(`w1') beta(-1(.2)1 lb)
Analysis : DMP (2022) Number of obs = 2,036
Beta(short) = 1.925
Treatment : tye_tfe890_500kNI_100_l6 Beta(medium) = 2.055
Outome : avgrep2000to2016 R2(short) = 0.033
R2(medium) = 0.105
Var(Y) = 101.739
Var(X) = 0.901
Hypothesis : Beta > Beta(Hypothesis) Var(X_Residual) = 0.882
Other Params : cbar = 1, rybar = +inf
--------------------------------------------------------------------------------
Beta(Hypothesis) rxbar(Breakdown)
--------------------------------------------------------------------------------
-1.000 89.1 %
-0.800 87.9 %
-0.600 86.5 %
-0.400 84.8 %
-0.200 82.8 %
0.000 80.4 %
0.200 77.4 %
0.400 73.8 %
0.600 69.5 %
0.800 64.1 %
1.000 57.5 %
--------------------------------------------------------------------------------
当然,我们也可以通过指定 beta(b ub)
来检验 ub
表示 upper bound。例如,我们想检验对于一系列的
. regsensitivity breakdown `y' `x' `w', compare(`w1') cbar(0(.1)1) beta(4 ub)
Analysis : DMP (2022) Number of obs = 2,036
Beta(short) = 1.925
Treatment : tye_tfe890_500kNI_100_l6 Beta(medium) = 2.055
Outome : avgrep2000to2016 R2(short) = 0.033
R2(medium) = 0.105
Var(Y) = 101.739
Var(X) = 0.901
Hypothesis : Beta < 4 Var(X_Residual) = 0.882
Other Params : rybar = +inf
--------------------------------------------------------------------------------
cbar rxbar(Breakdown)
--------------------------------------------------------------------------------
0.000 128.1%
0.100 114.0%
0.200 103.6%
0.300 95.7 %
0.400 89.6 %
0.500 85.0 %
0.600 81.7 %
0.700 79.5 %
0.800 78.8 %
0.900 78.8 %
1.000 78.8 %
--------------------------------------------------------------------------------
为了将上表可视化,我们可以结合 regsensitivity bounds
的命令绘制可识别集合和 4 的交点。
无论是 regsensitivity bounds
还是