RDD：断点回归的非参数估计及Stata实现
温馨提示： 定期 清理浏览器缓存，可以获得最佳浏览体验。
作者： 崔颖（中央财经大学）
Source: Non-Parametric Regression Discontinuity (Francis, 2013)
目录
断点回归 (RDD) 的应用日益广泛。在 「连享会-RDD 专辑1」 和 「连享会-RDD 专辑2」 中，我们对 RDD 的原理、平滑性检验方法等主题进行了介绍。
本文介绍 RDD 的非参数估计在 Stata 中的实现方法，主要采用 rdrobust 命令。该命令是由哥伦比亚大学 Sebastian Calonico 教授、普林斯顿大学 Matias D. Cattaneo 教授及众合作者共同开发。Google 的网页 RD software package 提供了丰富的学习资料，包括许多相关论文的原始数据及复制结果代码。
温馨提示： 文中链接在微信中无法生效。请点击底部「阅读原文」。
1. 命令安装与方法介绍
net install rdrobust, from(http://www-personal.umich.edu/~cattaneo/rdrobust)
RD 可以用来识别自然实验或结构性政策变化附近的局部处理效应。
例如，如果你关心政府奖金对大学入学情况有怎样的影响，你可能会想要将那些获得政府奖金的学生和未获得政府奖金的学生进行比较。但这种方法是存在问题的，因为获得政府奖金的低收入家庭学生与未获得政府奖金的学生可能在多方面均存在差异。
应用 RD 方法的前提条件是个人不能通过合理低报收入水平而获得政府奖金，那些在断点附近的人自报收入分布情况应该和非断点附近的人基本上保持一致。
如果政府奖金资格确定的收入线是未知的，那么，此前提条件可能是合理的。即使学生会系统性地低报他们的收入，但因他们并不知道实际确认资格的收入分界线，可以认为那些收入在断点上下的学生随机抽取自相同的池子，仅是否收到政府奖金这一项差异。
缺乏实验数据的计量经济学识别方法往往需要建立在外生性假定基础之上。也就是说，x 对 y 的影响与误差项 u 不相关。在外生变量直接导致被解释变量变化的情况下，回归识别因果效应才是充分有效的。
在上述例子中，显然，不能简单地将 y (GPA、出勤率、毕业率等)的变化归结为政府奖金的功劳，因为那些收到奖金和未收到奖金的学生存在多方面差异。然而，由于确认资格的收入分界线是未知的，在断点两侧小邻域内的个体可以被视为是相同的。因此，我们有理由认为未知的收入线外生随机地将断点附近的个体分成了两组，一组收到了政府奖金，一组未收到。
2. 模拟生成非线性相关数据
这里，我们假设被解释变量与收入的关系是非线性的 (线性相关性的举例和分析可以参见 Sharp Regression Discontinuity Example and Limitations)。Stata 随机生成一些非线性相关的自变量收入 income 和因变量学习表现 perfo 并绘制二者散点关系图。
温馨提示： 文中链接在微信中无法生效。请点击底部「阅读原文」。
clear
set obs 10000
gen income=3^((runiform()-0.75)*4)
label var income "Reported Income"
sum income
Variable | Obs Mean Std. Dev. Min Max
-------------+---------------------------------------------------------
income | 10,000 .6789349 .7606786 .0370671 2.999232
gen perfo=ln(income)+sin((income-r(min))/r(max)*4*_pi)/3+3
label var perfo "Performance Index - Base"
scatter perfo income
下图展示了自变量与因变量的非线性关系：
现在，让我们加入一些随机扰动。
gen perf1=perfo+rnormal()*0.5
label var perf1 "Performance Index - with noise"
scatter perf1 income
接着，使用命令rcspline，可以将局部平均的学习表现视作收入的三次样条 (Cubic Spline) 函数。
ssc install rcspline
rcspline perf1 income, nknots(7) showknots title(Cubic Spline)
此时，样条曲线是平滑的。接下来，让我们在 0.5 处设置一个断点。
gen grant=income<0.5
sum grant
Variable | Obs Mean Std. Dev. Min Max
-------------+---------------------------------------------------------
grant | 10,000 .5921 .491469 0 1
*样本中大约有59%低收入学生是具备获得政府奖金资格的
*现在加入政府奖金对学生表现的正向效应
*首先生成以政府奖金资格确认收入线为中心的收入变量
gen income_center=income-0.5
gen perf2=perf1+0.5*grant-0.1*income_center*grant
*这样政府奖金对低收入学生将更加有效
label var perf2 "Observed Performance"
3. 进行非参数估计
首先，让我们尝试使用传统的普通最小二乘回归。
reg perf2 income grant
Source | SS df MS Number of obs = 10,000
-------------+---------------------------------- F(2, 9997) = 5734.77
Model | 7041.8845 2 3520.94225 Prob > F = 0.0000
Residual | 6137.80314 9,997 .613964504 R-squared = 0.5343
-------------+---------------------------------- Adj R-squared = 0.5342
Total | 13179.6876 9,999 1.31810057 Root MSE = .78356
------------------------------------------------------------------------------
perf2 | Coef. Std. Err. t P>|t| [95% Conf. Interval]
-------------+----------------------------------------------------------------
income | .9003367 .0168801 53.34 0.000 .8672482 .9334251
grant | -.3767682 .0261265 -14.42 0.000 -.4279814 -.325555
_cons | 1.924569 .0267009 72.08 0.000 1.87223 1.976909
------------------------------------------------------------------------------
显然，估计结果是错误的。政府奖金 grant 的估计系数为负，表明政府奖金会阻碍学习表现，这显然与常理相违背。
接下来，我们尝试使用 RD 命令rdrobust。
*-默认断点在 0 点处，因此我们使用中心化后的变量 income_centered
rdrobust perf2 income_center
Sharp RD estimates using local polynomial regression.
Cutoff c = 0 | Left of c Right of c Number of obs = 10000
-------------------+---------------------- BW type = mserd
Number of obs | 5921 4079 Kernel = Triangular
Eff. Number of obs | 683 530 VCE method = NN
Order est. (p) | 1 1
Order bias (q) | 2 2
BW est. (h) | 0.129 0.129
BW bias (b) | 0.197 0.197
rho (h/b) | 0.652 0.652
Outcome: perf2. Running variable: income_center.
--------------------------------------------------------------------------------
Method | Coef. Std. Err. z P>|z| [95% Conf. Interval]
-------------------+------------------------------------------------------------
Conventional | -.48486 .06467 -7.4971 0.000 -.611614 -.358102
Robust | - - -6.1641 0.000 -.633493 -.327828
--------------------------------------------------------------------------------
回归结果中估计系数为负，原因是 RD 方法通常默认断点处处置变量由 0 变为 1，与本案例中「收入高于收入线，获得政府奖金的资格取消」正好相反。因此，我们需要改变 RD 估计量的符号方向。这可以通过将收入取相反数来实现。
gen nincome_center=income_center*(-1)
rdrobust perf2 nincome_center
Sharp RD estimates using local polynomial regression.
Cutoff c = 0 | Left of c Right of c Number of obs = 10000
-------------------+---------------------- BW type = mserd
Number of obs | 4079 5921 Kernel = Triangular
Eff. Number of obs | 530 683 VCE method = NN
Order est. (p) | 1 1
Order bias (q) | 2 2
BW est. (h) | 0.129 0.129
BW bias (b) | 0.197 0.197
rho (h/b) | 0.652 0.652
Outcome: perf2. Running variable: nincome_center.
--------------------------------------------------------------------------------
Method | Coef. Std. Err. z P>|z| [95% Conf. Interval]
-------------------+------------------------------------------------------------
Conventional | .48486 .06467 7.4971 0.000 .358102 .611614
Robust | - - 6.1641 0.000 .327828 .633493
--------------------------------------------------------------------------------
同时，我们可以将 rdrobust 新命令的回归结果与 Stata 传统 RD 回归命令 rd 的结果相比较。
rd perf2 nincome_center
Two variables specified; treatment is
assumed to jump from zero to one at Z=0.
Assignment variable Z is nincome_center
Treatment variable X_T unspecified
Outcome variable y is perf2
Estimating for bandwidth .1832282339354582
Estimating for bandwidth .0916141169677291
Estimating for bandwidth .3664564678709164
------------------------------------------------------------------------------
perf2 | Coef. Std. Err. z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
lwald | .4860547 .0539727 9.01 0.000 .3802703 .5918392
lwald50 | .4929885 .0744274 6.62 0.000 .3471136 .6388635
lwald200 | .5737225 .0377565 15.20 0.000 .4997212 .6477238
------------------------------------------------------------------------------
rd 命令的好处在于它可以同时汇报出不同带宽下的估计结果。此处，带宽 (100 50 200) 分别代表最小化 MSE (mean squared error) 的最优带宽 (100%) 及其一半与两倍的带宽 (50%, 200%)。
我们可以将不同带宽的估计结果绘制在一张图上。
gen effect_est = .
label var effect_est "Estimated Effect"
gen band_scale = .
label var band_scale "Bandwidth as a Scale Factor of Bandwidth that Minimizes MSE"
forv i = 1/16 {
rd perf2 nincome_center, mbw(100 `=`i'*25')
if `i' ~= 4 replace effect_est = _b[lwald`=`i'*25'] if _n==`i'
if `i' == 4 replace effect_est = _b[lwald] if _n==`i'
replace band_scale = `=`i'*25' if _n==`i'
}
gen true_effect = .5
label var true_effect "True effect"
two (scatter effect_est band_scale) (line true_effect band_scale)
从图上可以看出，最小化 MSE (即 100%) 带宽估计出的结果最为准确，且估计系数在其附近区间内也相对稳定。
总结
本文介绍了 Stata 实现断点回归的最新命令 rdrobust。选取政府奖学金影响学生学业表现的案例，通过数值模拟随机生成观测数据，分别运用 rdrobust 和 rd 命令对模拟数据进行回归分析并比较回归结果。
参考资料
1. RD Software Packages (https://sites.google.com/site/rdpackages/).
2. Calonico S, Cattaneo M D, Farrell M H, et al. rdrobust: Software for regression-discontinuity designs[J]. The Stata Journal, 2017, 17(2): 372-404.PDF
相关课程
连享会-直播课 上线了！
http://lianxh.duanshu.com
免费公开课：
直击面板数据模型 - 连玉君，时长：1小时40分钟 Stata 33 讲 - 连玉君, 每讲 15 分钟. 部分直播课 课程资料下载 (PPT，dofiles等)
课程一览
支持回看，所有课程可以随时购买观看。
| 专题 | 嘉宾 | 直播/回看视频 |
|---|---|---|
| ⭐ 最新专题 ⭐ | DSGE, 因果推断, 空间计量等 | |
| ⭕ Stata数据清洗 | 游万海 | 直播, 2 小时，已上线 |
| 研究设计 | 连玉君 | 我的特斯拉-实证研究设计，-幻灯片- |
| 面板模型 | 连玉君 | 动态面板模型，-幻灯片- |
| 面板模型 | 连玉君 | 直击面板数据模型 [免费公开课，2小时] |
Note: 部分课程的资料，PPT 等可以前往 连享会-直播课 主页查看，下载。
关于我们
Stata连享会 由中山大学连玉君老师团队创办，定期分享实证分析经验。直播间 有很多视频课程，可以随时观看。 连享会-主页 和 知乎专栏，300+ 推文，实证分析不再抓狂。 公众号推文分类： 计量专题 | 分类推文 | 资源工具。推文分成 内生性 | 空间计量 | 时序面板 | 结果输出 | 交乘调节 五类，主流方法介绍一目了然：DID, RDD, IV, GMM, FE, Probit 等。 公众号关键词搜索/回复 功能已经上线。大家可以在公众号左下角点击键盘图标，输入简要关键词，以便快速呈现历史推文，获取工具软件和数据下载。常见关键词： 课程, 直播, 视频, 客服, 模型设定, 研究设计, stata, plus, 绘图, 编程, 面板, 论文重现, 可视化, RDD, DID, PSM, 合成控制法等
连享会小程序：扫一扫，看推文，看视频……
扫码加入连享会微信群，提问交流更方便
✏ 连享会学习群-常见问题解答汇总：
✨ https://gitee.com/arlionn/WD
热门资讯
