👇 连享会 · 推文导航 | www.lianxh.cn
🍎 Stata：Stata基础 | Stata绘图 | Stata程序 | Stata新命令 📘 论文：数据处理 | 结果输出 | 论文写作 | 数据分享 💹 计量：回归分析 | 交乘项-调节 | IV-GMM | 时间序列 | 面板数据 | 空间计量 | Probit-Logit | 分位数回归 ⛳ 专题：SFA-DEA | 生存分析 | 爬虫 | 机器学习 | 文本分析 🔃 因果：DID | RDD | 因果推断 | 合成控制法 | PSM-Matching 🔨 工具：工具软件 | Markdown | Python-R-Stata 🎧 课程：最新专题 | 计量专题 | 关于连享会
🍓 课程推荐：2024 机器学习与因果推断专题
主讲老师：司继春 (上海对外经贸大学) ；张宏亮（浙江大学）
课程时间：2024 年 11 月 9-10 日 ；16-17日
课程咨询：王老师 18903405450（微信）
课程特色 · 2024机器学习与因果推断：
懂原理、会应用。本次课程邀请了两位老师合作讲授，目的在于最大限度地实现理论与应用的有机结合。为期四天的课程，分成两个部分：第一部分讲解常用的机器学习算法和适用条件，以及文本分析和大语言模型；第二部分通过精讲 4-6 篇发表于 Top 期刊的论文，帮助大家理解各类机器学习算法的应用场景，以及它们与传统因果推断方法的巧妙结合。 以 Top 期刊论文为范例。目前多数人的困惑是不清楚如何将传统因果推断方法与机器学习结合起来。事实上，即便是 MIT 和 Harvard 的大牛们也都在「摸着石头过河」。为此，通过论文精讲和复现来学习这部分内容或许是目前最有效的方式了。张宏亮老师此前在浙江大学按照这一模式教授了「因果推断和机器学习」课程，效果甚佳：学生们能够逐渐建立起研究设计的理念，并在构造识别策略时适当地嵌入机器学习方法。
作者： 曹昊煜 (兰州大学)
邮箱： caohy19@lzu.edu.cn
编者按：本文部分内容摘译自下文，特此致谢！
Source：Cattaneo M D, Idrobo N, Titiunik R. A practical introduction to regression discontinuity designs: Foundations[M]. Cambridge University Press, 2019. -PDF-
温馨提示： 文中链接在微信中无法生效。请点击底部「阅读原文」。或直接长按/扫描如下二维码，直达原文：
目录
1. 简介
2. 精确断点回归设计
2.1 RD 相关理论介绍
2.2 RD 效应的局部性质
3. RD 绘图
3.1 图形的绘制
3.2 箱体的设置
4. 基于连续性假定的断点回归分析
4.1 局部多项式的点估计
4.2 点估计的实现
4.3 局部多项式的推断
4.4 统计推断的实现
4.5 一些拓展
5. 断点回归的有效性及其证伪
5.1 前定协变量和安慰剂结果变量
5.2 驱动变量的密度
5.3 安慰剂断点
5.4 断点附近观测值的敏感性
5.5 带宽选择的敏感性
6. 参考文献
7. 相关推文
1. 简介
社会科学家主要对因果关系感兴趣。这种关系在随机试验中很容易得到，但在社会科学中，随机试验会受到道德和实践等诸多因素的限制。在缺少随机试验的前提下，我们往往诉诸于严格的准自然实验干预。其中，RD 方法就是一种最为可信的因果推断分析方法。
在 RD 方法中，所有单位都拥有一个 “得分”，并且得分值高于某个断点的个体接受处理。这一设计的最大特点在于接受处理的概率在门槛值处急剧变化。具体来看，RD 的特征如下：
RD 具有三个明显的特征：得分、断点和处理； RD 的主要特征是处理分组在断点处为得分的分段函数，这一条件是可以直接验证的，除此之外，还有很多证伪检验和相关经验来提高可信度； 除了上述特征，一般的解释、估计和推断也需要一些额外的假设。首先，我们需要定义参数及其识别条件；其次，必须施加假设以保证参数可估计，主要的框架有两类，分别是连续性假设和局部随机化假设。
2. 精确断点回归设计
2.1 RD 相关理论介绍
在 RD 中，每个个体拥有一个得分 (或称为驱动变量 running variable，或指标 index)，并且处理取决于得分是否高于断点。
引入一些符号，假设有 个个体，使用 标注，每个个体的得分表示为 ， 是一个已知的断点。当 时，个体接受处理，否则不接受处理。处理状态可以表示为 ，该形式为示性函数。因此，知道个体得分，就知道个体处理状态，即处理的概率是得分的函数。
但是满足处理条件并不意味着实际受到了处理。当处理条件与实际受到处理的情形一致时，我们称之为精确断点回归 (Sharp RD) 设计，否则称之为模糊断点回归 (Fuzzy RD)。图 1 描述了精确断点回归设计的处理规则，即在断点附近，受到处理的概率从零跳跃至 1。
在使用潜在结果框架分析处理效应时，通常假设个体拥有两个潜在结果 和 。在这一框架下，处理效应定义为潜在结果特征或者分布的比较，例如均值、方差或者分位数，但二者中往往只有一个可以被观测。其中，可观测结果变量表示为：
图 2 描述了精确断点回归中的处理效应。在该图中，绘制了给定得分水平时的平均潜在结果，表示为 和 。
从图中可以看出，平均处理效应是给定任意得分时，回归曲线的垂直距离，但由于无法观测到同一得分值处的两条曲线，所以这一距离无法直接估计。但在 点处是个例外，我们可以同时观测到两条曲线。此时，平均处理效应可以定义为：
需要注意的是，我们捕获的信息是 处的处理效应。换句话说， 只能视为局部处理效应。
对得分值相近但分布在断点两侧的个体施加可比性假定，是 RD 设计的基础，这一思想主要由 Hahn 等 (2001) 的连续性假定刻画。Hahn 等 (2001) 认为，如果关于 的回归函数 和 在 处连续，那么在清晰断点中，有：
2.2 RD 效应的局部性质
与其他方法相比，精确断点回归仅仅在驱动变量的一个点上计算了平均差异，也就是局部的因果关系。换言之，其处理效应只包含有限的外部有效性，即精确断点回归设计的结果对远离断点的个体没有良好的代表性。
一般来讲，断点回归估计结果的外部有效性或代表性取决于特殊的应用。举例来说，在图 3 的左图中， 和 之间的垂直距离在 处明显为正，并且平均潜在结果的差异基本上处处为正。但右图在断点处的差异为 0，并且处理效应有正有负。
提升断点回归的外部有效性是一个重要的话题，除了现有的方法外，更重要的是需要更多的假设。可能包括以下几个方面：
断点附近的回归函数； 局部独立假定； 利用设计的特定特征，如不完美的依从性； 多断点回归。
3. RD 绘图
3.1 图形的绘制
首先，可以通过绘制结果变量和驱动变量的散点图，观察断点前后的变化。但由于很难从中发现 “跳跃”，这种方法很少使用。不过，在这里我们仍然使用 Meyersson (2014) 的原始数据进行一个简单的可视化分析。
关于 Meyersson (2014) 文章介绍，详见「伊斯兰政党上位对女性的影响，我们了解多少？」。
* 数据和代码地址：https://gitee.com/arlionn/data/tree/master/data01/RDD-Stata-R
* 数据初步整理
use "https://file.lianxh.cn/data/RDD_Cattaneo_regdata0.dta", clear
* save RDD_Cattaneo_regdata0.dta, replace // 保存到本地
* use RDD_Cattaneo_regdata0.dta, clear
rename hischshr1520f Y
rename i94 T
gen X = vi1*100
replace Y = Y*100
cls
* 对应的Stata代码
twoway (scatter Y X, mcolor(black) msize(vsmall) xline(0, lcolor(black))), ///
graphregion(color(white)) ytitle(Outcome) xtitle(Score)
# 对应的R代码
plot(X, Y, xlab = "Score", ylab = "Outcome", col = 1, pch = 20,
+ cex.axis = 1.5, cex.lab = 1.5) abline(v = 0)
图 4 为散点图的绘制结果。尽管该结果对分析数据的分布有很大帮助，但我们很难从中发现有关 “跳跃” 的信息。因此，基于原始数据的散点图，即使是在大样本条件下，也很难清晰地得到有效结论。
一个更有效的方法是在绘图前加总或平滑数据。典型的 RD 图示应当提供两方面的信息：
一是总体的多项式拟合，用实线表示； 二是使用散点表示简单均值。
前者是一个简单的平滑近似，可以使用四阶或者五阶多项式在断点左右分别拟合结果变量和得分的关系。结果变量的简单均值则需要在得分的不相交区间或者 “箱体” 内计算，然后使用散点的形式绘制在图上。这两种图示均可视为对未知方程的 (非) 平滑近似。最重要的是，在标准的 RD 图示中，全局多项式是使用全样本估计的，而非区间样本。
图 5 绘制的是总体拟合和局部均值图。该图可以使用 rdplot 命令直接实现，在后文中会有介绍。在这张图中，我们可以看出回归函数很可能是非线性的，同时可以在断点处发现一个明显的向上跳跃，表明当伊斯兰政党的领导人当选时，当地的女性受教育水平可能会有所提高。
3.2 箱体的设置
选择箱体的形式
在绘制 RD 图示过程中有两种不同类型的箱体，即等长度箱体和等样本箱体，我们分别称之为 “等间隔 (evently-spaced, ES)” 和 “等数量 (quantile-spaced, QS)”。
为了更细致的定义这些箱体，我们假设驱动变量在 之间取值。我们使用正号和负号下标表示处理组和控制组，分别构造不相交的区间，并使用 ， 表示左右箱体的总数。定义：
上述的定义方式说明 划分了驱动变量的范围 ，其中心为断点处。
在进一步定义区间之前，我们首先约定 和 表示控制组和处理组子区间内的样本个数， 为取整函数。接下来我们可以分别定义两种类型的箱体。
Evenly-spaced（ES）箱体：区分驱动变量取值范围的不相交区间，所有区间长度相同；
Quantile-spaced（QS）箱体：区分驱动变量取值范围的不相交区间，每个区间大致包含相同数量的观测个数。
两种箱体具有一定的差异性。尽管 ES 箱体具有相同的长度，但每个区间内的样本数是不同的。如果样本分布并不均匀，那么区间内均值计算结果的精度可能出现或多或少的差异。而 QS 箱体在每个区间内包含大致相同的样本数，因此所有区间上的均值结果是可比较的。除此之外，QS 箱体还具有一个优势，那就是提供了观测值密度的大致分布。
我们使用 rdplot 命令生成 RD 图示，并比较不同箱体之间的差异。其中，箱体的选择主要使用的选项是 binselect。具体来看：
将断点两侧各划分为 20 个等长的区间，使用的选项是 nbins；使用 QS 箱体绘制，只需要将 binselect中的参数改为qs。
* ES箱体对应的 Stata 代码
rdplot Y X, nbins(20 20) binselect(es) graph_options(graphregion(color(white)) ///
xtitle(Score) ytitle(Outcome))
# ES箱体对应的 R 代码
out = rdplot(Y, X, nbins = c(20, 20), binselect = "es", y.lim = c(0,
+ 25), cex.axis = 1.5, cex.lab = 1.5)
* QS箱体对应的 Stata 代码
rdplot Y X, nbins(20 20) binselect(qs) graph_options(graphregion(color(white)) ///
xtitle(Score) ytitle(Outcome))
# QS箱体对应的 R 代码
out = rdplot(Y, X, nbins = c(20, 20), binselect = "qs", y.lim = c(0,
+ 25), cex.axis = 1.5, cex.lab = 1.5)
选择箱体的数量
当箱体的形式以 QS 或 ES 的方式选定后，还需要确定断点两侧的箱体数量。接下来我们介绍两种在给定箱体类型后，以数据驱动，自动进行的选择 和 方法。
第一个方法通过最小化局部均值的积分均方误差 (Integrated mean-squared error, IMSE) 来选择断点两侧的箱体数量。积分均方误差定义为方差和平方误差的期望和。如果我们选择较大的箱体个数，偏误将会很小，因为区间范围更小并且局部常数拟合也会更好，这一方式的缺陷在于每个箱体内的观测值非常少，导致箱体内的可变性非常大。IMSE 方法的思路正在于平衡误差平方和变异。
选择 IMSE 最优的箱体将得到 “追踪” 潜在回归函数的均值，这意味着 IMSE 最优的箱体更有利于刻画回归函数的形状。然而，IMSE 方法通常在局部均值与全局多项式拟合几乎重合时会绘制出非常平滑的图片，在在断点附近难以捕捉数据的局部变异。
IMSE 最优的 和 分别可以做如下的定义：
这里的 是观测值得总数， 为取顶函数， 和 取决于 ES 或者 QS 箱体类型。实践中，未知常数 和 使用原始的、客观的数据驱动方法估计。
在软件中，可以使用 IMSE 方法绘制 ES 箱体的图像，使用的命令仍然是 rdplot 和选项 binselect(es) ，但是注意命令中没有了 nbins(20 20) 。以下命令是在 IMSE 最优条件下绘制的 ES 箱体。
* ES箱体对应的 Stata 代码
rdplot Y X, binselect(es) graph_options(graphregion(color(white)) ///
xtitle(Score) ytitle(Outcome))
# ES箱体对应的 R 代码
out = rdplot(Y, X, binselect = "es", x.lim = c(-100, 100), y.lim = c(0,
+ 25), cex.axis = 1.5, cex.lab = 1.5)
在 ES 箱体中，每个箱体的长度仍然是相同的，但是 IMES 标准导致了断点两侧不同的箱体个数，因此断点两侧的区间长度也不相同。
最后，我们使用 IMSE 标准考察 QS 箱体的形式，使用 binselect(qs) 选项并忽略箱体数量选项来实现这一过程。
* QS箱体对应的 Stata 代码
rdplot Y X, binselect(qs) graph_options(graphregion(color(white)) ///
xtitle(Score) ytitle(Outcome))
# QS箱体对应的 R 代码
out = rdplot(Y, X, binselect = "qs", x.lim = c(-100, 100), y.lim = c(0,
+ 25), cex.axis = 1.5, cex.lab = 1.5)
在 QS 情形下，箱体的数量要略多于 QS 箱体，并且断点两侧的区间长度有所不同，其原因在于算法需要保证每个箱体内拥有一样多的观测值。
第二种选择 和 的方法称为 Minmicking Variance 方法，该方法中箱体数量的选择依据是使区间均值的总体变化等于原始数据散点图的变化。MV 方法选择的箱体数量定义为：
公式中的 为样本数量，常数的确定取决于 ES 或者 QS 箱体类型，并且其大小与 IMSE 最优的选择不同。一般而言， 并且有 ，在软件中使用如下方式使用 MV 方法。
* ES箱体对应的 Stata 代码
rdplot Y X, binselect(esmv) graph_options(graphregion(color(white)) ///
xtitle(Score) ytitle(Outcome))
# ES箱体对应的 R 代码
out = rdplot(Y, X, binselect = "esmv", cex.axis = 1.5, cex.lab = 1.5)
* QS箱体对应的 Stata 代码
rdplot Y X, binselect(qsmv) graph_options(graphregion(color(white)) ///
xtitle(Score) ytitle(Outcome))
# QS箱体对应的 R 代码
out = rdplot(Y, X, binselect = "qsmv", cex.axis = 1.5, cex.lab = 1.5)
4. 基于连续性假定的断点回归分析
本小节中讨论的方法主要基于连续性假定，或者说：
在这一框架下，将 定义为感兴趣的参数，估计方法是在断点两侧使用多项式方法近似回归函数。实践中，我们往往使用最小二乘法拟合结果变量关于驱动变量的多项式函数。当我们使用所有观测值时，这种多项式拟合是全局的，但是，如果我们仅适用断点附近的观测值拟合多项式时，使用的是局部方法或非参数方法。接下来的讨论中，我们主要使用的是局部方法。
4.1 局部多项式的点估计
局部多项式方法采用线性回归，但仅仅使用断点附近的观测值。具体而言，该方法关注的观测值处于 和 区间内，其中带宽 。在带宽内，通常使用加权方法保证临近断点的观测值相对于其他观测具有更高的权重，这一过程主要通过核函数实现 。局部多项式估计由以下步骤组成：
选择多项式阶数 和核函数 ； 选择带宽 ； 对断点右侧的观测值，结果变量 使用加权最小二乘法拟合。解释变量包括常数项和 ，其中 为选定的多项式阶数，权重为 。截距的估计结果 正是 的估计:
对断点左侧的观测值使用同样的方式，截距的估计结果 正是 的估计：
计算精确断点的点估计：。
选择核函数和多项式阶数
核函数 为每一个观测值指派非负的权重，计算依据是观测值对应的驱动变量值与断点的距离 。一般推荐使用三角核函数进行加权，其函数形式为：
其原因是，在最优均方误差带宽条件下，三角核函数可以得到均方误差意义下的最优点估计。虽然三角核函数具有优良的性质，但有时我们仍然可能使用均匀核函数，其函数形式为 ，带宽以外的权重仍为 0，但区间内的权重全部相同，此时的加权最小二乘法等同于线性回归。
除了核函数外，多项式阶数的选取也会影响估计结果。首先，如果选择 0 阶多项式，也就是对常数拟合，在边界处无法得到优良的理论性质。其次，对于给定带宽，增加多项式的阶数不但可以改变近似的程度，也可以增加处理效应的可变性。实践中，0 阶多项式和线性多项式拟合几乎相同，但后者拥有更小的渐进偏误。最后，高阶多项式可能产生过拟合导致结果不可靠。
带宽选择及其实现
带宽 用于控制局部多项式回归使用的样本情况，带宽选择是断点回归分析的基础，因为其直接影响者估计和推断的性质，并且经验分析对带宽选择非常敏感。
选择一个更小的带宽可以减少模型误设的偏误，但同时有可能提高参数估计的方差，因为带宽越窄，使用的样本越少。反过来，使用更大的带宽可以提高估计的精度，但可能提高模型设定的误差。因此，带宽选择又称为 “偏误-方差权衡”。
由于带宽选择非常重要，所以我们需要使用一些数据驱动的方法来避免数据挖掘。多数带宽选择方法正是试图平衡偏误和方法。实践中使用最为广泛的方式即为最小化局部多项式估计的均方误差 (MSE)。断点回归处理效应渐进均方误差的通常形式为：
进一步定义渐进偏误和估计方差为：
标量 和 分别表示处理效应点估计的偏误和方差，其中使用了样本数量和带宽做调整。尽管省略了很多技术细节，我们展示了偏误和方差的一般形式以澄清局部多项式估计中包含最优均方误差选择的关键权衡。
偏误的一般形式由 决定，同时：
以上的两个极限与未知函数的斜率有关，常数 和 与核函数和多项式阶数有关，该计算方法假设了共同的带宽 ，但该表达式也可以扩展到断点左右的不同带宽。
偏误项取决于回归函数的 阶导数。当我们近似 和 时，该近似必存在一个误差。例如我们使用局部线性模型逼近时，误差的主要部分取决于未知函数的二阶导数。
方差 取决于样本大小和带宽，也可做如下的定义：
该定义与结果变量在断点处的条件变动有关， 为断点处驱动变量的密度函数，常数 和 与和函数和多项式阶数有关。
为得到最优的 MSE 点估计 我们选择最小化 MSE 的带宽：
依据该准则，我们可以得到：
该公式很好地表现了偏误和方差之间的权衡，带宽与 成比例，并且当存在较大的渐进方差时，会出现一个较大的带宽，这在直觉上很容易理解，因为带宽较大时会包含更多的观测值，从而降低估计的方差。
有时在断点两侧使用不同的带宽是非常合理的，因为断点回归的处理效应 是两侧的差，使用不同的带宽可以更好地考虑断点两侧的 MSE 近似。实践中，我们可以定义两侧带宽：
当处理组和控制组的偏误或方差很大时，该方法在实践中非常有用。
最优点估计
给定多项式阶数和核函数，局部多项式的断点回归点估计 依赖于一个选定的带宽 ，该带宽在断点两侧允许不同，并由此可以得到一致的和 MSE 最优的点估计。进一步可以证明三角核函数是 MSE 最优的选择。因此，现代的断点回归设计往往使用数据驱动的带宽选择和三角核函数。
4.2 点估计的实现
现在使用 Meyersson (2014) 的数据演示基于局部多项式的断点回归点估计。我们先设定带宽 ，随后再进行最优带宽的选择。核函数暂时使用均匀核函数，也就是说在带宽内的所有观测值权重相同。
* 对应的 Stata 代码
reg Y X if X < 0 & X >= -20
matrix coef_left = e(b)
local intercept_left = coef_left[1, 2]
reg Y X if X >= 0 & X <= 20
matrix coef_right = e(b)
local intercept_right = coef_right[1, 2]
local difference = `intercept_right' - `intercept_left'
dis "The RD estimator is `difference'"
# 对应的 R 代码
out = lm(Y[X < 0 & X >= -20] ~ X[X < 0 & X >= -20])
left_intercept = out$coefficients[1]
out = lm(Y[X >= 0 & X <= 20] ~ X[X >= 0 & X <= 20])
right_intercept = out$coefficients[1]
difference = right_intercept - left_intercept
print(paste("The RD estimator is", difference, sep = " "))
. reg Y X if X < 0 & X >= -20
Source | SS df MS Number of obs = 610
-------------+---------------------------------- F(1, 608) = 12.65
Model | 1114.58439 1 1114.58439 Prob > F = 0.0004
Residual | 53577.4144 608 88.1207473 R-squared = 0.0204
-------------+---------------------------------- Adj R-squared = 0.0188
Total | 54691.9988 609 89.8062377 Root MSE = 9.3873
------------------------------------------------------------------------------
Y | Coef. Std. Err. t P>|t| [95% Conf. Interval]
-------------+----------------------------------------------------------------
X | -.2398405 .0674381 -3.56 0.000 -.3722804 -.1074007
_cons | 12.73323 .7757093 16.41 0.000 11.20983 14.25662
------------------------------------------------------------------------------
. reg Y X if X >= 0 & X <= 20
Source | SS df MS Number of obs = 282
-------------+---------------------------------- F(1, 280) = 0.93
Model | 78.2903403 1 78.2903403 Prob > F = 0.3361
Residual | 23607.8344 280 84.3136941 R-squared = 0.0033
-------------+---------------------------------- Adj R-squared = -0.0003
Total | 23686.1247 281 84.2922587 Root MSE = 9.1822
------------------------------------------------------------------------------
Y | Coef. Std. Err. t P>|t| [95% Conf. Interval]
-------------+----------------------------------------------------------------
X | -.0989301 .1026652 -0.96 0.336 -.3010237 .1031636
_cons | 15.29072 .9384934 16.29 0.000 13.44332 17.13812
------------------------------------------------------------------------------
The RD estimator is 2.557493537948039
结果显示在带宽为 20 时，断点回归的点估计为 2.56，该结果与原文有出入，但总体差距不大。断点两侧的截距估计分别为 15.29 和 12.73。上述的分步回归无法得到点估计的推断，使用交互项的等价形式即可得到处理效应的显著性水平。
* 对应的 Stata 代码
gen T_X = X * T
reg Y X T T_X if X >= -20 & X <= 20
# 对应的 R 代码
T_X = X * T
out = lm(Y[X >= -20 & X <= 20] ~ X[X >= -20 & X <= 20] + T[X >= -20 & X <= 20] + T_X[X >= -20 & X <= 20])
. reg Y X T T_X if X >= -20 & X <= 20
Source | SS df MS Number of obs = 890
-------------+---------------------------------- F(3, 886) = 4.92
Model | 1280.06282 3 426.687606 Prob > F = 0.0022
Residual | 76884.5593 886 86.777155 R-squared = 0.0164
-------------+---------------------------------- Adj R-squared = 0.0130
Total | 78164.6222 889 87.9242094 Root MSE = 9.3154
------------------------------------------------------------------------------
Y | Coef. Std. Err. t P>|t| [95% Conf. Interval]
-------------+----------------------------------------------------------------
X | -.2398405 .066922 -3.58 0.000 -.3711846 -.1084964
T | 2.816382 1.232233 2.29 0.023 .3979463 5.234818
T_X | .1178932 .1244169 0.95 0.344 -.126293 .3620793
_cons | 12.73323 .7697729 16.54 0.000 11.22244 14.24402
------------------------------------------------------------------------------
在使用交互项的回归方法中，断点回归的点估计为 2.816 ，与之前的结果仍然有一定差异，但几乎处于同一水平。其原因可能是使用的样本存在微小差异。此时的估计使用的是均匀核函数，正如前文所述，在估计中应用三角核函数是更好的选择，其差异在于使用的权重形式不同。
* 对应的 Stata 代码
gen weights = .
replace weights = (1 - abs(X / 20)) if X < 0 & X >= -20
replace weights = (1 - abs(X / 20)) if X >= 0 & X <= 20
reg Y X [aw = weights] if X < 0 & X >= -20
matrix coef_left = e(b)
local intercept_left = coef_left[1, 2]
reg Y X [aw = weights] if X >= 0 & X <= 20
matrix coef_right = e(b)
local intercept_right = coef_right[1, 2]
local difference = `intercept_right' - `intercept_left'
dis "The RD estimator is `difference'"
# 对应的 R 代码
w = NA
w[X < 0 & X >= -20] = 1 - abs(X[X < 0 & X >= -20]/20)
w[X >= 0 & X <= 20] = 1 - abs(X[X >= 0 & X <= 20]/20)
out = lm(Y[X < 0] ~ X[X < 0], weights = w[X < 0])
left_intercept = out$coefficients[1]
out = lm(Y[X >= 0] ~ X[X >= 0], weights = w[X >= 0])
right_intercept = out$coefficients[1]
difference = right_intercept - left_intercept
print(paste("The RD estimator is", difference, sep = " "))
. reg Y X [aw = weights] if X < 0 & X >= -20
(sum of wgt is 304.1711338609457)
Source | SS df MS Number of obs = 610
-------------+---------------------------------- F(1, 608) = 7.74
Model | 653.698673 1 653.698673 Prob > F = 0.0056
Residual | 51378.5978 608 84.5042726 R-squared = 0.0126
-------------+---------------------------------- Adj R-squared = 0.0109
Total | 52032.2964 609 85.4389104 Root MSE = 9.1926
------------------------------------------------------------------------------
Y | Coef. Std. Err. t P>|t| [95% Conf. Interval]
-------------+----------------------------------------------------------------
X | -.2224576 .079983 -2.78 0.006 -.3795341 -.0653811
_cons | 12.85216 .6618279 19.42 0.000 11.55242 14.15191
------------------------------------------------------------------------------
. reg Y X [aw = weights] if X >= 0 & X <= 20
(sum of wgt is 177.2444166690111)
Source | SS df MS Number of obs = 282
-------------+---------------------------------- F(1, 280) = 0.51
Model | 45.0201751 1 45.0201751 Prob > F = 0.4737
Residual | 24492.1672 280 87.4720257 R-squared = 0.0018
-------------+---------------------------------- Adj R-squared = -0.0017
Total | 24537.1874 281 87.3209515 Root MSE = 9.3526
------------------------------------------------------------------------------
Y | Coef. Std. Err. t P>|t| [95% Conf. Interval]
-------------+----------------------------------------------------------------
X | -.0957842 .1335135 -0.72 0.474 -.3586019 .1670334
_cons | 15.27445 .8872344 17.22 0.000 13.52795 17.02095
------------------------------------------------------------------------------
The RD estimator is 2.422284803520951
在使用三角核函数的方法中，断点回归的点估计量变为 2.42，这也说明估计结果相对依赖于核函数的选择。尽管使用 reg 命令可以很好地澄清断点回归的本质，但该方法计算的置信区间很可能是无效的。因此我们可以使用 rdrobust 命令，该命令可以控制带宽选择、多项式阶数等特征。
* 对应的 Stata 代码
rdrobust Y X, kernel(uniform) p(1) h(20)
# 对应的 R 代码
out = rdrobust(Y, X, kernel = "uniform", p = 1, h = 20)
. rdrobust Y X, kernel(uniform) p(1) h(20)
Sharp RD estimates using local polynomial regression.
Cutoff c = 0 | Left of c Right of c Number of obs = 2634
-------------------+---------------------- BW type = Manual
Number of obs | 2317 317 Kernel = Uniform
Eff. Number of obs | 610 282 VCE method = NN
Order est. (p) | 1 1
Order bias (q) | 2 2
BW est. (h) | 20.000 20.000
BW bias (b) | 20.000 20.000
rho (h/b) | 1.000 1.000
Outcome: Y. Running variable: X.
--------------------------------------------------------------------------------
Method | Coef. Std. Err. z P>|z| [95% Conf. Interval]
-------------------+------------------------------------------------------------
Conventional | 2.5575 1.2253 2.0873 0.037 .156007 4.95898
Robust | - - 1.2480 0.212 -1.26048 5.67978
--------------------------------------------------------------------------------
在初步使用 rdrobust 命令时，可以使用 kernel 设定核函数形式， p 设定多项式阶数，h 设定带宽。结果中包含了很多细节，右上角显示，回归中使用的样本个数为 2634 ，核函数为均匀核函数，手动选择带宽大小，未使用稳健标准误。左侧的结果类似于 rdplot 的结果，分别报告了断点左右的样本个数、估计中使用的样本数量、带宽大小和多项式阶数。最后报告的是点估计结果，从中可知 与分组回归的结果非常近似。
使用 rdrobust 可以非常轻易地使用三角核函数，只需要改变选项中的参数即可。
* 对应的 Stata 代码
rdrobust Y X, kernel(Triangular) p(1) h(20)
# 对应的 R 代码
out = rdrobust(Y, X, kernel = "Triangular", p = 1, h = 20)
. rdrobust Y X, kernel(triangular) p(1) h(20)
Sharp RD estimates using local polynomial regression.
Cutoff c = 0 | Left of c Right of c Number of obs = 2634
-------------------+---------------------- BW type = Manual
Number of obs | 2317 317 Kernel = Triangular
Eff. Number of obs | 610 282 VCE method = NN
Order est. (p) | 1 1
Order bias (q) | 2 2
BW est. (h) | 20.000 20.000
BW bias (b) | 20.000 20.000
rho (h/b) | 1.000 1.000
Outcome: Y. Running variable: X.
--------------------------------------------------------------------------------
Method | Coef. Std. Err. z P>|z| [95% Conf. Interval]
-------------------+------------------------------------------------------------
Conventional | 2.4223 1.3315 1.8192 0.069 -.187419 5.03199
Robust | - - 0.9210 0.357 -1.95833 5.43033
--------------------------------------------------------------------------------
断点回归的点估计结果为 0.2422 与加权最小二乘法的估计结果基本一致。最后，为了减少估计中的近似误差，使用二阶多项式拟合，同样仅仅需要设定 p(2) 即可。
* 对应的 Stata 代码
rdrobust Y X, kernel(Triangular) p(2) h(20)
# 对应的 R 代码
out = rdrobust(Y, X, kernel = "Triangular", p = 2, h = 20)
. rdrobust Y X, kernel(triangular) p(2) h(20)
Sharp RD estimates using local polynomial regression.
Cutoff c = 0 | Left of c Right of c Number of obs = 2634
-------------------+---------------------- BW type = Manual
Number of obs | 2317 317 Kernel = Triangular
Eff. Number of obs | 610 282 VCE method = NN
Order est. (p) | 2 2
Order bias (q) | 3 3
BW est. (h) | 20.000 20.000
BW bias (b) | 20.000 20.000
rho (h/b) | 1.000 1.000
Outcome: Y. Running variable: X.
--------------------------------------------------------------------------------
Method | Coef. Std. Err. z P>|z| [95% Conf. Interval]
-------------------+------------------------------------------------------------
Conventional | 1.736 1.8849 0.9210 0.357 -1.95833 5.43033
Robust | - - -0.0978 0.922 -5.12828 4.64067
--------------------------------------------------------------------------------
需要注意的是，此时的点估计结果出现了大幅度下降，并且不再显著，这一现象在断点回归中经常出现。除非高阶项的近似确实全部为 0 ，否则合并这些项会减少近似误差但也会改变估计结果。
以上的分析均是在假定 的前提下进行的，但我们并不知道带宽等于 20 在误差和方差中意味着什么，或者说它在推断中是否合理。为此，我们使用 rdwselect 命令选择最优带宽。
* 对应的 Stata 代码
rdbwselect Y X, kernel(triangular) p(1) bwselect(mserd)
# 对应的 R 代码
out = rdbwselect(Y, X, kernel = "triangular", p = 1, bwselect = "mserd")
. rdbwselect Y X, kernel(triangular) p(1) bwselect(mserd)
Bandwidth estimators for sharp RD local polynomial regression.
Cutoff c = | Left of c Right of c Number of obs = 2634
-------------------+---------------------- Kernel = Triangular
Number of obs | 2317 317 VCE method = NN
Min of X | -100.000 0.000
Max of X | -0.098 99.051
Order est. (p) | 1 1
Order bias (q) | 2 2
Outcome: Y. Running variable: X.
--------------------------------------------------------------------------------
| BW est. (h) | BW bias (b)
Method | Left of c Right of c | Left of c Right of c
-------------------+------------------------------+-----------------------------
mserd | 18.034 18.034 | 30.330 30.330
--------------------------------------------------------------------------------
MSE 最优的带宽选择取决于多项式阶数和核函数，因此在 rdwselect 中必须设定二者的形式。结果中说明了带宽的形式，本例中是 mserd，即断点两侧同样长度的带宽。结果下方是带宽的选择结果，带宽 是用于点估计的带宽我们有时称之为 “主带宽”，另一个带宽 将用于推断中需要的误差项估计。结果中用于局部线性回归的带宽为 18.034，在断点两侧相同，有时我们也需要在断点两侧确定不同长度的带宽。
* 对应的 Stata 代码
rdbwselect Y X, kernel(triangular) p(1) bwselect(msetwo)
# 对应的 R 代码
out = rdbwselect(Y, X, kernel = "triangular", p = 1, bwselect = "msetwo")
. rdbwselect Y X, kernel(triangular) p(1) bwselect(msetwo)
Bandwidth estimators for sharp RD local polynomial regression.
Cutoff c = | Left of c Right of c Number of obs = 2634
-------------------+---------------------- Kernel = Triangular
Number of obs | 2317 317 VCE method = NN
Min of X | -100.000 0.000
Max of X | -0.098 99.051
Order est. (p) | 1 1
Order bias (q) | 2 2
Outcome: Y. Running variable: X.
--------------------------------------------------------------------------------
| BW est. (h) | BW bias (b)
Method | Left of c Right of c | Left of c Right of c
-------------------+------------------------------+-----------------------------
msetwo | 19.371 15.830 | 31.017 26.432
--------------------------------------------------------------------------------
选项 msetwo 可以允许断点两侧出现不同的带宽。此时的结果中断点左侧的带宽大于右侧的带宽。一旦我们完成带宽选择时，我可们以将之以参数形式传递到 rdrobust 中，但实际上，rdrobust 也包含了带宽选择的选项。
* 对应的 Stata 代码
rdrobust Y X, kernel(triangular) p(1) bwselect(mserd)
# 对应的 R 代码
out = rdrobust(Y, X, kernel = "triangular", p = 1, bwselect = "mserd")
. rdrobust Y X, kernel(triangular) p(1) bwselect(mserd)
Sharp RD estimates using local polynomial regression.
Cutoff c = 0 | Left of c Right of c Number of obs = 2634
-------------------+---------------------- BW type = mserd
Number of obs | 2317 317 Kernel = Triangular
Eff. Number of obs | 551 272 VCE method = NN
Order est. (p) | 1 1
Order bias (q) | 2 2
BW est. (h) | 18.034 18.034
BW bias (b) | 30.330 30.330
rho (h/b) | 0.595 0.595
Outcome: Y. Running variable: X.
--------------------------------------------------------------------------------
Method | Coef. Std. Err. z P>|z| [95% Conf. Interval]
-------------------+------------------------------------------------------------
Conventional | 2.451 1.3868 1.7674 0.077 -.26702 5.16894
Robust | - - 1.4162 0.157 -.881458 5.47327
使用 ereturn list 可以分别得到处理组和控制组的截距估计，二者之差即为断点回归点估计结果。在本例中，这一关系可以表示为 e(tau_cl_r) - e(tau_cl_l) = 15.279 - 12.828 = 2.451，也就是说相对于控制组，伊斯兰领导人的当选大约可以提高 的女性受教育水平。
* 对应的 Stata 代码
rdrobust Y X, p(1) kernel(triangular) bwselect(mserd)
local bandwidth = e(h_l)
rdplot Y X if abs(X) <= `bandwidth', p(1) h(`bandwidth') kernel(triangular)
# 对应的 R 代码
bandwidth = rdrobust(Y, X, kernel = "triangular", p = 1, bwselect = "mserd")$bws[1,1]
out = rdplot(Y[abs(X) <= bandwidth], X[abs(X) <= bandwidth], p = 1, kernel = "triangular", cex.axis = 1.5, cex.lab = 1.5)
4.3 局部多项式的推断
为了更好地考察点估计的性质，我们需要进行假设检验并构造置信区间。直觉上我们应该使用最小二乘法的置信区间，但这时候需要假定模型设定的正确性。所以更加理想的方式是选择一个与 “偏误-方差” 权衡相联系的带宽。
出于这一考虑，有效地推断需要考虑模型误设的影响。实践中，我们选择的 MSE 最优带宽可以用于点估计，但如果将之用于推断则会产生一些问题。其原因在于 MSE 最优带宽没有小到可以消除用于统计推断的标准分布中的主要偏误项。出现这一问题的根源在于这种带宽的选择目的是点估计，因此，使用 OLS 方法在区间 内进行推断通常是无效的。
在断点回归的推断中，通常使用两种方法：一是将带宽 同时用于估计和推断，但是在推断中校正 T 统计量；二是重新选择一个不同的带宽用于推断。
使用 带宽进行估计和推断
我们首先讨论在带宽为 的情形下如何构造有效的推断，局部多项式回归点估计的大样本近似分布为：
其中的 和 为渐进偏误和方差。这一构造方式再次强调了局部估计中带宽选择和模型误设的偏误，方差项 可以使用 (加权) 最小二乘法计算，也可以在其中处理异方差和相关问题。给定局部多项式估计量的渐进分布后，95% 的置信区间可以给定为：
该区间包含了未知的偏误项，其原因在于局部多项式估计是非参数近似，而非假定回归函数就是 阶多项式。如果我们将偏误项忽略，将导致错误的推断，除非该偏误非常之小，或者说局部线性回归与真实的模型非常接近。
接下来，介绍一些用于推断的不同策略，并解释为什么其中的部分方式是无效的。
传统推断方法
第一种方法是传统推断方法，即在使用 MSE 带宽的情况下忽略模型设定不正确的偏误。基于该方法的推断不但是无效的，而且从方法论的角度出发也不合理。
这种过于简单的方式将局部多项式回归视为了一种参数方法，因此忽略了偏误项。在这种情况下，无效的推断几乎必然出现，除非近似偏误非常小以至于可以忽略。此时的 95% 置信区间可以定义为：
这一区间与参数估计的置信区间相同，我们称之为传统的置信区间，记为 ，使用该区间等同于假定估计中选择的多项式与真实的函数 和 非常近似。但由于真实的函数形式未知，所以这一假定很难验证，如果在实践中使用该方法，可能会过度拒绝原假设。
第二种方法是标准偏误校正，该方法的思路是手动估计 MSE 最优带宽下的偏误，然后构造包含偏误估计的置信区间：
标准偏误方法可以允许更大的带宽，但该方法通常很少使用，因为在偏差估计步骤中引入的可变性并没有包含在所使用的方差项中，与 一样， 忽略了方差的可变性。
第三种方法是稳健偏误校正，该方法提供了更为高级的推断过程。即使我们使用 MSE 最优带宽 (非欠平滑)，或者 (违反标准偏误校正的假定)，稳健的偏误校正推断仍然是有效的。稳健的偏误校正方法包含一个新的渐进方差 ，此时的置信区间为：
这一置信区间由断点回归估计量和偏误估计之差以及新的渐进方差构造，在多数情况下，基于该置信区间的推断都是有效的。
| Centered at | Standard Eror | |
|---|---|---|
| Conventional: | ||
| Bias Corrected: | ||
| Robust Bias Corrected: |
上表总结了不同推断方式的差异。传统的 OLS 置信区间 忽略了偏误项，使用的仅仅是传统标准误 。偏误校正置信区间 从点估计中减去了偏误的估计，但忽略了偏误校正过程中引入的变异，也就是说仍然使用的是 。最后，稳健偏误校正置信区间 同时考虑了偏误和标准误的变动，采用了一个同样条件下较大的标准误 。从实践角度出发，稳健偏误校正的重要特征是可以使用 作为估计带宽，而此时传统置信区间是无效的。
使用与估计不同的带宽推断
理论上，传统置信区间的无效性源于在推断中使用了出于点估计目标确定的带宽。虽然 在估计中可以得到一致并且最小化均方误差的估计量，但是使用其构造置信区间则存在极大的挑战。因为 MSE 最优的带宽选择并非总能得到良好的分布近似，有时甚至难以得到有效的分布。
因此，一个可以选择的方式是分离点估计和推断的目标，为每一项工作分别赋予不同的带宽。那么此时我们分别使用 估计处理效应，使用 进行统计推断，前者最小化了均方误差 MSE，后者最小化了覆盖误差 CER，构造的置信区间仍然是稳健偏误校正的置信区间。
4.4 统计推断的实现
本小节仍然使用一阶多项式和三角核函数，演示断点回归中的统计推断。
* 对应的 Stata 代码
rdrobust Y X, kernel(triangular) p(1) bwselect(mserd)
# 对应的 R 代码
out = rdrobust(Y, X, kernel = "triangular", p = 1, bwselect = "mserd")
. rdrobust Y X, kernel(triangular) p(1) bwselect(mserd)
Sharp RD estimates using local polynomial regression.
Cutoff c = 0 | Left of c Right of c Number of obs = 2634
-------------------+---------------------- BW type = mserd
Number of obs | 2317 317 Kernel = Triangular
Eff. Number of obs | 551 272 VCE method = NN
Order est. (p) | 1 1
Order bias (q) | 2 2
BW est. (h) | 18.034 18.034
BW bias (b) | 30.330 30.330
rho (h/b) | 0.595 0.595
Outcome: Y. Running variable: X.
--------------------------------------------------------------------------------
Method | Coef. Std. Err. z P>|z| [95% Conf. Interval]
-------------------+------------------------------------------------------------
Conventional | 2.451 1.3868 1.7674 0.077 -.26702 5.16894
Robust | - - 1.4162 0.157 -.881458 5.47327
--------------------------------------------------------------------------------
在估计结果中，断点回归的处理效应为 2.451，带宽为 18.034，最后两行结果报告了两种不同的统计推断方法。传统的置信区间 可以使用定义计算：
下面一行报告了稳健偏误校正的置信区间，偏误的估计和调整后的标准误均没有报告在结果中，但我们可以发现稳健偏误调整的置信区间长度一般大于传统的置信区间，虽然这一结果在估计和推断使用不同的带宽时可能不成立。
在结果中不报告偏误校正的点估计结果有其原因：当我们使用 MSE 带宽进行点估计时，偏误校正的点估计是次优的，因此在实践中，我们往往报告的是 MSE 最优的点估计，随后在讨论中使用稳健偏误校正的置信区间。但如果我们想要考察不同情形下的估计结果，可以附加 all 选项。
* 对应的 Stata 代码
rdrobust Y X, kernel(triangular) p(1) bwselect(mserd) all
# 对应的 R 代码
out = rdrobust(Y, X, kernel = "triangular", p = 1, bwselect = "mserd", all = TRUE)
. rdrobust Y X, kernel(triangular) p(1) bwselect(mserd) all
Sharp RD estimates using local polynomial regression.
Cutoff c = 0 | Left of c Right of c Number of obs = 2634
-------------------+---------------------- BW type = mserd
Number of obs | 2317 317 Kernel = Triangular
Eff. Number of obs | 551 272 VCE method = NN
Order est. (p) | 1 1
Order bias (q) | 2 2
BW est. (h) | 18.034 18.034
BW bias (b) | 30.330 30.330
rho (h/b) | 0.595 0.595
Outcome: Y. Running variable: X.
--------------------------------------------------------------------------------
Method | Coef. Std. Err. z P>|z| [95% Conf. Interval]
-------------------+------------------------------------------------------------
Conventional | 2.451 1.3868 1.7674 0.077 -.26702 5.16894
Bias-corrected | 2.2959 1.3868 1.6556 0.098 -.422077 5.01389
Robust | 2.2959 1.6211 1.4162 0.157 -.881458 5.47327
--------------------------------------------------------------------------------
结果中出现了三行内容，前两行是传统的置信区间和偏误调整的置信区间，二者使用的标准误相同，而第三行为稳健偏误调整的置信区间，其标准误略大于前两种推断方法。由于仅在传统传统情形中没有减去误差，所以可以推得 。
最后，我们探讨使用 带宽的稳健偏误校正推断，此时需要在命令中改变带宽选择方式。
* 对应的 Stata 代码
rdrobust Y X, kernel(triangular) p(1) bwselect(cerrd)
# 对应的 R 代码
out = rdrobust(Y, X, kernel = "triangular", p = 1, bwselect = "cerrd")
. rdrobust Y X, kernel(triangular) p(1) bwselect(cerrd)
Sharp RD estimates using local polynomial regression.
Cutoff c = 0 | Left of c Right of c Number of obs = 2634
-------------------+---------------------- BW type = cerrd
Number of obs | 2317 317 Kernel = Triangular
Eff. Number of obs | 375 225 VCE method = NN
Order est. (p) | 1 1
Order bias (q) | 2 2
BW est. (h) | 12.164 12.164
BW bias (b) | 30.330 30.330
rho (h/b) | 0.401 0.401
Outcome: Y. Running variable: X.
--------------------------------------------------------------------------------
Method | Coef. Std. Err. z P>|z| [95% Conf. Interval]
-------------------+------------------------------------------------------------
Conventional | 1.7897 1.6255 1.1011 0.271 -1.3961 4.97556
Robust | - - 0.9758 0.329 -1.72607 5.14899
--------------------------------------------------------------------------------
此时的带宽选择为 12.164，小于 MSE 的最优带宽，因此点估计结果更小，但 p 值却更大，估计结果可以被视为一个偏误更小但变异更大的估计结果。在实践中，选择 MSE 带宽或 CER 带宽，以及在断点两侧选择相同的带宽或者不同的带宽，都可能会产生很大的差异，我们可以使用 rdbwselect 考察。
* 对应的 Stata 代码
rdbwselect Y X, kernel(triangular) p(1) all
# 对应的 R 代码
out = rdbwselect(Y, X, kernel = "triangular", p = 1, all = TRUE)
. rdbwselect Y X, kernel(triangular) p(1) all
Bandwidth estimators for sharp RD local polynomial regression.
Cutoff c = | Left of c Right of c Number of obs = 2634
-------------------+---------------------- Kernel = Triangular
Number of obs | 2317 317 VCE method = NN
Min of X | -100.000 0.000
Max of X | -0.098 99.051
Order est. (p) | 1 1
Order bias (q) | 2 2
Outcome: Y. Running variable: X.
--------------------------------------------------------------------------------
| BW est. (h) | BW bias (b)
Method | Left of c Right of c | Left of c Right of c
-------------------+------------------------------+-----------------------------
mserd | 18.034 18.034 | 30.330 30.330
msetwo | 19.371 15.830 | 31.017 26.432
msesum | 16.163 16.163 | 26.628 26.628
msecomb1 | 16.163 16.163 | 26.628 26.628
msecomb2 | 18.034 16.163 | 30.330 26.628
-------------------+------------------------------+-----------------------------
cerrd | 12.164 12.164 | 30.330 30.330
certwo | 13.065 10.677 | 31.017 26.432
cersum | 10.902 10.902 | 26.628 26.628
cercomb1 | 10.902 10.902 | 26.628 26.628
cercomb2 | 12.164 10.902 | 30.330 26.628
--------------------------------------------------------------------------------
前五个结果均为 MSE 最优的带宽，mserd 表示在断点两侧使用长度相同的带宽，mestwo 表示在断点两侧使用长度不同的带宽，msesum 表示最小化回归系数之和的 MSE 的带宽，msecomb1 是最小化 mserd 和 msesum 的带宽，msecomb2 取 msetwo、mserd 和 msesum 的中位数。其余结果的解释与上述内容相同，只不过选取准则从 MSE 变成 CER。
4.5 一些拓展
到目前为止，我们的解释变量仅仅包含驱动变量，也就是假定了所有样本观测值是独立的，本小节我们讨论增加控制变量和聚类问题。关于具体代码实现过程，详见本文「数据文件」。
在分析中加入协变量
最简单的断点回归分析是使用结果变量对驱动变量拟合，但有时需要在其中加入协变量。除非我们需要使用参数假设或者重新定义感兴趣的参数，否则断点两侧的协变量必须平衡，否则协变量的调整不能改进标准断点回归的设计，还可能影响连续性假定。
聚类稳健标准误
另一个问题是组内的聚类问题。当观测值处于某一组内并且怀疑误差项相关时，使用聚类稳健标准误可能更为合理。使用最小二乘法估计和推断，同样可以在局部多项式回归中进行聚类调整，此时最优带宽也会发生变化，因此，使用聚类方差不但会改变推断结果，也会改变点估计的结果。
5. 断点回归的有效性及其证伪
我们在此讨论的有效性检验主要涉及五个方面：
断点回归对于协变量和安慰剂结果变量的处理效应为零； 断点附近得分密度的连续性； 人工断点处的处理效应； 排除断点附近的观测样本； 带宽选择的敏感性。
5.1 前定协变量和安慰剂结果变量
一个最为关键的检验即为断点附近，处理组和控制组的相似性。其基本思想为，如果个体缺乏精确操纵得分的能力，那么获得相同得分的个体将不存在系统性的差异。这种样本特征可以分为两类：一类是断点设计前就已经确定的变量，称之为前定协变量，第二类是断点设计后决定的变量，但由于大量信息的存在，并未受到干预的影响，我们称之为安慰剂结果变量。
重要的是，前定协变量总是可以明确定义，但安慰剂结果变量在应用中有所不同。在断点设计前就已经确定的变量均为前定协变量，相反，某一变量是否为安慰剂结果变量取决于实际的情况。举例来说，如果要研究清洁水与儿童死亡率的关系，那么由车祸导致的儿童死亡率即可作为安慰剂结果变量。
实施前定协变量和安慰剂结果变量最为重要的原则是适用与点估计相同的分析方法，也就是说，对于每一个前定协变量和安慰剂结果变量，我们必须首先选择最优带宽然后使用局部多项式回归估计 “处理效应”，并使用有效推断讨论显著性问题。此时，由于前定协变量或安慰剂结果变量不受干预的影响，因此不应拒绝原假设，否则支持断点回归有效性的连续性假定可能难以满足。
当我们在基于连续性的方法中进行分析时，这种证伪检验使用第四部分讨论的局部多项式技术，检验是否存在针对前定协变量的处理效应。首先仍然使用循环语句和 rdplot 命令图示断点处的情况。
* 对应的 Stata 命令
foreach v in lpop1994 partycount vshr_islam1994 buyuk merkezp merkezi{
rdplot `v' X
}
图示的分析结果并未发现断点处的跳跃，但是统计分析仍然是必要的。为了进行该检验，对每个协变量都要选择最优带宽。从图中可以看出每个协变量使用的带宽与之前都不相同，也具有不同的函数形式，因此对于每一个协变量的估计和检验结果都不相同。这意味着我们使用 rdrobust 进行正式的统计检验时，需要对每一个协变量分别执行同样的证伪检验。
* 对应的 Stata 代码
foreach v in lpop1994 partycount vshr_islam1994 buyuk merkezp merkezi{
rdrobust `v' X
}
在分析这一过程的结果时可以关注两个方面的信息，一是点估计的结果是否都接近于 0，并且是否均不显著的。在令人信服的断点回归分析中，每个前定协变量都需要进行证伪检验，在此提供一套代码将证伪检验的结果以矩阵形式表现出来：
* 对应的 Stata 代码
mat A = J(6,6,.)
local j = 1
local i = 1
foreach v in lpop1994 partycount vshr_islam1994 buyuk merkezp merkezi{
qui rdrobust `v' X
mat A[`i',`j'] = e(h_l)
mat A[`i',`j'+1] = e(tau_cl)
mat A[`i',`j'+2] = e(pv_cl)
mat A[`i',`j'+3] = e(ci_l_cl)
mat A[`i',`j'+4] = e(ci_r_cl)
mat A[`i',`j'+5] = e(N_h_l) + e(N_h_r)
local j = 1
local i = `i' + 1
}
mat rownames A = lpop1994 partycount vshr_islam1994 buyuk merkezp merkezi
mat colnames A = Bw Est p-value CI_l CI_r Obs
matlist A
在本例的结果中，所有协变量的最优带宽均不相同，点估计结果非常接近于 0，从 值来看，所有协变量在断点处均不存在明显的跳跃。需要注意的是，由于此时我们并不关注点估计的结果，因此使用 CER 带宽可能更为合理，此时只需要在上述代码中附加 bwselect(cerrd) 选项即可。除此之外，还建议使用 rdplot 命令绘制断点两侧局部线性回归的结果。
5.2 驱动变量的密度
第二类证伪检验主要用于验证略低于断点的观测值个数是否明显地区别于略高于断点的观测值。其潜在假定是，如果个体没有精确操纵得分的能力，那么略高于断点的观测值个数应当与略低于断点的观测值非常接近，换句话说，当驱动变量没有受到精确的操纵时，驱动变量的密度函数应当是连续的。对于这一假定，我们通常使用驱动变量的直方图对比断点左右的分布情况。
除了图形的方式检验以外，我们还可以进一步使用统计检验验证一般假设，通常称之为密度检验。一类可能的策略即为在断点附近选择一个很小的区间，并在该区间内使用简单的贝努利检验说明 “成功” 的概率是 0.5。该策略检验的是当个体以 50% 的概率被分配在处理组中时，选定区间内受到处理的观测样本的数量与之是否可比。
* 对应的 Stata 代码
. bitest 100 53, 1/2
# 对应的 R 代码
binom.test(53, 100, 1/2)
. bitesti 100 53 1/2
N Observed k Expected k Assumed p Observed p
------------------------------------------------------------
100 53 50 0.50000 0.53000
Pr(k >= 53) = 0.308650 (one-sided test)
Pr(k <= 53) = 0.757941 (one-sided test)
Pr(k <= 47 or k >= 53) = 0.617299 (two-sided test)
结果中的 p 值为 0.53，这意味着这一简单的检验在断点附近没有发现处理组和控制组的操纵行为。然而，在基于连续性的方法中，如何选择合适的区间来进行贝努利检验往往没有明确的准则，因此最好的方式是选择多种区间进行这一假定。
另一种在连续性框架下可以使用的方式是在断点两侧的某一区间内估计观测值的密度函数，估计方法为局部多项式的密度估计，由于该方法不必进行预先的箱体选择，因此更加可信。该方法的原假设为断点两侧不存在对密度函数的操纵，所以不拒绝原假设意味着断点回归设计是有效的。
* 对应的 Stata 代码
rddensity X
# 对应的 R 代码
out = rddensity(X)
. rddensity X ,plot
Computing data-driven bandwidth selectors.
Point estimates and standard errors have been adjusted for repeated observations.
(Use option nomasspoints to suppress this adjustment.)
RD Manipulation test using local polynomial density estimation.
c = 0.000 | Left of c Right of c Number of obs = 2664
-------------------+---------------------- Model = unrestricted
Number of obs | 2334 330 BW method = comb
Eff. Number of obs | 742 318 Kernel = triangular
Order est. (p) | 2 2 VCE method = jackknife
Order bias (q) | 3 3
BW est. (h) | 23.975 30.652
Running variable: X.
------------------------------------------
Method | T P>|T|
-------------------+----------------------
Robust | -0.8958 0.3703
------------------------------------------
P-values of binomial tests. (H0: prob = .5)
-----------------------------------------------------
Window Length / 2 | <c >=c | P>|T|
-------------------+----------------------+----------
0.378 | 9 11 | 0.8238
0.756 | 18 27 | 0.2327
1.135 | 27 34 | 0.4426
1.513 | 40 46 | 0.5900
1.891 | 47 53 | 0.6173
2.269 | 61 61 | 1.0000
2.647 | 74 68 | 0.6749
3.025 | 80 77 | 0.8732
3.404 | 99 84 | 0.3007
3.782 | 118 90 | 0.0609
-----------------------------------------------------
检验结果的 T 值为 -0.8958，相应的 P 值为 0.1633。这意味着在连续性方法中，我们无法拒绝密度函数连续的原假设。
5.3 安慰剂断点
另一种有用的证伪检验是考察在 “虚假的” 或安慰剂断点处是否仍有显著的处理效应。使用这一方法的动机在于，断点回归的关键识别假定是当不存在处理的情况下，处理组和控制组的回归函数是连续的。虽然这一假定在断点处很难检验，但研究者可以分析在其他位置的回归函数是否连续。当然，非断点处的连续性检验对于断点处的连续性而言既不必要也不充分，但可以应对一部分的质疑。
这一检验通过将断点值改变为其他处理状态未发生改变的点，并估计处理效应实现，预期的结果是在安慰剂断点处并不存在显著的处理效应。使用 Meyersson (2014) 的数据，仍然执行 rdrobust 命令以检验安慰剂断点处的结果。为避免混淆，当选择的安慰剂断点在真实断点左侧时，仅使用控制组，反之仅使用处理组。我们以 为例，说明检验的结果。
* 对应的 Stata 代码
rdrobust Y X if X >= 0, c(1)
# 对应的 R 代码
out = rdrobust(Y[X >= 0], X[X >= 0], c = 1)
. rdrobust Y X if X >= 0, c(1)
Sharp RD estimates using local polynomial regression.
Cutoff c = 1 | Left of c Right of c Number of obs = 317
-------------------+---------------------- BW type = mserd
Number of obs | 32 285 Kernel = Triangular
Eff. Number of obs | 32 51 VCE method = NN
Order est. (p) | 1 1
Order bias (q) | 2 2
BW est. (h) | 2.568 2.568
BW bias (b) | 3.535 3.535
rho (h/b) | 0.726 0.726
Outcome: Y. Running variable: X.
--------------------------------------------------------------------------------
Method | Coef. Std. Err. z P>|z| [95% Conf. Interval]
-------------------+------------------------------------------------------------
Conventional | -2.7759 4.069 -0.6822 0.495 -10.7509 5.19907
Robust | - - 0.2920 0.770 -9.38529 12.6707
--------------------------------------------------------------------------------
无论是传统的推断还是稳健偏误调整的推断均无法拒绝不存在处理效应的原假设，在该安慰剂断点处，回归函数满足连续性假设。除此之外，我们仍然可以使用多种安慰剂断点值，使用下列命令可以结果组织成矩阵形式。
* 对应的 Stata 代码
mat A = J(7,8,.)
local i = 1
forvalues c = -3/3{
if `c'>0{
local condition "if X > 0"
}
else if `c'== 0{
local condition ""
}
else{
local condition "if X < 0"
}
dis "此时选择的断点为 c = `c'"
qui rdrobust Y X,c(`c')
mat A[`i',1] = `c'
mat A[`i',2] = e(h_l)
mat A[`i',3] = e(tau_cl)
mat A[`i',4] = e(pv_cl)
mat A[`i',5] = e(ci_l_cl)
mat A[`i',6] = e(ci_r_cl)
mat A[`i',7] = e(N_h_l)
mat A[`i++',8] = e(N_h_r)
}
mat colnames A = Cutoff Bw Est p-value CI_l CI_r Obs_left Obs_right
matlist A
| Cutoff Bw Est p-value CI_l CI_r Obs_left Obs_right
-------------+----------------------------------------------------------------------------------------
r1 | -3 15.4055 1.776339 .1616983 -.7115872 4.264265 487 303
r2 | -2 11.50318 .329291 .8272219 -2.627645 3.286227 359 234
r3 | -1 14.85515 1.957432 .1683084 -.8273605 4.742224 462 266
r4 | 0 18.03426 2.450962 .0771591 -.2670203 5.168944 551 272
r5 | 1 16.63668 1.849264 .2192406 -1.100955 4.799482 515 239
r6 | 2 15.58472 1.19124 .4709407 -2.04725 4.429729 464 215
r7 | 3 13.93194 1.204791 .5192256 -2.458826 4.868407 413 193
在这个使用了六个不同的断点值的结果中，仅当 是真实断点处的估计结果，作为基准，其他结果均为安慰剂检验结果。可以看到安慰剂断点的结果小于真实断点，并且都不显著，可以作为断点回归有效性的证据。
5.4 断点附近观测值的敏感性
第四种证伪方法用于检验估计结果对距离断点非常近的个体的敏感性。如果存在系统的得分操纵，那么可以很自然地认为，最接近断点的个体极有可能是经过操纵的。该检验的思路在于剔除这些个体后再进行估计和推断，这一方法有时称为 “甜甜圈洞” 方法。即使我们并不怀疑存在操纵得分的个体，同样可以使用该方法剔除不可避免的外推性，因为断点附近的观测值可能会对局部多项式回归的结果有更为明显的影响。
我们同样使用 Meyersson (2014) 的数据，首先使用 的数据进行一次检验，再选取多个半径进行检验。
* 对应的 Stata 代码
rdrobust Y X if abs(X) >= 0.3
# 对应的 R 代码
out = rdrobust(Y[abs(X) >= 0.3], X[abs(X) >= 0.3])
. rdrobust Y X if abs(X) >= 0.3
Sharp RD estimates using local polynomial regression.
Cutoff c = 0 | Left of c Right of c Number of obs = 2619
-------------------+---------------------- BW type = mserd
Number of obs | 2310 309 Kernel = Triangular
Eff. Number of obs | 481 248 VCE method = NN
Order est. (p) | 1 1
Order bias (q) | 2 2
BW est. (h) | 15.901 15.901
BW bias (b) | 27.231 27.231
rho (h/b) | 0.584 0.584
Outcome: Y. Running variable: X.
--------------------------------------------------------------------------------
Method | Coef. Std. Err. z P>|z| [95% Conf. Interval]
-------------------+------------------------------------------------------------
Conventional | 3.1911 1.5139 2.1078 0.035 .223786 6.15834
Robust | - - 1.7741 0.076 -.33272 6.68399
--------------------------------------------------------------------------------
结果表示当我们剔除了断点附近半径为 0.3 的观测值后，结果仍然是稳健的。需要注意的是，虽然剔除部分样本后估计中使用的观测值通常会减少，但其也会受到带宽的影响，因此在估计中使用的观测值可能会增加也可能会减少。与之前的检验相同，我们同样可以使用不同的 “甜甜圈洞” 半径重复上述检验。
* 对应的 Stata 代码
mat A = J(6,7,.)
local i = 1
forvalues r = 0(0.1)0.5{
qui rdrobust Y X if abs(X) >= `r'
mat A[`i',1] = `r'
mat A[`i',2] = e(h_l)
mat A[`i',3] = e(tau_cl)
mat A[`i',4] = e(pv_cl)
mat A[`i',5] = e(ci_l_cl)
mat A[`i',6] = e(ci_r_cl)
mat A[`i++',7] = e(N_h_l)+ e(N_h_r)
}
mat colnames A = Dount-Radius Bw Est p-value CI_l CI_r Obs
matlist A
| Dount-R~s Bw Est p-value CI_l CI_r Obs
-------------+-----------------------------------------------------------------------------
r1 | 0 18.03426 2.450962 .0771591 -.2670203 5.168944 823
r2 | .1 17.79619 2.902978 .0387416 .1502204 5.655735 814
r3 | .2 16.53243 3.136307 .0323339 .2642017 6.008412 764
r4 | .3 15.90105 3.191064 .03505 .2237857 6.158342 729
r5 | .4 17.06812 3.074018 .0324783 .2566045 5.891431 774
r6 | .5 15.47231 3.57846 .016862 .6435131 6.513407 699
5.5 带宽选择的敏感性
最后一种证伪检验方法主要涉及对最优带宽的敏感性。与前一种方法相比，该方法研究的是通过增加或减少区间 “末端” 的点来验证断点回归的有效性。在连续性框架下，该检验方法通过改变用于局部多项式回归的带宽来实现，以说明带宽对结果的影响。
当我们的主要目的是解释点估计结果时，带宽选择的变动最好基于 MSE 最优带宽的小范围变动，否则估计结果可能会严重依赖于估计和推断方法的统计性质。换句话说，带宽大于 MSE 最优带宽过多会导致极大的偏误，而过小则可能降低估计的精度。与之对应的是，如果我们更加关注结果的推断，则带宽的选择应当围绕 CER 最优带宽进行。
* 对应的 Stata 代码
qui rdrobust Y X
local h2 = 2*e(h_l)
local h3 = 3*e(h_l)
rdrobust Y X
rdrobust Y X,h(`h2' `h2')
rdrobust Y X,h(`h3' `h3')
qui rdrobust Y X,bwselect(cerrd)
local h2 = 2*e(h_l)
local h3 = 3*e(h_l)
rdrobust Y X,bwselect(cerrd)
rdrobust Y X,h(`h2' `h2') bwselect(cerrd)
rdrobust Y X,h(`h3' `h3') bwselect(cerrd)
基于 MSE 准则的最优带宽为 ，点估计结果为 2.451，以此为基准，分别估计两倍带宽和三倍带宽下的估计结果，分别为 2.448 和 1.539。如果我们更关注推断结果，则可使用 及其倍数。在实际中，可以将不同带宽下的估计结果绘制在图中，以展示断点回归设计的有效性。
6. 参考文献
Meyersson E. Islamic Rule and the Empowerment of the Poor and Pious[J]. Econometrica, 2014, 82(1): 229-269. -PDF- Cattaneo M D, Idrobo N, Titiunik R. A practical introduction to regression discontinuity designs: Foundations[M]. Cambridge University Press, 2019. -PDF-
7. 相关推文
Note：产生如下推文列表的 Stata 命令为：
lianxh rdd, m
安装最新版lianxh命令：
ssc install lianxh, replace
温馨提示： 文中链接在微信中无法生效。请点击底部「阅读原文」。
专题：断点回归RDD RDD：离散变量可以作为断点回归的分配变量吗？ rddensity, lpdensity无法安装？那就手动安装 RDD：断点回归可以加入控制变量吗？ 断点回归RDD：样本少时如何做？ Stata：断点回归分析-RDD-文献和命令 Stata：两本断点回归分析-RDD-易懂教程 Stata：RDD-中可以加入控制变量 Stata：时间断点回归RDD的几个要点 Stata：断点回归分析-(RDD)-文献和命令 Stata：断点回归RDD简明教程 RDD：断点回归的非参数估计及Stata实现 Stata: 两本断点回归分析 (RDD) 易懂教程 Stata: 断点回归 (RDD) 中的平滑性检验 Stata 新命令：多断点 RDD 分析 - rdmc RDD 最新进展：多断点 RDD、多分配变量 RDD 专题：内生性-因果推断 Abadie新作：简明IV,DID,RDD教程和综述
🍓 课程推荐：2024 机器学习与因果推断专题
主讲老师：司继春 (上海对外经贸大学) ；张宏亮（浙江大学）
课程时间：2024 年 11 月 9-10 日 ；16-17日
课程咨询：王老师 18903405450（微信）
尊敬的老师 / 亲爱的同学们：
连享会致力于不断优化和丰富课程内容，以确保每位学员都能获得最有价值的学习体验。为了更精准地满足您的学习需求，我们诚挚地邀请您参与到我们的课程规划中来。 请您在下面的问卷中，分享您 感兴趣的学习主题或您希望深入了解的知识领域 。您的每一条建议都是我们宝贵的资源，将直接影响到我们课程的改进和创新。 我们期待您的反馈，因为您的参与和支持是我们不断前进的动力。感谢您抽出宝贵时间，与我们共同塑造更加精彩的学习旅程！https://www.wjx.cn/vm/YgPfdsJ.aspx# 再次感谢大家宝贵的意见！
New！ Stata 搜索神器：
lianxh和songblGIF 动图介绍
搜： 推文、数据分享、期刊论文、重现代码 ……
👉 安装：
. ssc install lianxh
. ssc install songbl
👉 使用：
. lianxh DID 倍分法
. songbl all
🍏 关于我们
连享会 ( www.lianxh.cn，推文列表) 由中山大学连玉君老师团队创办，定期分享实证分析经验。 直通车： 👉【百度一下： 连享会】即可直达连享会主页。亦可进一步添加 「知乎」,「b 站」,「面板数据」,「公开课」 等关键词细化搜索。
