RDD断点回归设计的详细步骤和代码指南

编辑整理: @计量经济圈; 来源:熊彼特的厨房 公众号

RDD断点回归设计越来越来越流行，在因果识别中的地位很高，因为是在一个断点处的两端进行比较，那么这个treatment效应能够较好的识别，下面这篇文章将会领着计量经济圈的圈友去认识RDD并且知道操作的步骤和Stata代码。经典论文完整版do文件，请进入计量社群领取，我们统一在社群发布(文后)。

1.实证研究中的规定动作

第1步检查配置变量（assignment variable，又叫running variable、forcing variable）是否被操纵。

这里的配置变量，其实就是RD中决定是否进入实验的分数（Score），是否被操纵的意思就是，是否存在某种跳跃性的变化。在实际操作中有两种方式来检验，一是画出配置变量的分布图。最直接的方法，是使用一定数量的箱体（bin），画出配置变量的历史直方图（histogrm）。为了观察出分布的总体形状，箱体的宽度要尽量小。频数（frequencies）在箱体间的跳跃式变化，能就断点处的跳跃是否正常给我们一些启发。从这个角度来说，最好利用核密度估计做出一个光滑的函数曲线。二是利用McCrary（2008）的核密度函数检验。（命令是DCdensity，介绍见陈强编著的《高级计量经济学及Stata应用》（第二版）第569页）, Frandsen (2013)提出了一种新的检验方法，但目前被使用的并不多。

第2步画因变量均值对配置变量的散点图，并选择带宽（bandwidth selection）。

首先，挑选出一定数目的箱体，求因变量在每个箱体内的均值，画出均值对箱体中间点的散点图。一定要画每个箱体平均值的图。如果直接画原始数据的散点图，那么噪音太大，看不出潜在函数的形状。不要画非参数估计的连续统，因为这个方法自然地倾向于给出存在断点的印象，尽管总体中本来不存在这样的断点。

然后，选择第三步骤中需要的带宽。Lee和Lemieux（2010）介绍了两种确定最优带宽的方法：拇指规则法（rule of thumb）和交叉验证法（CV）。还有另外两种比较受关注的方法：IK法和CCT法。IK法以Imbens和Kalyanaraman两个人命名，对应着论文Imbens和Kalyanaraman（2012）。这篇论文发表在Review of Economic Studies，Lee和Lemieux（2010）文中提到过此文2009年的NBER工作论文版。CCT法以Calonico、Cattaneo和Titiunik三个人命名，对应着论文Calonico、Cattaneo和Titiunik（2014a）。用非参数法做断点回归估计时的stata命令rd，就是用IK发确定最优带宽。stata命令rdrobust、rdbwselect，提供CV、IK、CCT三种不同的最优带宽计算方法选项。但是实际上rdrobust中已经更新了IK带宽选择函数，更新的算法与IK算法的区别有待考证，后续会补充。实际操作中一般是两种算法都会采纳，并汇报参数估计对带宽选择是不敏感的。

第3步估计，又分为参数估计和非（半）参数估计两种方案。

在Sharp RD情形，参数法将Y在每个箱体内的均值作为因变量，用处理变量、配置变量的多次项作为自变量，在断点两边分别跑回归，得到断点左右两边因变量的拟合值，两个拟合值的差值便是我们想估计的实验对因变量的因果效应。将这些拟合值画在第2步的图中，并用光滑的曲线连接起来。在推文人读过的RD论文中，多次项一般都使用1到4次项，但没有论文解释为什么只用到4次项。半参数的方法便是用非参数估计的方法替代断点两边估计因变量的拟合。对于Fuzzy的情形，参数估计意味着将配置变量（score）以及配置变量与是否超过断点的乘积(score*above_cutoff)作为实验变量的工具变量来进行两阶段最小二乘估计，实际应用中往往联合使用score, score*above_cutoff,score^2,score^2* above_cutoff作为估计的工具变量，见第二部分的例子。非参数估计的做法是，利用核密度函数局部现性回归来代替2SLS里面一般线性回归，rdrobust命令可以直接实现这种估计。

第4步检验前定变量在断点处是否跳跃，这一步的目的是证明不存在跳跃，否则就麻烦了。

前定变量指的是那些在实验之前已经确定的变量，例如，发生在2008年的实验，那些2007年的观测值便是前定的，理论上这些变量是不应该在断点出跳跃的。此步和第1步是RD方法的适用性检验。此步的检验包括两项内容：1.像前三步那样画前定变量的图。无论参数还是非参数，RD研究都要大把的图！这些图在正式发表的论文中都必不可少！原文中说了这么句话：用RD做的论文，如果缺乏相关的图，十有八九是因为图显示的结果不好，作者故意不报告。2.将前定变量作为因变量，将常数项、处理变量、配置变量多次项、处理变量和配置变量多次项的交互项作为自变量，跑回归。一个前定变量有一个回归，看所有回归中处理变量的系数估计是否都为。检验这种跨方程的假设，需要用似不相关回归（Seemingly Unrelated Regression, SUR）（命令是sureg，用法见陈强编著的《高级计量经济学及Stata应用》（第二版）第471-474页）。在推文人读过的RD实证论文中（尤其是AER2015-2016年所有用RD做的论文中），均没用SUR，只是简单的看每个回归中处理变量的系数估计均为。

第5步检验结果对不同带宽、不同多项式次数的稳健性。

尝试的其它带宽，一般是最优带宽的一半和两倍。挑选多项式的最优次数，可用赤池信息准则（Akaike's Information Criterion，AIC）。在我们尝试的包含配置变量1次方、2次方、……N次方的众多方程中，AIC取值最小的那个就是我们想要的。实操时，试到多少次为好？Gelman和Imbens（2014）的NBER工作论文说，试到N次的做法要不得，最多只能搞到2次。

第6步检验结果对加入前定变量的稳健性。

如上所述，如果不能操控配置变量的假设成立，那么无论前定变量与因变量的相关性有多高，模型中加入前定变量都不应该影响处理效应的估计结果。如果加入前定变量导致处理效应的估计结果变化较大，那么配置变量可能存在排序现象，前定变量在断点处也很可能存在跳跃。实操时在确定多项式的次数后，直接在回归方程中加入前定变量。如果这导致处理效应估计值大幅变化或者导致标准误大幅增加，那么可能意味着函数中多项式的次数不正确。另外一个检验是残差化，看相同次数的多项式模型对残差的拟合好不好。

Pinotti, Paolo(2017) AER论文的实际操作

Replication大法好！否则就会“读了那么多文献还是做不好研究。这篇论文的研究问题是移民获得合法身份后的犯罪是否会减少。其中，实验变量是是否获得合法身份，因变量是移民申请人在申请合法身份后的第一年（2008年）是否有犯罪记录，这里面的选择偏误是合法身份并不是随机发放给移民申请人的，那些预期犯罪更少的移民更有机会（有雇主帮助申请）和动力（花更多的时间和精力去准备申请）来申请合法身份，直接对比犯罪率会夸大合法身份的作用（使负向系数更小）。为了克服这个选择偏误，作者利用了意大利的自然实验，意大利的移民身份申请是先到先得的，也就是在系统开放后，申请时间越晚中签率会越低，作者发现申请递交时间timing上有一个断点，当申请人晚于这个断点提交申请书时，会导致其中签率跳跃式下滑，但不至于完全为零，于是作者找到了一个Fuzzy断点情景。下面我们看看作者是怎么操作的。其实这个情景很像上海的机动车抽签系统。

2.1配置与处置变量的散点图

在实际操作中，作者Pinotti在描述政策背景的时候直接汇报了配置与处置变量的散点图（5min作为箱体），如下图

然后作者使用Andrews(1993)检验的方法检验存在结构性断点，但是作者没有在正文与附录中汇报该检验的结果。

2.2确定断点

作者针对每一个“摇号点”，用Andrew(1993)的方案找出“the most likely break points”，同样作者没有汇报cutoff point的详细数据以及寻找过程。仅仅在附件中汇报了下图。

2.3全局2SLS估计

作者将样本人为限制在cutoff point附近半个小时以内样本。

（1）首先汇报了因变量（2008年犯罪率）与前定变量（2007年犯罪率）与配置变量的散点图及其置信区间。

上图的实际实现分为两步，第一步是用是否大于cutoff point的dummy与score的四次项多项式回归，得出不同分数的拟合值与置信区间，第二步是画出散点图与拟合曲线图。以2007年为例，stata实现为：

bys bin: egen mean=mean(serious07)

reg serious07 ontime mindelay mindelay2 ontimexmindelay ontimexmindelay2, robust

predict fit

predict fitsd, stdp

gen upfit=fit+1.645*fitsd

gen downfit=fit-1.645*fitsd

preserve

twoway (rarea upfit downfit mindelay, sort fcolor(gs12) lcolor(gs12)) ///

(line fit mindelay if mindelay0, sort lcolor(red) lwidth(thick)) (scatter mean midbin, msize(large) mcolor(black) msymbol(circle_hollow)), ///

ytitle("") xtitle("Timing of the application, X (cutoff: X=0)") xline(0, lcolor(black)) legend(off) xlabel(-30(10)30) title("2007, all applicants")

graph copy all2007, replace

restore

drop *fit* mean

（2）然后,作者在一张表中汇报了2sls以及前定变量的ols结果

通过上表，我们可以得出在前定变量方面，断点两边的差异是不显著，在因变量方面显著，而且显著性来自type-A样本。在回归中作者使用了四次项多项式，而且只用了ontime这一个工具变量。以2008年为例，其实stata现方式为:

*** reduced form, 2008 ***

reg serious08 ontime mindelay mindelay2 ontimexmindelay ontimexmindelay2, robust

reg serious08 ontime mindelay mindelay2 ontimexmindelay ontimexmindelay2 if type=="A", robust

reg serious08 ontime mindelay mindelay2 ontimexmindelay ontimexmindelay2 if type=="B", robust

xtreg serious08 ontime mindelay mindelay2 ontimexmindelay ontimexmindelay2 age age2, robust i(lotterycode) fe cluster(lotterycode)

xtreg serious08 ontime mindelay mindelay2 ontimexmindelay ontimexmindelay2 age age2 if type=="A", robust i(lotterycode) fe cluster(lotterycode)

xtreg serious08 ontime mindelay mindelay2 ontimexmindelay ontimexmindelay2 age age2 if type=="B", robust i(lotterycode) fe cluster(lotterycode)

*** 2SLS, 2008 ***

ivreg2 serious08 (permit2007=ontime) mindelay mindelay2 ontimexmindelay ontimexmindelay2, robust first

ivreg2 serious08 (permit2007=ontime) mindelay mindelay2 ontimexmindelay ontimexmindelay2 if type=="A", robust first

ivreg2 serious08 (permit2007=ontime) mindelay mindelay2 ontimexmindelay ontimexmindelay2 if type=="B", robust first

xtivreg2 serious08 (permit2007=ontime) mindelay mindelay2 ontimexmindelay ontimexmindelay2 age age2, robust first i(lotterycode) fe cluster(lotterycode)

xtivreg2 serious08 (permit2007=ontime) mindelay mindelay2 ontimexmindelay ontimexmindelay2 age age2 if type=="A", robust first i(lotterycode) fe cluster(lotterycode)

xtivreg2 serious08 (permit2007=ontime) mindelay mindelay2 ontimexmindelay ontimexmindelay2 age age2 if type=="B", robust first i(lotterycode) fe cluster(lotterycode)

2.4其他前定变量的稳健性

作者分别汇报了Age,来源国，来源国收入水平等维度上的配置变量散点图来显示cutoff point两边的无差异性，做法与图三相同。

2.5稳健性检验

（1）局部现性2SLS估计

为了体现论文的稳健性，作者接着使用了局部线性估计，结果如下：

在实际操作中，作者汇报了IK2012与CCT2014两种带宽选择方案的结果，其中IK的样本大约为40%左右，CCT为30%。

表4的实现方式如下，注意，作者并没有使用Covariates,我的理解是实际操作中Cov的影响很小（毕竟在断点附近COV几乎没有差异），而且离散型是的Cov会导致结果不收敛。

*** reduced form, 2008 ***

rdrobust serious08 mindelay, bwselect(IK)

rdrobust serious08 mindelay if type=="A", bwselect(IK)*新的rdrobust命令中bwselect算法已经更新

rdrobust serious08 mindelay if type=="B", bwselect(IK)

rdrobust serious08 mindelay

rdrobust serious08 mindelay if type=="A"

rdrobust serious08 mindelay if type=="B"

*** 2SLS, 2008 ***

rdrobust serious08 mindelay, fuzzy(permit2007) bwselect(IK)

rdrobust serious08 mindelay if type=="A", fuzzy(permit2007) bwselect(IK)

rdrobust serious08 mindelay if type=="B", fuzzy(permit2007) bwselect(IK)

rdrobust serious08 mindelay, fuzzy(permit2007)

rdrobust serious08 mindelay if type=="A", fuzzy(permit2007)

rdrobust serious08 mindelay if type=="B", fuzzy(permit2007)

（2）作者使用不同的多项式次数与带宽选择进行估计，汇报结果如下

实际实现为：

*** generate polynomials grade 3-6 ***

forvalues i=3/6 {

gen mindelay$_i=mindelay^$_i

gen ontimexmindelay$_i=ontime*mindelay$_i

}

*** PARAMETRIC ESTIMATES (top graphs) ***

foreach t in "A" "B" {

preserve

keep if type=="$_t"

qui xtreg serious08 ontime age age2, robust i(lotterycode) fe

display "type $_t| polynomial degree 0| coeff. "_b[ontime] "| std.err. " _se[ontime]

qui xtreg serious08 ontime mindelay ontimexmindelay age age2, robust i(lotterycode) fe

display "type $_t| polynomial degree 1| coeff. "_b[ontime] "| std.err. " _se[ontime]

qui xtreg serious08 ontime mindelay mindelay2 ontimexmindelay ontimexmindelay2 age age2, robust i(lotterycode) fe

display "type $_t| polynomial degree 2| coeff. "_b[ontime] "| std.err. " _se[ontime]

qui xtreg serious08 ontime mindelay mindelay2 mindelay3 ontimexmindelay ontimexmindelay2 ontimexmindelay3 age age2, robust i(lotterycode) fe

display "type $_t| polynomial degree 3| coeff. "_b[ontime] "| std.err. " _se[ontime]

qui xtreg serious08 ontime mindelay mindelay2 mindelay3 mindelay4 ontimexmindelay ontimexmindelay2 ontimexmindelay3 ontimexmindelay4 age age2, robust i(lotterycode) fe

display "type $_t| polynomial degree 4| coeff. "_b[ontime] "| std.err. " _se[ontime]

qui xtreg serious08 ontime mindelay mindelay2 mindelay3 mindelay4 mindelay5 ontimexmindelay ontimexmindelay2 ontimexmindelay3 ontimexmindelay4 ontimexmindelay5 age age2, robust i(lotterycode) fe

display "type $_t| polynomial degree 5| coeff. "_b[ontime] "| std.err. " _se[ontime]

qui xtreg serious08 ontime mindelay mindelay2 mindelay3 mindelay4 mindelay5 mindelay6 ontimexmindelay ontimexmindelay2 ontimexmindelay3 ontimexmindelay4 ontimexmindelay5 ontimexmindelay6 age age2, robust i(lotterycode) fe

display "type $_t| polynomial degree 6| coeff. "_b[ontime] "| std.err. " _se[ontime]

restore

}

*** NONPARAMETRIC ESTIMATES (bottom graphs) ***

foreach t in "A" "B" {

preserve

keep if type=="$_t"

forvalues h=1/30 {

qui rdrobust serious08 mindelay if type=="A", fuzzy(permit2007) h($_h)

display "type $_t| bandwidth $_h | `e(tau_F_cl)' | `e(se_F_cl)' | `e(N)'"

}

restore

}

作者同时在附录中汇报了前定变量的回归结果：

2.6安慰剂检验

（实际上就是P值的可视化，算是一个cutoff point的稳健性检验）

帮点击一下下面的小广告，谢谢支持！