R语言：断点回归（RD）学习手册（包含rdrobust命令详解、设计前提条件内生分组等显著性检验、全套标准动作）
R语言：断点回归（RD）学习手册（包含rdrobust命令详解、设计前提条件内生分组等显著性检验、全套标准动作）
断点回归由Thistlewaite and Campbell(1960)首次使用，但直到1990年代末才引起经济学家的重视。
Thistlethwaite、Campbell于1960年首次提出使用断点回归设计研究处理效应， 在该文中他们的目的是研究奖学金对于未来学业的影响, 学生是否获得奖学金取决于考试的分数。由于奖学金由学习成绩决定，故成绩刚好达到获奖标准与差一点达到的学生具有可比性。如果考试分数大于获奖标准分数, 则进入处理组；如果考试分数小于获奖标准分数, 则进入控制组。因此处理变量在获奖标准分数处形成了一个断点， 该研究设计的主要思想是可以利用靠近这一断点附近的样本来有效估计处理效应。
Angrist and Lavy(1999)在研究班级规模对成绩的影响时，利用以色列教育系统的一项制度进行断点回归；该制度限定班级规模的上限为40名学生，一旦超过40名学生（比如41名学生），则该班级被一分为二。
此后30年, 该方法并未引起学术界的重视，直到1990年以后, 断点回归设计开始被应用于各种领域，并且近年来成为因果分析和政策评估领域最重要的研究方法。
Hahn et al(2001)提供了断点回归在计量经济学理论基础。目前，断点回归在教育经济学、劳动经济学、健康经济学、政治经济学以及区域经济学的应用仍方兴未艾。参见Imbens and Lemieux(2008)，Van Der Klaauw(2008)以及Lee and Lemieux(2010)的文献综述。
断点回归设计是一种准自然实验, 其基本思想是存在一个连续变量, 该变量能决定个体在某一临界点两侧接受政策干预的概率, 由于X在该临界点两侧是连续的，因此个体针对X的取值落入该临界点任意一侧是随机发生的, 即不存在人为操控使得个体落入某一侧的概率更大, 则在临界值附近构成了一个准自然实验。一般将该连续变量X称为分组变量 (assignment variable) 。
一．断点回归操作动作建议
在进行断点回归（RD）设计时，一般有如下步骤：
1、参考变量分布连续性检验/检验内生分组
这里检验内生分组，即主要检验配置变量，其实就是RD中个体是否将自行进入断点两侧，决定是否进入实验的，并是否存在某种跳跃性的变化。如果存在内生分组，个体将自行进入实验，导致在断点两侧的分布不均匀，这样分组变量x的密度函数f（x）在x=c处不连续，出现左右极限不相等的情况。
McCrary（2008）提出了一种核密度函数的检验方法（命令是DCdensity，介绍见下述操作），将参考变量划分成不同的区间并计算各区间中的个体数量，如果个体能够操纵参考变量，我们将能观测到断点左右个体数量有较大差别，比如很多个体通过操纵到了断点的右侧，那么，在断点右侧的区间中个体数量可能将大大超过断点左侧区间中个体的数量，利用带宽选择和曲线拟合方法， 可以检验在断点处c是否存在跳跃 。
2、检查为精确断点回归还是模糊断点回归分析
检验处理变量是否完全由“某连续变量是否超过某一断点”所决定，如果个体被处理的概率从0跳跃为1，即为精确断点回归，如果个体被处理的概率从 a跳跃为 b，0<a<b<1，则为模糊断点回归。
3、图形分析
画出结果变量与参考变量之间的关系图，如果是模糊断点，再画出原因变量与参考变量的关系图，呈现结果变量和原因变量在断点处行为，为断点回归设计提供理论支撑。
4、检验结果对不同带宽、不同多项式次数的稳健性
设置不同带宽，通过选择最优带宽，再检验并选择相对应的模型。stata断点回归命令有相关的操作选项。另外还有图形选择（在最优带宽处画线），可以考虑加协变量进行选择。
5、检验其他影响结果变量的因素（协变量），在断点处是否存在跳跃
检验协变量在断点处是否存在跳跃，若是存在跳跃，说明该协变量的条件密度函数在断点处不是连续的，需要剔除。若将存在跳跃的协变量剔除。则需要重新选择最优带宽再重新进行断点回归分析。
6、显著性检验
模型估计完成后，可以进行下列模型设定检验，以判断估计结果的稳健性（见赵西亮编著的《基本有用的计量经济学》）
（1）协变量连续性检验，也称为伪结果检验（ pseudo outcome）。以协变量 作为伪结果，利用与前面相同的方法，检验相应的RDD估计量是否显著，如果 显著说明这些协变量不符合连续性假设，上文的RDD估计量可能存在问题。
（2）参考变量分布连续性检验，如果参考变量分布连续，意味着在断点处个体没有精确操纵参考变量的能力，局部随机化假设成立，从而保证断点附近左右样本能够代表断点处的总体。（此处与检验内生分组一致）
（3）伪断点检验（ pseudo cutoff point）。在参考变量的其他位置，比如断点 左右两侧中点位置作为伪断点，利用同样的方法估计RDD估计量，我们知道在 伪断点干预效应为零，如果发现伪断点的RDD估计量不为零，则说明我们的RDD设计可能有问题，可能混杂了其他未观测因素的影响，得到的因果效应可能是由其他未观测混杂的跳跃造成的，而不完全是干预的影响
（4）带宽选择的敏感性检验。选择不同的带宽对RDD估计量进行重新估 计，检验估计结果是否有较大的变量，如果差异较大，尤其是影响方向有变化说明RDD设计可能有问题。
上述显著性检验其实在前面进行分析时候已经部分有所提及需要进行检验的。
二 、断点回归命令
在Stata中，断点回归的基本命令是rd，另外，还有一些其他命令，例如rdrobust、rdlocrand、rddensity等等。
本文主要介绍在R语言中的rdrobust。rdrobust有两个配套命令:rdbwselect用于带宽选择，rdplot用于RD绘图(详细信息请参见Calonico、Cattaneo和Titiunik [2015a])。
1、rdbwselect
下载安装方法为：
install.packages('rdrobust')
library(rdrobust)
rdrobust语法格式为：
rdbwselect(y, x, c = NULL, fuzzy = NULL,
deriv = NULL, p = NULL, q = NULL,
covs = NULL, covs_drop = TRUE,
kernel = "tri", weights = NULL, bwselect = "mserd",
vce = "nn", cluster = NULL, nnmatch = 3,
scaleregul = 1, sharpbw = FALSE,
all = NULL, subset = NULL,
masspoints = "adjust", bwcheck = NULL,
bwrestrict = TRUE, stdvars = FALSE)
选项含义为：
y表示因变量。
x表示驱动变量(又称使动变量)。
c指定RD中断点位置;默认值是c = 0。
fuzzy指定用于实现模糊RD估计的处理状态变量，默认是Sharp RD设计，因此不使用此选项。
cov表示指定用于估计和推断的附加协变量。
kernel是用来构造局部多项式估计器的核函数。选项有三角形(默认选项)、epanechnikov和uniform。
weights是用于估计过程的可选加权的变量。单位权重乘以核函数。
bwselect指定要使用的带宽选择过程。
2、rdrobust
语法含义为：
rdrobust(y, x, c = NULL, fuzzy = NULL,
deriv = NULL, p = NULL, q = NULL,
h = NULL, b = NULL, rho = NULL, covs = NULL, covs_drop = TRUE,
kernel = "tri", weights = NULL, bwselect = "mserd",
vce = "nn", cluster = NULL,
nnmatch = 3, level = 95, scalepar = 1, scaleregul = 1,
sharpbw = FALSE, all = NULL, subset = NULL,
masspoints = "adjust", bwcheck = NULL,
bwrestrict = TRUE, stdvars = FALSE)
选项含义为：
y表示因变量。
x表示驱动变量(又称使动变量)。
c指定RD中断点位置;默认值是c = 0。
fuzzy指定用于实现模糊RD估计的处理状态变量，默认是Sharp RD设计，因此不使用此选项。
covs表示指定用于估计和推断的附加协变量。
p指定用于构造点估计器的局部多项式的阶数;默认值是p = 1(本地线性回归)。
q指定用于构造偏差校正的局部多项式的顺序;默认值是q = 2(局部二次回归)。
h指定用于构造RD点估计器的主要带宽。如果没有指定，带宽h将由同伴命令rdbwselect计算。如果指定了两个带宽，第一个带宽用于低于截止值的数据，第二个带宽用于高于截止值的数据。
b指定用于构造偏差校正估计器的偏差带宽。如果没有指定，带宽b将由同伴命令rdbwselect计算。如果指定了两个带宽，第一个带宽用于低于截止值的数据，第二个带宽用于高于截止值的数据。
3、rdplot
语法含义为
rdplot(y, x, c = 0, p = 4, nbins = NULL, binselect = "esmv",
scale = NULL, kernel = "uni", weights = NULL, h = NULL,
covs = NULL, covs_eval = 0, covs_drop = TRUE,
support = NULL, subset = NULL,
hide = FALSE, ci = NULL, shade = FALSE, title = NULL,
x.label = NULL, y.label = NULL, x.lim = NULL, y.lim = NULL,
col.dots = NULL, col.lines = NULL)
选项含义为
rdplot函数可用于RD结果的可视化，其语法与rdrobust或者rdrobust相同
y表示因变量。
x表示驱动变量(又称使动变量)。
c指定RD中断点位置;默认值是c = 0。
p指定用于近似控制单元和处理单元的总体条件平均函数的全局多项式的阶数;默认值是p = 4。
binselect指定选择容器数量的过程。
三、案例
rdrobust采用Calonico, Cattaneo和Titiunik (2014a)、Calonico, Cattaneo和Farrell(2018)、Calonico, Cattaneo, Farrell和Titiunik(2019)和Calonico, Cattaneo和Farrell(2020)开发的鲁棒偏差校正置信区间和推理程序，实现了局部多项式回归不连续(RD)点估计。它还计算可供选择的估计和推理程序，在文献中。
配套命令有:rdbwselect用于数据驱动带宽选择，rdplot用于数据驱动RD图(详见Calonico, Cattaneo and Titiunik (2015a))。
这个命令的详细介绍在Calonico, Cattaneo和Titiunik (2015b)和Calonico, Cattaneo, Farrell和Titiunik(2017)中给出。配套的Stata包在Calonico, Cattaneo和Titiunik (2014b)中有描述。
有关更多细节，以及用于分析RD设计的相关Stata和R包，请访问https://sites.google.com/site/rdpackages/
Cattaneo, Frandsen和Titiunik构建的数据集(2015)，其中包括1914-2010年期间美国参议院在任优势的衡量。
该数据包含以下两个变量的1390个观察值的数据框架。
其包含两个变量vote和margin，vote表示某次选举民主党在州参议院的席位占比，margin表示上次选举中获得相同参议院席位的边际收益，其中大于0表示民主党胜出，反之则为失败。
将vote作为被解释变量，margin作为解释变量，即可研究民主党赢得参议院席位对于在下次选举中获得相同席位的影响。具体操作代码如下
1导入数据
data(rdrobust_RDsenate)
2最优带宽
rdbwselect(y = rdrobust_RDsenate$vote, x = rdrobust_RDsenate$margin, all = TRUE)
summary()
结果为：
> rdbwselect(y = rdrobust_RDsenate$vote, x = rdrobust_RDsenate$margin, all = TRUE)
Call: rdbwselect
Number of Obs. 1297
BW type All
Kernel Triangular
VCE method NN
Number of Obs. 595 702
Order est. (p) 1 1
Order bias (q) 2 2
Unique Obs. 595 665
=======================================================
BW est. (h) BW bias (b)
Left of c Right of c Left of c Right of c
=======================================================
mserd 17.754 17.754 28.028 28.028
msetwo 16.170 18.126 27.104 29.344
msesum 18.365 18.365 31.319 31.319
msecomb1 17.754 17.754 28.028 28.028
msecomb2 17.754 18.126 28.028 29.344
cerrd 12.407 12.407 28.028 28.028
certwo 11.299 12.667 27.104 29.344
cersum 12.834 12.834 31.319 31.319
cercomb1 12.407 12.407 28.028 28.028
cercomb2 12.407 12.667 28.028 29.344
=======================================================
3 参数估计
rdrobust(y = rdrobust_RDsenate$vote, x = rdrobust_RDsenate$margin) %>% summary()
Call: rdrobust
Number of Obs. 1297
BW type mserd
Kernel Triangular
VCE method NN
Number of Obs. 595 702
Eff. Number of Obs. 360 323
Order est. (p) 1 1
Order bias (q) 2 2
BW est. (h) 17.754 17.754
BW bias (b) 28.028 28.028
rho (h/b) 0.633 0.633
Unique Obs. 595 665
=============================================================================
Method Coef. Std. Err. z P>|z| [ 95% C.I. ]
=============================================================================
Conventional 7.414 1.459 5.083 0.000 [4.555 , 10.273]
Bias-Corrected 7.507 1.459 5.146 0.000 [4.647 , 10.366]
Robust 7.507 1.741 4.311 0.000 [4.094 , 10.919]
=============================================================================
4 绘图
rdplot(rdrobust_RDsenate$vote, rdrobust_RDsenate$margin, p = 3, x.label = 'Margin', y.label = 'Vote')
四、稳健性检验
下面我们结合rdrobust这个命令来介绍如何进行断点回归的一个稳健性检验。
也就是rd的有效性检验，首先进行RD断点回归设计的时候呢，需要满足一些先决条件。
1、局部平滑性假设
除了结果变量，其他协变量，在断点处也是没有出现跳跃性的这样一个现象，满足平滑性，也就是连续性，我们可以假设将协变量作为结果变量，利用rdrobust函数来进行检验，协变量应该在统计学上是不显著的。
上述假设也称为协变量连续性检验，也称为伪结果检验（ pseudo outcome）。以协变量 作为伪结果，利用与前面相同的方法，检验相应的RDD估计量是否显著，如果 显著说明这些协变量不符合连续性假设，上文的RDD估计量可能存在问题。
具体代码如下
rdrobust(y = other variables, x = assignment variable) %>% summary()
2、驱动变量的有效性
也就是参考变量分布连续性检验，如果参考变量分布连续，意味着在断点处个体没有精确操纵参考变量的能力，局部随机化假设成立，从而保证断点附近左右样本能够代表断点处的总体。（此处与检验内生分组一致）
该假设可通过rddensity中的对应函数进行McCrary检验，具体代码如下
library(rdrobust)
library(rddensity)
data(rdrobust_RDsenate)
rdplotdensity(rddensity(rdrobust_RDsenate$margin), X = rdrobust_RDsenate$margin)
3、RDD中的稳健性检验
对于RDD中的稳健性检验，一般主要包括对伪断点检验、样本选择和带宽选择进行分析。
3.1、伪断点检验
伪断点检验（ pseudo cutoff point）。在参考变量的其他位置，比如断点 左右两侧中点位置作为伪断点，利用同样的方法估计RDD估计量，我们知道在 伪断点干预效应为零，如果发现伪断点的RDD估计量不为零，则说明我们的RDD设计可能有问题，可能混杂了其他未观测因素的影响，得到的因果效应可能是由其他未观测混杂的跳跃造成的，而不完全是干预的影响
具体可使用rdrobust进行检验，预期结果为系数不再显著。具体代码如下
## 除c以外的参数应当与benchmark一致
sen_cut_1 <- rdrobust(rdrobust_RDsenate$vote, rdrobust_RDsenate$margin, c = -2, p = 1, kernel = 'triangular', bwselect = 'msetwo')
sen_cut_2 <- rdrobust(rdrobust_RDsenate$vote, rdrobust_RDsenate$margin, c = -1, p = 1, kernel = 'triangular', bwselect = 'msetwo')
sen_cut_3 <- rdrobust(rdrobust_RDsenate$vote, rdrobust_RDsenate$margin, c = 0, p = 1, kernel = 'triangular', bwselect = 'msetwo')
sen_cut_4 <- rdrobust(rdrobust_RDsenate$vote, rdrobust_RDsenate$margin, c = 1, p = 1, kernel = 'triangular', bwselect = 'msetwo')
sen_cut_5 <- rdrobust(rdrobust_RDsenate$vote, rdrobust_RDsenate$margin, c = 2, p = 1, kernel = 'triangular', bwselect = 'msetwo')
## 创建含有coef和se的数据帧
coef <- c(sen_cut_1$coef[1], sen_cut_2$coef[1], sen_cut_3$coef[1], sen_cut_4$coef[1], sen_cut_5$coef[1])
se <- c(sen_cut_1$se[1], sen_cut_2$se[1], sen_cut_3$se[1], sen_cut_4$se[1], sen_cut_5$se[1])
robust_cutpoint <- data.frame(position = c(-2:2), coef, se)
## 利用ggplot2画图
(photo2 <- ggplot(data = robust_cutpoint, aes(x = position)) + geom_point(aes(y = coef), shape = 15)) + geom_errorbar(aes(ymin=coef-1.96*se, ymax=coef+1.96*se), width=.03) + geom_hline(yintercept = 0, colour = "red")
3.2 样本选择的稳健性检验
我们也可以随机去除断点附近的样本来证明结果的稳健性。若去除一定数量的样本后结果依然显著，则结果稳健。具体代码如下
## 生成去除30%、40%和50%样本个体的子样本
rdrobust_RD_hole_1 <- subset(rdrobust_RDsenate, abs(margin) > 0.3)
rdrobust_RD_hole_2 <- subset(rdrobust_RDsenate, abs(margin) > 0.4)
rdrobust_RD_hole_3 <- subset(rdrobust_RDsenate, abs(margin) > 0.5)
sen_hole_1 <- rdrobust(rdrobust_RD_hole_1$vote, rdrobust_RD_hole_1$margin, c = 0, p = 1, kernel = 'triangular', bwselect = 'msetwo')
sen_hole_2 <- rdrobust(rdrobust_RD_hole_2$vote, rdrobust_RD_hole_2$margin, c = 0, p = 1, kernel = 'triangular', bwselect = 'msetwo')
sen_hole_3 <- rdrobust(rdrobust_RD_hole_3$vote, rdrobust_RD_hole_3$margin, c = 0, p = 1, kernel = 'triangular', bwselect = 'msetwo')
3.3 带宽选择稳健性检验
带宽选择的敏感性检验。选择不同的带宽对RDD估计量进行重新估 计，检验估计结果是否有较大的变量，如果差异较大，尤其是影响方向有变化说明RDD设计可能有问题。
一般来说，可以尝试使用rdbwselect计算得出的若干种最优带宽。若回归系数依然显著，则可证明结果稳健性。具体代码如下
sen_bandwidth_1 <- rdrobust(rdrobust_RD_hole_1$vote, rdrobust_RD_hole_1$margin, c = 0, p = 1, kernel = 'triangular', bwselect = 'msetwo')
sen_bandwidth_2 <- rdrobust(rdrobust_RD_hole_2$vote, rdrobust_RD_hole_2$margin, c = 0, p = 1, kernel = 'triangular', bwselect = 'cerrd')
sen_bandwitch_3 <- rdrobust(rdrobust_RD_hole_3$vote, rdrobust_RD_hole_3$margin, c = 0, p = 1, kernel = 'triangular', bwselect = 'cercomb')