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

欢迎投稿(荐稿)计量经济圈，计量相关都行

邮箱：econometrics666@sina.cn

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

RDD断点回归设计越来越来越流行，在因果识别中的地位很高，因为是在一个断点处的两端进行比较，那么这个treatment效应能够较好的识别，下面这篇文章将会领着计量经济圈的圈友去认识RDD并且知道操作的步骤和Stata代码。想要进一步了解RDD这种方法以及其他因果推断技巧的，进入计量经济圈社群交流(文后)。

目前常见的用于RDD计算的R包分别有rdd, rdrobust, rddtools。Thoemmes等 (2017)的论文比较了这三个包的异同。

Thoemmes（2017）从实践的角度，进一步将RD的步骤分为三类：Assumption Checks，Estimation, Sensitivity Checks。我认为这里的简要归纳要比Lee和Lemieux(2010)更清楚简洁。

Assumption Checks实际上对应了RD的三个前提假设,分别是配置变量存在断点，对应着McCrary（2008）的断点检验；实验效应发生在我们假定的断点处，而不是其他断点处，这对应着安慰剂检验（Placebo Test），在上篇第二部分有展示；实验变量仅仅对我们关心的结果变量有影响，但是对其他前定变量无影响，其实就是检验断点附近前定变量的连续性。

Estimation是对因果效应的估计，又分为参数法（Global Regression）和非参数法（Local Linear Regression），后者对应着不同带宽的选择。

Sensitivity Checks是排除结果对的模型依赖性（Model Dependency）主要是对不同估计方法、带宽选择以及模型中多项式次数的鲁棒性检查，详见上篇推文的第二部分。

当前有三个R语言包可以实现RDD，Thoemmes详细比较了三个包在以上三方面的异同，并归纳如下表。

综合比较而言，rdd的上手成本是最低的也是最基本的选择，如果你熟悉AER，sandwich这几个包的话，那几乎零成本。 rdrobust 在带宽选择和画图方面具备优势，是唯一个输出各类带宽选择及其置信区间的包，当然也是因为这个包CCT开发的。rddtools的强项在于可以自动进行预设与敏感性检查，它可以输出不同带宽选择的估计结果，还可以自动进行安稳剂检验，这是其他包不具备的。如果对于R语言比较陌生的新手，可以选择rddapp这个带有界面的包进行操作。

作者使用一个数据实例展示了三个包联合使用的效果，这里简单归纳一下。

首先展示一下配置变量与结果变量散点图，然后是三步实际操作。

1 第一步：预设检查

（1）断点检验

（2）安慰剂检测

（3）其他前定变量检测

没有画图，文中直接汇报了检验的p值

2 第二步：估计

3 第三步：敏感测试

借着这个图多说一句，带宽的选择实际上是Bias与显著性的权衡取舍，BW越小越无偏，但是也意味着SE会更大，从而导致越发不显著

当然目前还有很多是这几个R包实现不了的，作者列举了两个，一是多个配置变量的估计，二是关于RDD统计解释力（Statistical Power）的估计。

上述步骤实现的R代码如下：

#### PREPERATION ####

library(ggplot2)

library(dplyr)

library(rdd)

library(rddtools)

library(rdrobust)

# set working directory to file location

setwd("C:\\Users\\fjt36-admin\\Dropbox\\WORK\\2016 Winter Cornell\\rdd software jebs")

# read Stata Dataset generated by "./ICPSR_04091/DS0001/04091-0001-Setup.do"

dat = foreign::read.dta(’ICPSR_04091/ICPSR_04091.dta’)

#distribution of treatment, and treatment effect from randomized experiment

xtabs( ~ DC_TRT, dat)

lm_rct1 = lm(SBIQ48 ~ DC_TRT, dat)

summary(lm_rct1)

#treatment effect of randomized experiment interacted with mother’s IQ at time 0

#IQ is centered to provide average treatment effect

lm_rct2 = lm(SBIQ48 ~ DC_TRT * scale(MOMWAIS0), dat)

summary(lm_rct2)

#subset data to create artifical sharp RDD

#here assuming that only mothers with an IQ lower than the median IQ (85) are

#selected for the treament

dat_s = filter(dat,

(MOMWAIS0 >= median(MOMWAIS0) & !DC_TRT) |

(MOMWAIS0 < median(MOMWAIS0) & DC_TRT),!is.na(SBIQ48))

#table to confirm RDD

table(dat_s$DC_TRT, dat_s$MOMWAIS0)

table(dat_s$DC_TRT)

#graphical inspection of all data and then data of only the sharp RDD, here using

#loess smoother with no binning

dat$DC_TRT_LAB <- factor(dat$DC_TRT, labels = c("Control", "Treatment"))

dat_s$DC_TRT_LAB <-

factor(dat_s$DC_TRT, labels = c("Control", "Treatment"))

p1 <- ggplot(dat) +

aes(

x = MOMWAIS0,

y = SBIQ48,

group = DC_TRT_LAB,

color = DC_TRT_LAB,

shape = DC_TRT_LAB

) +

geom_point(aes(shape = DC_TRT_LAB), size = 2) +

geom_smooth(method = ’loess’,

formula = y ~ x,

aes(linetype = DC_TRT_LAB)) + theme_bw() + ylab("Child IQ at age 2") +

xlab("Mother’s IQ\n\n(a)") +

xlim(c(60, 110)) + ylim(c(60, 130)) +

scale_colour_manual(values = c("darkgrey", "black")) +

theme(legend.position = "none")

p2 <- ggplot(dat_s) +

aes(

x = MOMWAIS0,

y = SBIQ48,

group = DC_TRT_LAB,

color = DC_TRT_LAB,

shape = DC_TRT_LAB

) +

geom_point(aes(shape = DC_TRT_LAB), size = 2) +

geom_smooth(method = ’loess’,

formula = y ~ x,

aes(linetype = DC_TRT_LAB)) + theme_bw() + ylab("Child IQ at age 2") +

xlab("Mother’s IQ\n\n(b)") + xlim(c(60, 110)) + ylim(c(60, 130)) +

scale_colour_manual(values = c("darkgrey", "black")) +

theme(legend.title = element_blank()) +

scale_shape_discrete(

name = "DC_TRT",

breaks = c("0", "1"),

labels = c("Treatment", "Control")

)

multiplot(p1, p2, cols = 2)

#naive treatment effect from RDD data

lm_naive = lm(SBIQ48 ~ DC_TRT, dat_s)

summary(lm_naive)

#RDD effect using global, parametric approach "by hand" without any package

lm_rct3 = lm(SBIQ48 ~ DC_TRT * scale(MOMWAIS0), dat_s)

summary(lm_rct3)

#### "rdd" package ####

#McCrary sorting test as implemented in RDD

rdd::DCdensity(dat_s$MOMWAIS0, median(dat$MOMWAIS0), ext.out = TRUE)

# M1, using defaults

lm_rdd = rdd::RDestimate(SBIQ48 ~ MOMWAIS0, dat_s,

cutpoint = median(dat$MOMWAIS0))

summary(lm_rdd)

#plot of RDD as implemented in RDD

plot(lm_rdd)

#single non-outcome covariate test

lm_rdd_nonoutcome = rdd::RDestimate(APGAR5 ~ MOMWAIS0, dat_s,

cutpoint = median(dat$MOMWAIS0))

summary(lm_rdd_nonoutcome)

#### "rdrobust" package ####

# M2

llm_rdrobust1 = rdrobust::rdrobust(dat_s$SBIQ48, dat_s$MOMWAIS0,

c = median(dat$MOMWAIS0))

print(llm_rdrobust1)

#plot of the RDD as implemented in rdrobust, with addition of confidence interval

rdrobust::rdplot(dat_s$SBIQ48,

dat_s$MOMWAIS0,

c = median(dat$MOMWAIS0),

ci = 95)

# M3

bw = rdrobust::rdbwselect_2014(dat_s$SBIQ48,

dat_s$MOMWAIS0,

c = 85,

bwselect = ’IK’)

llm_rdrobust2 = rdrobust::rdrobust(

dat_s$SBIQ48,

dat_s$MOMWAIS0,

c = median(dat$MOMWAIS0),

h = bw$bws[1],

b = bw$bws[2]

)

print(llm_rdrobust2)

# M4

llm_rdrobust3 = rdrobust::rdrobust(dat_s$SBIQ48,

dat_s$MOMWAIS0,

c = median(dat$MOMWAIS0),

h = lm_rdd$bw[1])

print(llm_rdrobust3)

#### "rddtools" package ####

dat_rddtools = rddtools::rdd_data(

y = SBIQ48,

x = MOMWAIS0,

data = dat_s,

cutpoint = median(dat$MOMWAIS0)

)

# M5

llm_rddtools = rddtools::rdd_reg_np(dat_rddtools)

print(llm_rddtools)

# rddtools uses the same density test and plot as rdd

# (call DCDensity from rdd package)

rddtools::dens_test(llm_rddtools)

#default plot of rddtools package

plot(llm_rddtools)

#the default behavior of plotPlacebo() triggers a error, due to NA returned from

# the recalculations of bw by rddtools::rdd_bw_ik()

rddtools::plotPlacebo(llm_rddtools,

same_bw = T,

from = .25,

to = .75)

#sensitivity to bandwidth choice

rddtools::plotSensi(llm_rddtools,

from = 50,

to = 2,

by = -1)

# M6

lm_rddtools1 = rddtools::rdd_reg_lm(dat_rddtools)

print(lm_rddtools1)

# M7

lm_rddtools2 = rddtools::rdd_reg_lm(dat_rddtools,

bw = rddtools::rdd_bw_ik(dat_rddtools))

print(lm_rddtools2)

# M8

lm_rddtools3 = rddtools::rdd_reg_lm(dat_rddtools, bw = lm_rdd$bw[1])

print(lm_rddtools3)

#### "rddapp" package ####

#M9

#this is assuming that the beta version of rddapp is installed - request from authors

lm_rddapp = rddapp::RDestimate(SBIQ48 ~ MOMWAIS0, dat_s,

cutpoint = median(dat$MOMWAIS0))

summary(lm_rddapp)

计量经济圈是中国计量第一大社区，我们致力于推动中国计量理论和实证技能的提升，圈子以海内外高校研究生和教师为主。我们已经建立了计量经济圈社群，受到海内外高校研究生群体的欢迎，大家在里面分享前沿计量资料和热烈讨论各种计量方法。只要能够沉下心来积淀一段时间，我们相信这是对你的研究工作是大有裨益的，欢迎的你加入到咱们这个大家庭（戳这里）。进去之后一定要看“群公告”，不然接收不了群息，也不知道怎么进入咱们的微信群和计量论坛。

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