RDD: 断点回归的非参数估计及Stata实现
作者：崔颖（中央财经大学）
Source: Non-Parametric Regression Discontinuity (Francis, 2013)
2019暑期Stata现场班 (连玉君+刘瑞明主讲)
特别说明
文中包含的链接在微信中无法生效。请点击本文底部左下角的【阅读原文】，转入本文【简书版】。
本篇推文介绍Stata方便实现断点回归 (Regression Discontinuity, RD) 的实用命令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) 的带宽及其一半与两倍的带宽。
我们可以将不同带宽的估计结果绘制在一张图上。
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
关于我们
【Stata 连享会(公众号：StataChina)】由中山大学连玉君老师团队创办，旨在定期与大家分享 Stata 应用的各种经验和技巧。
公众号推文同步发布于 CSDN-Stata连享会 、简书-Stata连享会 和 知乎-连玉君Stata专栏。可以在上述网站中搜索关键词
Stata或Stata连享会后关注我们。点击推文底部【阅读原文】可以查看推文中的链接并下载相关资料。
Stata连享会 精彩推文1 || 精彩推文2
联系我们
欢迎赐稿： 欢迎将您的文章或笔记投稿至
Stata连享会(公众号: StataChina)，我们会保留您的署名；录用稿件达五篇以上，即可免费获得 Stata 现场培训 (初级或高级选其一) 资格。意见和资料： 欢迎您的宝贵意见，您也可以来信索取推文中提及的程序和数据。
招募英才： 欢迎加入我们的团队，一起学习 Stata。合作编辑或撰写稿件五篇以上，即可免费获得 Stata 现场培训 (初级或高级选其一) 资格。
联系邮件： StataChina@163.com
往期精彩推文
Go to "Discover" > "Top Stories" > "Wow"