遗漏变量？敏感性分析！新命令sensemakr-T310

发布时间：2021-05-05 阅读 6062

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

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

New！ lianxh 命令发布了：
随时搜索推文、Stata 资源。安装命令如下：
. ssc install lianxh
详情参见帮助文件 (有惊喜)：
. help lianxh

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

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

⛳ Stata 系列推文：

全部 | Stata入门 | Stata教程 | Stata资源 | Stata命令
计量专题 | 论文写作 | 数据分享 | 专题课程
结果输出 | Stata绘图 | 数据处理 | Stata程序
回归分析 | 面板数据 | 交乘项-调节 | IV-GMM
内生性-因果推断 | 倍分法DID | 断点回归RDD | PSM-Matching | 合成控制法
Probit-Logit | 时间序列 | 空间计量 | 分位数回归 | 生存分析 | SFA-DEA
文本分析-爬虫 | Python-R-Matlab | 机器学习
Markdown | 工具软件 | 其它

作者： 邹恬华 (对外经济贸易大学)
邮箱： 201802060@uibe.edu.cn

Source: Cinelli, Carlos, Jeremy Ferwerda, and Chad Hazlett. 2020. “Sensemakr: Sensitivity Analysis Tools for OLS in R and Stata.”, -Link-

存在可能的遗漏变量问题且数据不可得时，可以使用敏感性分析进行稳健性检验。
即检验：遗漏变量需要达到多强才能使得之前的研究结果改变？
敏感性分析可以通过 Stata 命令 sensemakr 进行。
开始敏感性分析之旅吧！

1. 我们为什么要做敏感性分析？

1.1 一个小例子

在 2003 年和 2004 年，达尔富尔政府策划了一场针对平民的恐怖暴力运动，估计造成 20 万人丧生。我们期望评估，在此事件中直接受伤或致残，是会使人更具有复仇心理，从而不愿与施暴者和平相处，还是更有动力看到这种暴力结束，支持呼吁和平？ Darfur 数据集(Hazlett, 2019) 是针对乍得东部达尔富尔难民的一项调查数据，我们建立 OLS 回归模型如下：

P e a c e F a c t o r = α + β D i r e c t l y H a r m e d + γ F e m a l e + δ O t h e r C o n t r o l s + V i l l a g e + ε

PeaceFactor：和平意愿，被解释变量
DirectlyHarmed：是否被直接伤害的虚拟变量，我们关心的核心解释变量
Female：是否为女性的虚拟变量，控制变量之一
OtherControls：除 Female 外的控制变量，包括 Age 年龄， HerderDar 是否为达尔富尔牧民， FarmererDar 是否为达尔富尔农民， PastVoted 是否之前在选举中投票
Village：被调查者居住村庄的固定效应
$β$ ：我们关心的系数，即被直接伤害如何影响村民的和平意愿

回归结果表明，系数 $β$ 显著为正，即被直接伤害平均而言更会激发村民和平意愿。但有学者认为以上模型可能存在遗漏变量问题：

恐怖暴力活动可能更有可能击中村庄的中心，而那些在中心的人可能对和平持有不同的态度。
袭击者更易攻击具有某种个人特征的人，比如更高的财富水平。
个人(先前)的政治态度可能导致他们采取某种行动，从而使他们在袭击期间面临更大的风险。
更复杂的是，以上因素可能相互作用，或产生其他非线性影响。

问题在于上述可能的遗漏变量数据不可得。在此情况下，想要说明上述结果的稳健性，也可以采取以下思路：

当遗漏变量达到何种强度时，之前的结果才不可信？若这种要求显然过于苛刻，则说明之前的结果大致可信。
此种思路即为 「敏感性分析」。

注：“强度”指与核心解释变量相关性，下同。

1.2 敏感性分析解决的问题

结合上述的小例子，我们总结出敏感性分析想要解决的问题：

遗漏变量需要达到多强才能推翻之前的研究结果？
在最坏的情况下（遗漏变量解释被解释变量所有剩余方差），之前的研究结果有多稳健？
遗漏变量需要达到多强才能使得之前的研究结果改变一定的水平？
遗漏变量相比现有某个变量需要达到多强才能使得之前的研究结果改变？

2. 开始使用 sensemakr 吧！

2.1 命令安装与数据下载

sensemakr 是 Cinelli et al.(2020) 开发的针对 OLS 估计的敏感性分析的 Stata 命令，同时提供功能类似的 R 语言版本。可以通过向 Stata 命令行输入以下命令下载：

. ssc install sensemakr, replace all

前文例子中提到的 Darfur 数据集可以在作者的 Github 主页下载：https://github.com/resonance1/sensemakr-stata/ (Hazlett, 2019)。亦可直接执行如下命令，将数据下载到当前工作路径下：

. net get sensemakr.pkg

checking sensemakr consistency and verifying not already installed...

copying into current directory...
      copying  darfur.dta
ancillary files successfully copied.

打开数据文件：

use darfur.dta, clear

数据集中包含 1.1 节小例子中的全部数据。

2.2 基础语法结构

sensemakr depvar covar [if] [in], treat(varlist) benchmark(varlist)

depvar：被解释变量
covar：解释变量
treat：此处申明最为关心的核心解释变量，这也意味着填写 covar 时的顺序不重要。
benchmark：与可能存在的遗漏变量进行对比的变量，即最后得出的结果为遗漏变量需要达到对比变量强度的多少倍才能影响之前的研究结果。

在上述的 Darfur 例子中，该语句可以写为：

sensemakr peacefactor directlyharmed age farmer herder ///
          pastvoted hhsize female i.village_,          ///
          treat(directlyharmed) benchmark(female)

Stata 汇报结果如下：

. net get sensemakr.pkg  // 下载数据
. use darfur.dta, clear  // 导入数据

. sensemakr peacefactor directlyharmed age farmer herder ///
            pastvoted hhsize female i.village_,          ///
            treat(directlyharmed) benchmark(female)

                                                            DOF    =     783
                                                            q      =    1.00
                                                            alpha  =    0.05
                                                            reduce =    TRUE
                                                            H0     =       0

 Treatment      |     Coef.      S.E.      t(H0)    R2yd.x     RV_q    RV_qa
----------------+-----------------------------------------------------------
 directlyharmed |    0.0973    0.0233     4.1844    0.0219   0.1388   0.0763


 Bound         | R2dz.x  R2yz.dx    Coef.     S.E.   t(H0)  Lower CI Upper CI
---------------+--------------------------------------------------------------
 1.00x  female | 0.0092   0.1246   0.0752   0.0219  3.4389   0.0323   0.1182
 2.00x  female | 0.0183   0.2493   0.0529   0.0204  2.6002   0.0130   0.0929
 3.00x  female | 0.0275   0.3741   0.0304   0.0187  1.6281  -0.0063   0.0670

 Extreme Bound | R2dz.x   R2yz.dx     Coef.
---------------+----------------------------
 1.00x  female | 0.0092    1.0000    0.0347
 2.00x  female | 0.0183    1.0000    0.0084
 3.00x  female | 0.0275    1.0000   -0.0121

部分结果和统计量的含义及解读如下：

R2yd.x：Partial R2 of the treatment with the outcome (R2yd.x):
An extreme confounder (orthogonal to the covariates) that explains 100 percent of the residual variance of the outcome, would need to explain at least 2.19 percent of the residual variance of the treatment to fully account for the observed estimated effect.
RV_q：Robustness Value, q = 1.00 (RV_q):
Unobserved confounders (orthogonal to the covariates) that explain more than 13.88 percent of the residual variance of both the treatment and the outcome are strong enough to bring the point estimate to 0 (a bias of 100 percent of the original estimate). Conversely, unobserved confounders that do not explain more than 13.88 percent of the residual variance of both the treatment and the outcome are not strong enough to bring the point estimate to 0.
RV_qa：Robustness Value, q = 1.00, alpha = 0.05 (RV_qa):
Unobserved confounders (orthogonal to the covariates) that explain more than 7.63 percent of the residual variance of both the treatment and the outcome are strong enough to bring the estimate to a range where it is no longer 'statistically different' from 0 (a bias of 100 percent of the original estimate), at the significance level of alpha = 0.05. Conversely, unobserved confounders that do not explain more than 7.63 percent of the residual variance of both the treatment and the outcome are not strong enough to bring the estima te to a range where it is no longer 'statistically different' from 0, at the significance level of alpha = 0.05
Bounds on Omitted Variable Bias:
The table shows the maximum strength of unobserved confounders, bounded by a multiple of the observed explanatory power of the chosen benchmark covariate(s) with the treatment and the outcome.

2.2.1 无对比变量结果

我们逐个解释汇报指标的经济含义：

R2yd.x： $R_{Y \sim D ∣ X}^{2}$ ，在控制其它控制变量情况下，核心解释变量 D 对被解释变量 Y 的偏 $R^{2}$ 。此处等于 0.0219，意味着在最坏条件下（即遗漏变量解释被解释变量所有剩余方差），遗漏变量需要至少解释核心解释变量 2.19%的剩余方差，才能彻底去除之前观察到的核心解释变量引起的效应。
RV_q：使得估计系数恰好为零的稳健性值。此处等于 0.1388，意味着遗漏变量需要最少同时解释 13.88%的核心解释变量与被解释变量的剩余方差，才能彻底去除之前观察到的核心解释变量引起的效应。
RV_qa：使得估计系数不再在 95%水平上显著的稳健性值。此处等于 0.0763，意味着遗漏变量需要最少同时解释 7.63%的核心解释变量与被解释变量的剩余方差，才能使得 95%置信区间下限为零。

2.2.2 对比变量结果（假设遗漏变量强度与对比变量相同）

然而，上述指标不够直观，并且我们很难论证为什么遗漏变量确实能够达到上述要求，并以此说明我们模型的稳健性。为此我们引入对比变量，尝试论证当遗漏变量强度只要小于对比变量时，之前的估计结果是有效的。

在 Darfur 这个例子当中，我们将对比变量设置为 female，意味着我们尝试论证当遗漏变量强度只要小于 female 时，之前的估计结果是有效的。这在一定程度上是说得通的，因为当地针对特定性别的伤害行为要多于针对城市中心人群/富有人群/特定政治立场人群。

Stata 汇报结果的倒数第二个表格汇报了有关对比变量的信息：

R2yz.dx： $R_{Y \sim Z ∣ X, D}^{2}$ ，在控制所有控制变量情况下，遗漏变量 Z 对被解释变量 Y 的偏 $R^{2}$ 。此处等于 0.1246，小于RV_q的值。
R2dz.x： $R_{D \sim Z ∣ X}^{2}$ ，在控制其它控制变量情况下，遗漏变量 Z 对核心解释变量 D 的偏 $R^{2}$ 。此处等于 0.0092，小于RV_q的值。
R2yz.dx、R2dz.x均小于RV_q意味着，当遗漏变量强度小于 female 时，之前的估计结果是有效的。 R2dz.x 的值甚至低于 R2yd.x ，意味着最坏情况下（即遗漏变量解释被解释变量所有剩余方差），当遗漏变量强度小于 female 时，之前的估计结果是有效的。

2.2.3 对比变量结果（假设遗漏变量强度是对比变量的 n 倍）

为了使得上述对比变量结果更具有说服力，我们继续探讨，若遗漏变量的强度时对比变量 female 的 n 倍时，之前的研究结果是否仍然成立。

我们使用 sensemakr 命令中的 contourplot 与 tcontourplot 绘制系数、t 统计量等值线图：

sensemakr peacefactor directlyharmed age farmer herder pastvoted hhsize ///
female i.village_, treat(directlyharmed) benchmark(female) contourplot

sensemakr peacefactor directlyharmed age farmer herder pastvoted hhsize ///
female i.village_, treat(directlyharmed) benchmark(female) tcontourplot

绘图结果如下，上图为系数 $β$ 等值线图，下图为 t 统计量等值线图。

图中横纵轴坐标分别为 $R_{D \sim Z ∣ X}^{2}$ 与 $R_{Y \sim Z ∣ X, D}^{2}$ 。

上图等值线为回归系数 $β$ 值，红线为 $β = 0$ 。四个数值点分别指不加入遗漏变量、加入 female 相同强度的遗漏变量、加入female2 倍强度的遗漏变量、加入female3 倍强度的遗漏变量的情况。四个数值点的小括号中汇报了其对应计量模型的回归系数 $β$ 值，可以观察到均大于零，也即全部四个数值点位于红线左侧。这意味着即使加入female3 倍强度的遗漏变量，也不会使得原估计系数由正转负。

下图等值线为 t 统计量，红线为 1.96（95%置信区间临界值）。四个数值点与上图相同，分别指加入female0~3 倍强度的遗漏变量的情况。四个数值点的小括号中汇报了其对应计量模型的 t 统计量。加入female0~2 倍强度的三个数值点 t 统计量小于 1.96，在红线左侧，加入female3 倍强度的数值点 t 统计量大于 1.96，在红线右侧。这意味着即使加入female2 倍强度的遗漏变量，也不会使得原估计系数由显著变为不显著；但当加入female3 倍强度的遗漏变量，原估计系数变得不显著。

此外，sensemakr 命令还允许绘制最坏情况图，代码如下：

sensemakr peacefactor directlyharmed age farmer herder pastvoted hhsize ///
female i.village_, treat(directlyharmed) benchmark(female) extremeplot

绘图结果如下：

三条数值线分别指当遗漏变量解释剩余方差的 100%、75%和 50%的情况，即 $R_{Y \sim Z ∣ X, D}^{2}$ 分别等于 1、0.75 与 0.5 的情况。横轴的小黑点分别指遗漏变量是female1~3 倍强度的情况。

图中可以看到，即使是最坏情况（遗漏变量解释剩余方差的 100%，图中实线），遗漏变量都需要female2 倍以上强度才能推翻之前结论。而对于 75%与 50%的情况，遗漏变量则需要female3 倍以上强度。

2.3 其它进阶选项

sensemakr 还有其它选项可供选择，主要是修改上述分析过程中的一些参数。以下列出这些选项：

alpha(real)：设置显著性水平，默认值 0.05
gbenchmark(varlist)：添加单个对比变量使用前文的 benchmark(varlist) ，而当对比变量有一组时，使用gbenchmark(varlist)
gname(string)：给一组对比变量取名，默认值为用组内变量名用连字符连接
kd(numlist)与ky(numlist)：遗漏变量分别相比核心解释变量与被解释变量的强度，默认值为（1 2 3）。用户单独设定 kd 与 ky 其中一个时，另一个会自动被设置成相同值；若用户两个都设定，则必须设置成相同值。
latex(filename)：将结果输出到 filename.tex
noreduce：遗漏变量对回归系数 $β$ 的影响方向，默认值为向 $β = 0$ 方向影响
q(real)：指定估计的哪一部分需要被解释为有问题的。默认值为 1，这意味着将当前效果 100％有问题，核心解释变量的真实效果为 0。
r2yz(numlist)：默认值为 1，意味着遗漏变量解释剩余方差的 100%。如果输入(.5 .75)意味着输出遗漏变量分别解释剩余方差的 50%和 75%的结果。
surppress：不要汇报详细概念描述与解释 等值线图相关：
clines：等值线的条数，默认值为 7
clim(numlist)：横纵坐标的最小值和最大值，最大范围为(0 1) 最坏情况图相关：
r2yz：设置线条 $R_{Y \sim Z ∣ X, D}^{2}$ 取值，默认值为(1 0.75 0.5)，最多可以设置四条线
elim(numlist)：设置横坐标的最小值和最大值，最大范围为(0 1)。纵坐标范围自动选择，不可被设定。

3. 反思：敏感性分析应该用来做什么？

作者担忧敏感性分析可能被滥用，并且可能会将其用作应“通过”的“稳健性测试”，类似于目前显著性检验的滥用。
敏感性分析实际上告诉我们的是，为了维持目前的结论，我们选择相信一个怎样的假设。在本文的例子当中，这个我们选择相信的假设是，遗漏变量解释全部剩余方差且与核心解释变量相关的强度相比 female 强两倍以上，这个事件不太可能发生。
但敏感性分析不能给出遗漏变量实际上是否满足这一假设，不能告诉我们目前的结论是否正确。这一讨论最终还是要依赖实际研究领域的知识，而非统计学知识。
敏感性分析实际上给出了量化后的推翻目前结论的条件。这使得有序的学术辩论成为可能。

4. 参考文献

Cinelli, Carlos, Jeremy Ferwerda, and Chad Hazlett. 2020. “Sensemakr: Sensitivity Analysis Tools for OLS in R and Stata.”, -Link-

5. 相关推文

Note：产生如下推文列表的 Stata 命令为：
lianxh 稳健性检验安慰剂检验
安装最新版 lianxh 命令：
ssc install lianxh, replace

专题	嘉宾	直播/回看视频
⭐ 最新专题		文本分析、机器学习、效率专题、生存分析等
研究设计	连玉君	我的特斯拉-实证研究设计，-幻灯片-
面板模型	连玉君	动态面板模型，-幻灯片-
面板模型	连玉君	直击面板数据模型 [免费公开课，2小时]