断点回归RDD：样本少时如何做？
温馨提示： 定期 清理浏览器缓存，可以获得最佳浏览体验。
编译：李琼琼 (山东大学)
Email: lqqflora@163.com
Source: Cattaneo, Matias, D, et al. Inference in regression discontinuity designs under local randomization[J]. The Stata Journal, 2016. -PDF-
1. 背景
断点回归被认为是局部的随机实验，因此随机推断的方法也可以适用于断点回归，但原则上需要使用断点 (cutoff) 附近的观测值。另一方面，经典的断点回归需要较大的样本量，因为无论是参数估计还是非参数估计，均要求分配变量的密度在断点处是连续的，并且做非参数估计时，使用局部多项式做非参数统计推断，也依靠了大样本渐进理论。这就导致了如果带宽过小或者数据集本身就不大，平滑性检验通不过或者估计效果会变得不好的问题。
为了解决样本量有限的问题，美国普林斯顿大学的 Cattano 教授及其合作者于 2016 年开发了一套 rdlocrand Stata 命令包，运用了随机实验的方法来解决有效样本量 (带宽内的样本数量) 过少的问题。
rdlocrand package 包含了四个命令，本文结合 Cattaneo et al. (2015，[PDF]) 关于美国 1994-2010 年参议院选举的样本数据对这四个命令及其功能进行介绍。
准备 ：安装需要使用的命令包
ssc install st0435.pkg,from (http://www.stata-journal.com/software/sj16-2/)
2. 有限样本的断点回归
2.1 RDD 适用性检验和最优带宽的选择
Note: Cattaneo et al. (2015) 使用的参议院选举数据和 David S. Lee (2007，[PDF])使用选举数据在样本的时间段和变量选择上有区别，另外对"在位优势"的定义也不同，David S. Lee 认为 t 期胜出选票的份额对 t+1 期的选举结果有正向影响说明存在"在位优势"，而 Cattaneo et al. 则把 "在位优势" 定义为 t 期对 t+2 期的影响。
rdwinselect 命令有两个作用，(1) 利用前定变量(除配置变量以外的其他解释变量) 做 RDD 适用性检验; (2) 在适用性条件被满足的前提下，挑选最优带宽，适用于有限样本和大样本。基本的语法格式如下：
rdwinselect runvar [covariates] [if] [in], [, options]
其中，主要选项为：
runvar: 分配变量covariates: 协变量cutoff(): 分配变量临界值, 默认为 c(0)obsmin(): 设定左、右带宽包括最少的观测值的个数wmin(): 设定模型带宽长度的最小临界值wstep: 设定模型每次带宽增加的长度nwindows: 设定带宽的个数approximate: 使用大样本渐进理论，默认基于随机推断p(): 结果变量修正模型的阶数，默认 p(0)plot: 画出协变量检验的最小 p 值图reps(): 随机试验重复的次数，默认 reps(1000)
命令的使用
.sysuse rdlocrand_senate, clear //调用数据
.global covariates presdemvoteshlag1 population demvoteshlag1 demvoteshlag2 ///
.demwinprv1 demwinprv2 dopen dmidterm //把分配变量以外的所有解释变量设为全局变量 covariates
.rdwinselect demmv $covariates, cutoff(0) //为RD回归选择带宽， demmv是分配变量，代表民主党超过最大对手党派的投票份额（单位：%）
*--- 使用 rdwinselect 命令 table1
Window selection for RD under local randomization
Cutoff c = 0.00 | Left of c Right of c Number of obs = 1390
------------------+----------------------- Order of poly = 0
Number of obs | 640 750 Kernel type = uniform
1st percentile | 6 7 Reps = 1000
5th percentile | 32 37 Testing method = rdrandinf
10th percentile | 64 75 Balance test = ttest
20th percentile | 128 150
| Bal. test Var. name Bin. test
Window length /2 | p-value (min p-value) p-value Obs<c Obs>=c
------------------+-------------------------------------------------------------
0.529 | 0.210 demvoteshlag2 0.327 10 16
0.733 | 0.262 dopen 0.200 15 24
0.937 | 0.132 dopen 0.126 16 27
1.141 | 0.044 dopen 0.161 20 31
1.346 | 0.229 dmidterm 0.382 28 36
1.550 | 0.102 dmidterm 0.728 35 39
1.754 | 0.075 dmidterm 0.747 41 45
1.958 | 0.046 dmidterm 0.602 43 49
2.163 | 0.075 dmidterm 0.480 45 53
2.367 | 0.132 dopen 0.637 53 59
Variable used in binomial test (running variable): demmv
Covariates used in balance test: presdemvoteshlag1 population demvoteshlag1 demv
> oteshlag2 demwinprv1 demwinprv2 dopen dmidterm
Recommended window is [-.733; .733] with 39 observations (15 below, 24 above).
解读：rdwinselect 命令推荐的最优窗口是 [-.733; .733], 带宽为 0.733 - (-0.733) = 1.467，断点左右两边观测值个数分别为 15 和 24.
*--- 使用 rdwinselect 命令 figure1
.quietly rdwinselect demmv $covariates, wmin(.5) wstep(.125) ///
nwin(80) approximate plot //quietly 表示安静地执行命令不显示结果，只出现图
解读：协变量检验从带宽窗口 [0.5; 0.5] 开始以 0.125 的长度依次增长，可以看出第三个带宽窗口 [0.75, 0.75] 的最小 p 值超过了0.15 线，是符合协变量检验的最优带宽, 这和前面推荐的带宽窗口 [-.733; .733] 非常接近。
2.2 断点回归估计
rdrandinf 命令运用随机推断，对断点回归进行估计，基本的语法格式如下：
rdrandinf outvar runvar [if] [in], [, options]
其中，主要选项为：
runvar: 分配变量covariates: 协变量cutoff(): 断点，分配变量临界值, 默认为 c(0)wl(): 设定带宽窗口的左边界wr(): 设定带宽窗口的右边界covariates: 协变量kernel: 选择核函数, uniform (均匀核,默认), triangular (三角核), epanechnikov(二次核)fuzzy(): 设定模糊断点回归中的内生处理变量
命令的使用
*===: 使用 rdrandinf 命令 table2
*-处理效应估计，设定带宽为 [0.75; 0.75]
*-变量 demvoteshfor2 为结果变量，即 t+2 期的民主党获得投票份额
. rdrandinf demvoteshfor2 demmv, wl(-.75) wr(.75)
Selected window = [-.75 ; .75]
Running randomization-based test...
Cutoff c = 0.00 | Left of c Right of c Number of obs = 1390
------------------+----------------------- Order of poly = 0
Number of obs| 595 702 Kernel type = uniform
Eff. Number of obs| 15 22 Reps = 1000
Mean of outcome| 42.808 52.497 Window = set by user
S.D. of outcome| 7.042 7.742 H0: tau = 0.000
Window| -0.750 0.750 Randomization = fixed margins
Outcome: demvoteshfor2. Running variable: demmv.
--------------------------------------------------------------------------
| Finite sample Large sample
| --------------- -------------------------------
Statistic| T P>|T| P>|T| Power vs d = 3.52
---------------+----------------------------------------------------------
Diff. in means| 9.689 0.000 0.000 0.300
--------------------------------------------------------------------------
解读：处理效应为 9.689, 表明民主党于 t 期控制参议院导致了民主党于 t+2 期获得投票份额在断点提高了约 10%, 在 1% 的显著性水平下显著。 Finite sample 下面的 p 值是基于随机检验的结果，而 large sample 下面的 p 值则基于渐进分布理论。
2.3 敏感性分析 (一)
rdsensitivity 命令分析了对于不同的带宽, 处理效应估计值对应的 p 值和置信区间的敏感性，基本语法如下：
rdsensitivity outvar runvar [if] [in], [, options]
其中，主要选项为：
outvar: 结果变量runvar: 分配变量covariates: 协变量cutoff(): 断点，分配变量临界值, 默认为 c(0)wlist: 设定被估计带宽长度的范围tlist(): 处理效应估计值的范围(整数)saving(filename): 保存敏感性分析的等高线图nodots: 删除重复的点nodraw: 删除等高线verbose: 输出 p 值的矩阵结果reps: 随机试验重复的次数，默认 reps(1000)
命令的使用
*===: 使用 rdsensitivity 命令 table3
. rdsensitivity demvoteshfor2 demmv, ///
wlist(.75(.25)2) tlist(0(1)20) ///
nodots verbose
//此命令运行时间较长，要耐心等待
Running randomization-based test...
Randomization-based test complete.
| .75 1 1.25 1.5 1.75 2
----+------------------------------------------------
0 | 0 .001 0 0 0 0
1 | 0 .001 0 0 0 0
2 | .001 .001 .002 0 0 0
3 | .013 .004 .003 0 0 0
4 | .03 .009 .007 .002 .001 0
5 | .068 .023 .018 .013 .003 .005
6 | .147 .062 .042 .037 .039 .029
7 | .309 .144 .093 .088 .106 .092
8 | .518 .306 .173 .239 .233 .262
9 | .788 .534 .299 .427 .484 .569
10 | .918 .869 .497 .731 .83 .939
11 | .608 .844 .756 .907 .839 .668
12 | .378 .514 .969 .574 .496 .36
13 | .201 .268 .665 .323 .231 .134
14 | .102 .139 .428 .154 .09 .036
15 | .04 .051 .254 .064 .035 .007
16 | .019 .016 .13 .022 .009 .002
17 | .008 .006 .073 .004 .003 0
18 | .003 0 .032 .001 0 0
19 | .001 0 .01 .001 0 0
20 | 0 0 .002 0 0 0
解读:
列名称(.75 等)表示窗宽； 行名称 (0 等) 表示处理效应值，也是估计的原假设； 矩阵的值代表 p 值。
当带宽为 [-.75; 0.75] 时，在 95% 的显著性水平下，处理效应估计值为 5 到 14 之间， 故认为带宽为 [-.75; 0.75] 对应的处理效应值在 5 到 14 之间。同样地，带宽为 [2; 2] 时，处理效应值为 7 到 13 之间。
*===：使用 rdsensitivity 命令 figure2
. rdsensitivity demvoteshfor2 demmv, ///
wlist(.75(.25)10) tlist(0(1)20) nodots ///
saving(figure2)
//saving()可以允许后面对等高线图进行颜色外观等的修改
. use figure2, clear
. twoway contour pvalue t w, ///
ccuts(0(0.05)1) ///
ccolors(gray*0.01 gray*0.05 ///
gray*0.1 gray*0.15 gray*0.2 ///
gray*0.25 gray*0.3 gray*0.35 ///
gray*0.4 gray*0.5 gray*0.6 ///
gray*0.7 gray*0.8 gray*0.9 ///
gray ///
black*0.5 black*0.6 black*0.7 ///
black*0.8 black*0.9 black) ///
xlabel(.75(1.25)10) ///
ylabel(0(2)20, nogrid) ///
graphregion(fcolor(none))
解读: 由图像也可以看出带宽为 [-.75; 0.75]时，p 值超过 5% 对应的处理效应值在 5 到 14 之间。
2.4 敏感性分析 (二)
rdrbounds 命令是基于 Rosenbaum (2002) 的局部随机化理论，计算不同偏离程度下 p 值变化的范围。
rdrbounds demvoteshfor2 demmv, expgamma(1.5 2 3) wlist(.5 .75 1) reps(1000)
rdrbounds outvar runvar [if] [in], [, options]
其中，主要选项为：
runvar: 分配变量cutoff(): 断点，分配变量临界值, 默认为 c(0)expgamma(): 需要被估计的 exp(gamma) 的值, 默认为expgamma(1.5 2 2.5 3)reps: 随机试验重复的次数，默认reps(500)
命令的使用
*--- 使用 rdrbounds 命令 table4
.rdrbounds demvoteshfor2 demmv, expgamma(1.5 2 3) wlist(.5 .75 1) reps(1000)
//此命令运行时间也比较长
Calculating randomization p-values...
w = 0.500 0.750 1.000
--------------------------------------------------------------
Bernoulli p-value| 0.005 0.000 0.000
--------------------------------------------------------------
Running sensitivity analysis...
gamma exp(gamma) w = 0.500 0.750 1.000
--------------------------------------------------------------
0.41 1.50 lower bound| 0.004 0.000 0.000
upper bound| 0.024 0.005 0.002
-----------------------------+--------------------------------
0.69 2.00 lower bound| 0.005 0.000 0.000
upper bound| 0.052 0.025 0.008
-----------------------------+--------------------------------
1.10 3.00 lower bound| 0.005 0.000 0.000
upper bound| 0.194 0.145 0.058
--------------------------------------------------------------
解读: rdrbounds 命令的结果分为两部分，第一个是每个带宽对应 Bernoulli trials 的 p 值, 第二个是在不同的 Γ 和 窗宽下，p 值大小临界值的距离，如果 p 值大小临界值的距离越大，说明随机试验偏差推断越敏感。
第二个结果有两种分析方式，(1) 固定窗宽，看随机试验不同的偏离度如何影响 p 值；(2) 固定随机试验的偏离度，看不同的窗宽如何影响 p 值。对于模型推荐的带宽 [-0.75; 0.75], 在各种偏度下结果都较为稳健。
3. 结语
本文主要对 rdlocrand package 的四个命令 rdwinselect rdrandinf rdsensitivity rdrbounds 进行介绍，并通过美国参议院选举的数据来应用这些命令做 RDD 适用性检验和最优带宽的选择、有限样本的RDD 估计, 以及处理效应值和窗宽敏感性分析。
4. 参考文献
[1] Cattaneo M D , Frandsen B R , Titiunik, Rocío. Randomization Inference in the Regression Discontinuity Design: An Application to Party Advantages in the U.S. Senate[J]. Journal of Causal Inference, 2015, 3(1):1-24pdf.
[2] Cattaneo, Matias, D, et al. Inference in regression discontinuity designs under local randomization[J]. The stata journal, 2016pdf.
[3] Lee D S. Randomized experiments from non-random selection in US House elections[J]. Journal of Econometrics, 2008, 142(2): 675-697pdf.
连享会计量方法专题……，https://gitee.com/arlionn/Course
5. Appendix 本文涉及的 Stata 代码
sysuse rdlocrand_senate, clear //调用数据
global covariates presdemvoteshlag1 population demvoteshlag1 demvoteshlag2 ///
demwinprv1 demwinprv2 dopen dmidterm //把分配变量以外的所有解释变量设为全局变量 covariates
*--- 使用 rdwinselect 命令 table1
rdwinselect demmv $covariates, cutoff(0) //为RD回归选择带宽， demmv是分配变量，代表民主党超过最大对手党派的投票份额（单位：%）
*--- 使用 rdwinselect 命令 figure1
quietly rdwinselect demmv $covariates, wmin(.5) wstep(.125) ///
nwin(80) approximate plot //quietly 表示安静地执行命令不显示结果，只出现图
*--- 使用 rdrandinf 命令 table2
rdrandinf demvoteshfor2 demmv, wl(-.75) wr(.75) //处理效应估计，设定带宽为 [0.75; 0.75]
//变量 demvoteshfor2 为结果变量，即 t+2 期的民主党获得投票份额
*--- 使用 rdsensitivity 命令 table3
rdsensitivity demvoteshfor2 demmv, wlist(.75(.25)2) tlist(0(1)20) nodots verbose //命令运行时间较长
*--- 使用 rdsensitivity 命令 figure2
rdsensitivity demvoteshfor2 demmv, wlist(.75(.25)10) tlist(0(1)20) nodots ///
saving(figure2) //saving()可以允许后面对等高线图进行颜色外观等的修改
use figure2, clear
twoway contour pvalue t w, ccuts(0(0.05)1) ccolors(gray*0.01 gray*0.05 ///
gray*0.1 gray*0.15 gray*0.2 gray*0.25 gray*0.3 gray*0.35 ///
gray*0.4 gray*0.5 gray*0.6 gray*0.7 gray*0.8 gray*0.9 gray ///
black*0.5 black*0.6 black*0.7 black*0.8 black*0.9 black) ///
xlabel(.75(1.25)10) ylabel(0(2)20, nogrid) graphregion(fcolor(none))
*--- 使用 rdrbounds 命令 table4
rdrbounds demvoteshfor2 demmv, expgamma(1.5 2 3) wlist(.5 .75 1) reps(1000)
//此命令运行时间比较长
相关课程
连享会-直播课 上线了！
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
热门资讯
