敏感性分析B-Stata实操:控制变量内生时的系数敏感性分析-regsensitivity

发布时间:2022-07-07 阅读 1591

Stata连享会   主页 || 视频 || 推文 || 知乎 || Bilibili 站

温馨提示: 定期 清理浏览器缓存,可以获得最佳浏览体验。

New! lianxh 命令发布了:
随时搜索推文、Stata 资源。安装:
. ssc install lianxh
详情参见帮助文件 (有惊喜):
. help lianxh
连享会新命令:cnssc, ihelp, rdbalance, gitee, installpkg

课程详情 https://gitee.com/lianxh/Course

课程主页 https://gitee.com/lianxh/Course

⛳ Stata 系列推文:

PDF下载 - 推文合集

作者:陈卓然 (中山大学)
邮箱chenzhr25@mail2.sysu.edu.cn


目录


1. 引言

因果效应识别的最大挑战之一就是遗漏变量。在线性模型中,当一个控制变量与遗漏变量相关时,这个控制变量的系数将是有偏的。为看清这一点,不妨回忆一下最简单的二元线性回归模型。假定真实模型如下:

从而正确的估计结果为:

然而由于某种原因,我们实际估计的模型如下:

从而错误的估计结果为:

我们将式 (1) 和式 (2) 写成离差的形式:

将式 (3) 左右两边同时乘以 (x1ix¯1),并对 i 求和可得:

简单整理之后可得:

因此,

其中,δ~1 实际上是 x2 对 x1 回归的系数,刻画了 x2 和 x1 之间的相关系数,即:

因此遗漏变量偏误的表达式为:

等式 (4) 的推导利用到 β^1 和 β^2 的无偏性。此外 β~1 和 β^1 的方差的表达式为:

其中,σ~2=i=1n(yiβ~0β~1x1)2/(n2)σ^2=i=1n(yiβ^0β^1x1β^2x2)2/(n3)

这其中有几点值得注意,如果 x1 和 x2 高度相关,但是 x2 对 y 的影响很小,则 β~ 的偏误事实上很小,也就是说 β~ 和 β^ 之间很接近。此时 Var(β~1X)<Var(β^1X)。相反,如果 x1 和 x2 之间相关性很小,但是 x2 对 y 的影响很大,此时 β~ 的偏误事实上也很小,但是 Var(β~1X)>Var(β^1X)

当然对于上述遗漏变量偏误的推导也有一些别的方法,比如:

不过上述的推导过程尽管看似复杂,但是很容易推广到更多元的情形之下,因此我们在这里将其较为详细地列示出来,当然如果读者掌握更高级的矩阵知识的话,也可以有更为简便的办法,在这里就不再赘述了。

由于在实证分析中,我们很难将全部的解释变量都纳入进模型当中,因此将不可避免地受到遗漏变量的影响。尽管如此,我们仍可以分析结果对于遗漏变量偏误的敏感性。

一些经典的研究遗漏变量对于实证结果的影响的文献 (Altoji,2005;Oster,2019) 均要求外生控制变量假设,即遗漏变量必须与控制变量不相关,而这一假设通常是一个很强而且不太合理的假设。本文主要探讨在 Diegert 等 (2022) 提出的敏感性分析方法,该方法允许内生控制变量的存在。在这一框架下,有两个核心问题:

  • 在放松的假设下,βlong 的一致估计是什么?或者说在备择假设下,βlong 的边界是什么?
  • 在 βlong 的一个假设不被拒绝之前,我们可以在多大程度上放松外生性假设?也就是截断点 (breakdown point):在一个假设被推翻之前,基准假设的最大可放松程度。我们希望这个程度越大越好,越大表示我们的结论越稳健。

接下来,我们将使用命令 regsensitivity 对上述两种敏感性进行分析。

2. 方法初探

这里采用 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 ]
--------------------------------------------------------------------------------

事实上,上述结果已经回答了我们在引言中提到的两个问题:βlong 在一系列假设之下的边界在哪里;βlong 在何时会被推翻。下面我们分别对于这两点展开讨论。

3. 边界

考虑如下模型:

其中,(Y,X,W0,W1) 是可观测的,W2 是遗漏变量,可能与 (X,W0,W1) 相关。注意此处 βlong 的记号,我们下文还会提到 βshort 和 βmed。现将其含义总结于下表之中:

符号 含义
βlong Y 对 (1,X,W0,W1,W2) 的回归中 X 的系数
βmed Y 对 (1,X,W0,W1) 回归中 X 的系数
βshort Y 对 (1,X,W0) 回归中 X 的系数

(X,W0,W1,W2) 的联合分布约束是由三个标量敏感性系数 (r¯X,r¯Y,c¯) 所控制。给定可观测变量 (Y,X,W0,W1) 的联合分布以及敏感性参数值,Diegert 等 (2022) 介绍了如何计算 βlong 在可识别集合中的的上下边界 BI(r¯X,r¯Y,c¯)

所谓可识别集合是指在现有假设下给定可观测数据的分布,βlong 可以被一致估计的值。当 r¯X>0 且 r¯Y>0 时 βlong 无法被点估计,因此我们需要估计其边界。在 Stata 中要计算BI(r¯X,r¯Y,c¯) 可以采用 regsensitivity bounds 的子命令完成。

regsensitivity bounds depvar indepvar controls, options...

其中,depvar 是被解释变量 Yindepvar controls 是解释变量 (X,W0,W1)。值得注意的是,在 regsensitivity 中,解释变量的顺序至关重要。第一个解释变量 X 是我们感兴趣的变量,也是要进行敏感性分析的变量,而控制变量则是一些我们并不关心的变量。

默认情况下,regsensitivity bounds 计算在 c¯ 和 r¯Y 保持不变的前提下,对于 r¯X 的一系列值,对应的边界。默认情况下 c¯=1r¯Y=+,为指定不同的 c¯ 可以使用 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) 的选项,这一选项是在指明哪些控制变量被放入 W1 而不是 W0. 这些控制变量被称作比较控制变量 (comparison controls), 他们是用来校准敏感性参数 (r¯X,r¯Y,c¯) 的。通过将更多的变量加入比较控制变量集中后,在给定敏感性参数时,βlong 的可识别集合通常而言会变大。而当这一选项被省略掉之后,全部的控制变量都将被包括进 W1 中,比如说

. 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

不难发现,当全部的控制变量都被包含进 W1 和 r¯X=1.336 时,可识别集合就变成了全体实数集 R, 而当 W1 中去除州层面的固定效应之后,只有当 r¯X=4.063 时,可识别集合才会变为全体实数集 R。为直接比较两种 W1 选择下的边界,我们可以手动地设置 rxbar 值,使其在两种情形下是相同的。

由于 regsensitivity bounds 的回归结果被存储在 e(idset_table) 中,因此我们可以提取其中的 rxbar 的值,然后将其作为选项添加到从 W1 中去除州层面的控制变量之后的敏感性回归命令中,即:

. 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 ]
--------------------------------------------------------------------------------

不难看出,当全部的控制变量都被包括进 W1 时,每一个相同的 r¯X 对应的边界要更窄。为了比较 c¯ 的不同值,我们可以将一个数列放进 cbar 的选项当中。此时 regsensitivity bounds 命令将会绘制一张图,而不是以表格的形式呈现在命令窗口

. regsensitivity bounds `y' `x' `w',compare(`w1') cbar(0(.2)1)

在默认情况下,regsensitivity plot 将绘制到对于每一个 c¯ 可识别集合变为全体实数集 r¯X。当然我们也可以限制 r¯X 的取值范围来考察可识别集合与 0 相交的点。

4. 截断前沿

regsensitivity bounds 的结果显示了对于一个给定假设下,参数 βlong 的截断点。对于一个假设 βlongBR,截断点是敏感性参数 r¯X 集合的下界。在 r¯X 这个集合中,βlongBR 这一假设对于可识别集合中的每一个 β 都不成立,正式的定义如下:

regsensitivity 可以处理的假设形式为:对于任意的 bβlongb。默认情况下,假设形式为 sign(βlong)=signβmed。其中 βmed 是 Y 对 (1,X,W0,w1) 回归中 X 的系数。在这个例子中由于 βmed>0,因此默认情况下的假设为 β>0

regsensitivity bounds 显示的是去除掉州层面的固定效应的 W1r¯Xbp(.1,+;(,0])=1.195。为进一步看清截断点是如何随着敏感性参数 c¯ 的变化而变化的,我们可以使用 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

为了检验给定某个 b 时,βlong>b 这样的假设,我们可以指明 beta(b,lb)。其中,lb 是 lower bound 的意思。选项 beta 也可以接受一个数列以检验一系列的假设。例如,对于一系列的 b 值,我们想检验 βlong>b,可以输入以下代码:

. 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)来检验 β<b。其中,ub 表示 upper bound。例如,我们想检验对于一系列的 c¯, β<4,可以输入以下代码:

. 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 的交点。

5. 回归结果的解读

无论是 regsensitivity bounds 还是