Stata：敏感性分析-rcr
👇 连享会 · 推文导航 | www.lianxh.cn
🍎 Stata：Stata基础 | Stata绘图 | Stata程序 | Stata新命令 📘 论文：数据处理 | 结果输出 | 论文写作 | 数据分享 💹 计量：回归分析 | 交乘项-调节 | IV-GMM | 时间序列 | 面板数据 | 空间计量 | Probit-Logit | 分位数回归 ⛳ 专题：SFA-DEA | 生存分析 | 爬虫 | 机器学习 | 文本分析 🔃 因果：DID | RDD | 因果推断 | 合成控制法 | PSM-Matching 🔨 工具：工具软件 | Markdown | Python-R-Stata 🎧 课程：公开课-直播 | 计量专题 | 关于连享会
连享会 · 文本分析 | 爬虫 | 机器学习
作者：李适源 (北京大学)
邮箱：shiyuanli@pku.edu.cn
编者按：本文主要摘译自下文，特此致谢！
Source：Krauth B. Bounding a linear causal effect using relative correlation restrictions[J]. Journal of Econometric Methods, 2016, 5(1): 117-141. -PDF-
目录
1. 敏感性分析的基本原理
1.1 外生性假定下的因果识别
1.2 从点估计值到边界分析
1.3 如何设定内生性偏误的取值范围
2. Stata 操作
2.1 示例数据简介
2.2 OLS 与 RCR 比较
3. 总结与扩展
4. 参考资料
5. 相关推文
温馨提示： 文中链接在微信中无法生效。请点击底部「阅读原文」。或直接长按/扫描如下二维码，直达原文：
1. 敏感性分析的基本原理
1.1 外生性假定下的因果识别
在使用传统计量方法进行因果推断时，我们常常会设定如下的线性模型：
其中， 是我们感兴趣的核心解释变量 (或称处理变量、原因变量)。 代表解释变量 对于结果变量 的因果效应。是一组控制变量。
需要注意的是：出于简便起见，本文论述基于传统计量框架下的 “同质性因果效应”，没有引入 “潜在结果框架” 来严格定义因果效应 (或称处理效应)。另外，本文使用的符号来自原文，与通常介绍因果推断的文献有所不同。
为识别总体因果效应 ，我们通常施加 “条件均值独立” 假定，即在给定控制变量 之后，误差项 均值独立于解释变量 。更进一步，可以假定 均值独立于控制变量 。在传统计量框架中，以下等式就是大家熟悉的外生性假定 (exogeneity assumption)。
注意上式第二个等号对应的假定并不是强假定。这是因为，如果 包含与 有关的因素，我们总能把这些因素从 中提出来，与 进行合并，这只影响到控制变量对应系数的估计是否存在偏误。详细论证可以参考 Stock 和 Waston (2015) 的计量教材 (Section 7.5)。
由于 “均值独立” 可以推出 “不相关”，上述假定条件可以推出误差项 与解释变量 不相关，误差项 与控制变量 也不相关。
此时，使用 对 和 做 OLS 回归，即可得到因果效应系数 的一致估计。
1.2 从点估计值到边界分析
为什么我们有这么强烈的信念，认定解释变量 与误差项 一定不相关呢？在观测性研究中 (observational study)，解释变量的取值并不是随机分配的，往往是个体自己选择的结果。即便引入了可观测的控制变量，我们仍然难以保证误差项 “干净无暇”。换言之，内生性问题仍然可能存在 (尤其是遗漏变量带来的偏误)。
相比于直接假定 与 并不存在相关关系，一个更为 “温和” 的假定是：误差项 当中可能仍存在一些遗漏变量——这些变量既影响 ，又与 相关，但却未能被放入回归中加以控制，因此将会导致因果效应系数的估计偏误。
更重要的是，我们设定误差项与解释变量之间相关性强度的取值范围 (即遗漏变量可能带来的威胁最小是多少，最大是多少)，然后在因果推断存在瑕疵 (内生性) 的情况下，重新估计 对应的系数。此时，我们可以得到系数 的一个边界 (bound)，即 。
其中， 是指在允许误差项与解释变量存在相关性的强度范围内，估计出的 对应系数的最小值；而 是指在允许误差项与解释变量存在相关性的强度范围内，估计出的 对应系数的最大值。
边界 可以告诉我们，因果效应系数的点估计值对于遗漏变量偏误的敏感程度 (sensitivity)。如果对遗漏变量不敏感，那么即使我们将误差项与解释变量之间相关性强度的取值范围设定得较宽，得到的边界也仍然远离 0 值。需要注意的是，这里的边界是指在我们设定的偏误范围内， 所有可能的点估计值集合，它并不是置信区间，故不能直接基于此进行统计检验。
另外，我们也可以求出边界值 对应的 95% 置信下界，以及边界值 对应的 95% 置信上界。从而，将最小可能取值的置信下界和最大可能取值的置信上界拼到一起，就构成了敏感性分析的置信区间 (原作者将之称作 “保守的置信区间”)。
这个置信区间可以近似解读为，在我们允许的误差项与解释变量的相关性强度范围内，该区间有 95% 的可能性覆盖住真实的因果效应系数。此时，可以直接考察这个置信区间是否包含 0 值。如果已经包含了 0 值，则意味着，在我们设定的误差项与解释变量的相关性强度范围之内，因果效应系数的估计结果在 5% 水平上不再显著，即系数估计结果对遗漏变量问题较为敏感 (不够稳健)。
1.3 如何设定内生性偏误的取值范围
了解敏感性分析的基本原理之后，一个关键问题是，应该怎样设定 的取值范围？原文作者提供的方案是，我们无需武断地直接设定 的取值，而是将解释变量 与控制变量 的相关性作为参照，设定 与 比值 的取值范围。
具体来说， 可解读为误差项与解释变量相关性 (内生性偏误幅度) 达到了解释变量与控制变量相关性强度 的 倍。
原文作者建议，可以将 的取值范围设定在 。当 取值为 时，意味着外生性假定成立，误差项与解释变量不存在相关性；当 取值为 时，意味着外生性假定遭到违反，误差项与解释变量的相关性强度，刚好等于解释变量与控制变量的相关性强度。换言之，此时被遗漏在误差项中的变量，对因果识别的重要程度，类似于已被纳入到回归方程中的控制变量。
总结起来，作者将偏误系数 定义为两个相关系数的比值，使得偏误系数的具体取值度量了内生性偏误的相对大小，即内生性偏误达到了解释变量与控制变量相关性的多少倍。因此，作者将这种敏感性分析的思路称作 “相对意义下的相关性限制” (relative correlation restrictions，rcr)。
2. Stata 操作
2.1 示例数据简介
首先，安装 rcr 命令：
. net from http://www.sfu.ca/~bkrauth/code
. net install rcr.pkg, replace
. net get rcr.pkg, replace
然后，通过 help rcr 查看命令语法：
rcr depvar treatmentvar controlvars [if] [in] [weight][, options]
其中，[options] 中的 lambda() 选项可以手动设定偏误系数的取值范围，默认取值范围为 lambda(0 1)。[options] 中的 cluster() 选项可以计算聚类稳健标准误。
我们使用 rcr_example.dta 数据进行演示。该数据取自田纳西州的师生比例研究 (Tennessee Student Teacher Achievement Ratio Study)，简称 STAR 项目。该项目主要用来评估班级规模对学生学业表现的影响。如果大家阅读过 Stock 和 Waston (2014) 的计量教材，那么一定对这份数据的研究情景非常熟悉。
以下分析是基于 Krueger (1999) 的研究。79 所小学被非随机地选中参加了 STAR 项目。在每一所小学当中，1985 年进入幼儿园的学生，被随机分配到三组：小班组、普通班组或有全职教师助理的普通班组。为简化论述，我们将上述三组学生合并为两组：小班组和普通班组。研究目标是考察，小班教学 (相比普通班) 对于学生成绩的因果效应。
. use "http://www.sfu.ca/~bkrauth/code/rcr_example.dta", clear // 导入数据
. describe, fullnames // 查看变量信息和样本量
Contains data from http://www.sfu.ca/~bkrauth/code/rcr_example.dta
Observations: 5,839 Project STAR data, kindergarten
Variables: 11 28 May 2010 09:12
-----------------------------------------------------------------------------
Variable Storage Display Value
name type format label Variable label
-----------------------------------------------------------------------------
SCHID long %9.0g SCHID
TCHID long %9.0g TCHID
SAT double %9.0g SAT
Small_Class double %9.0g Small.Class
Regular_Aide_Class
double %9.0g Regular.Aide.Class
White_Asian double %9.0g White.Asian
Girl double %9.0g Girl
Free_Lunch double %9.0g Free.Lunch
White_Teacher double %9.0g White.Teacher
Teacher_Experience
double %9.0g Teacher.Experience
Masters_Degree double %9.0g Masters.Degree
-----------------------------------------------------------------------------
可以看到，这份示例数据集包含 5839 个学生 (observations)，共 11 个变量。其中，结果变量 SAT 表示样本学生的平均百分位成绩 (average percentile score)。核心解释变量 (处理变量) Small_Class 是一个二值变量，取 1 表示学生所在班级规模为小班，课堂容量在 13-17 人之间；取 0 表示学生所在班级为普通班级规模，课堂容量在 22-25 人之间。
虽然这项研究本身是实验研究，但是 Krueger (1999) 注意到分班的随机化并不完美：有一些学生原本被分配到普通班，但通过家长的努力，成功转到了小班当中；另外，被分配到普通班的学生，有更大概率选择转校 (可能是想去其他学校的小班)。这些都是不完美随机化的证据。因此有必要在回归分析中控制以下变量：
SCHID：表示学生所在学校代码 (可以使用 i.SCHID来设置学校固定效应)。之所以要设置学校固定效应，是由于每所学校是否参与 STAR 项目是非随机的，而随机分班是在每一所学校内进行的；White_Asian：表示学生所在班级中白人和亚裔占比； Girl：表示学生所在班级中女孩占比； Free_Lunch：表示表示学生所在班级中，享受免费午餐人数占比； White_Teacher：表示学生所在班级中，白人教师占比； Teacher_Experience：表示学生所在班级中教师 (平均) 教龄； Masters_Degree：表示学生所在班级中拥有硕士学历的教师占比。
另外，TCHID 类似于学生所在班级的代码。考虑到同一班级内学生的扰动项可能存在聚类相关，因此在下列回归中使用 cluster(TCHID) 来估计聚类稳健的标准误。
2.2 OLS 与 RCR 比较
首先使用 OLS，引入上述控制变量 (并设置学校固定效应)。我们使用 areg 命令进行估计，利用 absorb(SCHID) 来吸收学校固定效应。估计结果如下：
. global Controls "White_Asian Girl Free_Lunch White_Teacher Teacher_Experience Masters_Degree"
. areg SAT Small_Class $Controls, absorb(SCHID) cluster(TCHID)
Linear regression, absorbing indicators Number of obs = 5,839
Absorbed variable: SCHID No. of categories = 79
F(7, 322) = 61.41
Prob > F = 0.0000
R-squared = 0.0927
Adj R-squared = 0.0793
Root MSE = 22.3517
(Std. err. adjusted for 323 clusters in TCHID)
------------------------------------------------------------------------------------
| Robust
SAT | Coefficient std. err. t P>|t| [95% conf. interval]
-------------------+----------------------------------------------------------------
Small_Class | 5.202 1.038 5.01 0.000 3.159 7.244
White_Asian | 8.389 1.363 6.15 0.000 5.707 11.071
Girl | 4.382 0.630 6.96 0.000 3.143 5.621
Free_Lunch | -13.077 0.774 -16.89 0.000 -14.600 -11.555
White_Teacher | -1.090 2.185 -0.50 0.618 -5.388 3.209
Teacher_Experience | 0.265 0.105 2.53 0.012 0.059 0.471
Masters_Degree | -0.600 1.061 -0.57 0.572 -2.687 1.487
_cons | 47.096 2.422 19.45 0.000 42.332 51.860
------------------------------------------------------------------------------------
OLS 估计结果显示，小班教学对应系数为 5.202，在 1% 的统计水平上显著。这在一定程度上说明，平均意义上，小班教学 (相比于普通班级规模) 提升了学生的标准化分位数成绩约 5.2 个百分点。
OLS 估计结果是否对未观测到的遗漏变量敏感？下面使用 rcr 命令来进行敏感性分析。命令格式和 reg 非常相似，不同之处是在选项中加入偏误系数 的取值范围 lambda(0 1)。估计结果如下：
. rcr SAT Small_Class $Controls,lambda(0 1) cluster(TCHID)
RCR model Number of obs = 5839
Lower bound on lambda = 0
Upper bound on lambda = 1
(Std. err. adjusted for 323 clusters in TCHID)
------------------------------------------------------------------------------
| Robust
SAT | Coefficient std. err. z P>|z| [95% conf. interval]
-------------+----------------------------------------------------------------
lambdaInf | 12.311 8.317 1.48 0.139 -3.990 28.611
betaxInf | 8.170 43.646 0.19 0.852 -77.375 93.714
lambda0 | 28.935 145.177 0.20 0.842 -255.607 313.478
betaxL | 5.135 1.359 3.78 0.000 2.472 7.798
betaxH | 5.202 1.031 5.05 0.000 3.182 7.221
-------------+----------------------------------------------------------------
betax | (Conservative confidence interval) 2.472053 7.221435
------------------------------------------------------------------------------
Treatment Variable: Small_Class
Control Variables : White_Asian Girl Free_Lunch White_Teacher
Teacher_Experience Masters_Degree
上表中比较关键的是最后三行系数。倒数第三行 betaxL 对应 的估计值，也就是在我们设定的内生性偏误范围内，回归系数 的最小可能取值。倒数第二行 betaxH 对应 的估计值，即在我们设定的内生性偏误范围内，回归系数 的最大可能取值。
可以看出， 和 这两个边界取值相当接近，也很接近 OLS 回归系数。因此从点估计值来看，在我们允许的偏误范围内 (误差项和解释变量相关性强度最多不超过控制变量和解释变量相关性)，回归系数的点估计值变化不大，非常接近原始的估计结果 5.202 (即施加了外生性假定的 OLS 估计值)。
最后一行为保守置信区间 (也就是对偏误不敏感的置信区间)。这个置信区间就是将 的 95% 置信区间下限，和 的 95% 置信区间上限拼接到一起得到的。保守置信区间约为 ，可以说明在我们允许的偏误范围内，处理变量 Small_Class 的系数估计结果在 5% 水平上仍然显著地区别于 0 值。这再次佐证了估计结果对可能存在的遗漏变量偏误并不敏感。
3. 总结与扩展
本文介绍的相对偏误与边界估计方法，主要适用于同质性因果效应模型，而且一般依赖于线性模型的参数化设定。但是在一般性的因果推断框架中，因果效应往往是因人而异，而且结果变量的条件期望函数可能与控制变量也并非线性关系。此时可能需要借助一些半参数方法或者非参数方法来估计总体平均因果效应 (rcr 命令可能不再适用)。
为此，大家可参考以往文献提供的敏感性分析方法 (Rosenbaum 和 Rubin，1984；Rosenbaum，1995；Vanderweele 和 Arah，2011)。他们使用更具一般性的潜在结果框架来定义因果效应，在该框架下识别和估计因果效应系数，并展开敏感性分析。大家也可以阅读计量界两位大佬的一篇综述 (Imbens 和 Wooldridge，2009)，他们在第 6 章的 6.1 节和 6.2 节详细回顾了敏感性分析 (Sensitivity Analaysis) 的经典方法。
此外，如果要针对工具变量回归的估计结果展开敏感性分析，可以参考 Conley 等 (2012) 提供的分析框架。
4. 参考资料
Conley T G, Hansen C B, Rossi P E. Plausibly exogenous[J]. Review of Economics and Statistics, 2012, 94(1): 260-272. -PDF- Imbens G W, Wooldridge J M. Recent developments in the econometrics of program evaluation[J]. Journal of economic literature, 2009, 47(1): 5-86. -PDF- Krauth B. Bounding a linear causal effect using relative correlation restrictions[J]. Journal of Econometric Methods, 2016, 5(1): 117-141. -PDF- Krueger A B. Experimental estimates of education production functions[J]. The quarterly journal of economics, 1999, 114(2): 497-532. -PDF- Rosenbaum P R, Rubin D B. Reducing bias in observational studies using subclassification on the propensity score[J]. Journal of the American statistical Association, 1984, 79(387): 516-524. -PDF- Rosenbaum, Paul R. 1995(1st edition),2002 (2nd edition). Observational Studies. New York: Springer. Stock, James H., and Mark W. Watson. 2015 (3rd edition). Introduction to econometrics. New Jersey: Pearson. VanderWeele T J, Arah O A. Unmeasured confounding for general outcomes, treatments, and confounders: bias formulas for sensitivity analysis[J]. Epidemiology (Cambridge, Mass.), 2011, 22(1): 42. -PDF-
5. 相关推文
Note：产生如下推文列表的 Stata 命令为：
lianxh 敏感性 不外生 遗漏变量, m
安装最新版lianxh命令：
ssc install lianxh, replace
专题：回归分析 遗漏变量？敏感性分析！新命令sensemakr-T310 专题：IV-GMM IV-估计：工具变量不外生时也可以用！ 专题：内生性-因果推断 因果推断：未测量混杂因素的敏感性分析-T249 Selection Ratio：帮你解决头疼的遗漏变量偏误
New！ Stata 搜索神器：
lianxh和songblGIF 动图介绍
搜： 推文、数据分享、期刊论文、重现代码 ……
👉 安装：
. ssc install lianxh
. ssc install songbl
👉 使用：
. lianxh DID 倍分法
. songbl all
🍏 关于我们
连享会 ( www.lianxh.cn，推文列表) 由中山大学连玉君老师团队创办，定期分享实证分析经验。 直通车： 👉【**百度一下：**连享会】即可直达连享会主页。亦可进一步添加 「知乎」,「b 站」,「面板数据」,「公开课」 等关键词细化搜索。