Stata+R：一文读懂精确断点回归-RDD

发布时间：2021-03-28 阅读 17204

Stata连享会主页 || 视频 || 推文 || 知乎

温馨提示： 定期清理浏览器缓存，可以获得最佳浏览体验。

New！ lianxh 命令发布了：
随时搜索推文、Stata 资源。安装命令如下：
. ssc install lianxh
详情参见帮助文件 (有惊喜)：
. help lianxh

✌ 课程详情： https://gitee.com/lianxh/Course

⛳ 课程主页： https://gitee.com/lianxh/Course

⛳ Stata 系列推文：

全部 | Stata入门 | Stata教程 | Stata资源 | Stata命令
计量专题 | 论文写作 | 数据分享 | 专题课程
结果输出 | Stata绘图 | 数据处理 | Stata程序
回归分析 | 面板数据 | 交乘项-调节 | IV-GMM
内生性-因果推断 | 倍分法DID | 断点回归RDD | PSM-Matching | 合成控制法
Probit-Logit | 时间序列 | 空间计量 | 分位数回归 | 生存分析 | SFA-DEA
文本分析-爬虫 | Python-R-Matlab | 机器学习
Markdown | 工具软件 | 其它

作者： 曹昊煜 (兰州大学)
邮箱： 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. 简介

社会科学家主要对因果关系感兴趣。这种关系在随机试验中很容易得到，但在社会科学中，随机试验会受到道德和实践等诸多因素的限制。在缺少随机试验的前提下，我们往往诉诸于严格的准自然实验干预。其中，RD 方法就是一种最为可信的因果推断分析方法。

在 RD 方法中，所有单位都拥有一个 “得分”，并且得分值高于某个断点的个体接受处理。这一设计的最大特点在于接受处理的概率在门槛值处急剧变化。具体来看，RD 的特征如下：

RD 具有三个明显的特征：得分、断点和处理；
RD 的主要特征是处理分组在断点处为得分的分段函数，这一条件是可以直接验证的，除此之外，还有很多证伪检验和相关经验来提高可信度；
除了上述特征，一般的解释、估计和推断也需要一些额外的假设。首先，我们需要定义参数及其识别条件；其次，必须施加假设以保证参数可估计，主要的框架有两类，分别是连续性假设和局部随机化假设。

2. 精确断点回归设计

2.1 RD 相关理论介绍

在 RD 中，每个个体拥有一个得分 (或称为驱动变量 running variable，或指标 index)，并且处理取决于得分是否高于断点。

引入一些符号，假设有 $n$ 个个体，使用 $i = 1, 2, 3, \dots, n$ 标注，每个个体的得分表示为 $X_{i}$ ， $c$ 是一个已知的断点。当 $X_{i} \geq c$ 时，个体接受处理，否则不接受处理。处理状态可以表示为 $T_{i} = I (X_{i} \geq c)$ ，该形式为示性函数。因此，知道个体得分，就知道个体处理状态，即处理的概率是得分的函数。

但是满足处理条件并不意味着实际受到了处理。当处理条件与实际受到处理的情形一致时，我们称之为精确断点回归 (Sharp RD) 设计，否则称之为模糊断点回归 (Fuzzy RD)。图 1 描述了精确断点回归设计的处理规则，即在断点附近，受到处理的概率从零跳跃至 1。

在使用潜在结果框架分析处理效应时，通常假设个体拥有两个潜在结果 $Y_{i} (1)$ 和 $Y_{i} (0)$ 。在这一框架下，处理效应定义为潜在结果特征或者分布的比较，例如均值、方差或者分位数，但二者中往往只有一个可以被观测。其中，可观测结果变量表示为：

Y_{i} = (1 - T_{i}) \cdot Y_{i} (0) + T_{i} \cdot Y_{i} (1) = {\begin{aligned} Y_{i} (0) & X_{i} < c \\ Y_{i} (1) & X_{i} \geq c \end{aligned}

图 2 描述了精确断点回归中的处理效应。在该图中，绘制了给定得分水平时的平均潜在结果，表示为 $E (Y_{i} (1) | X_{i} = x)$ 和 $E (Y_{i} (0) | X_{i} = x)$ 。

从图中可以看出，平均处理效应是给定任意得分时，回归曲线的垂直距离，但由于无法观测到同一得分值处的两条曲线，所以这一距离无法直接估计。但在 $c$ 点处是个例外，我们可以同时观测到两条曲线。此时，平均处理效应可以定义为：

τ_{S R D} = E [Y_{i} (1) - Y_{i} (0) | X_{i} = c]

需要注意的是，我们捕获的信息是 $X_{i} = c$ 处的处理效应。换句话说， $τ_{S R D}$ 只能视为局部处理效应。

对得分值相近但分布在断点两侧的个体施加可比性假定，是 RD 设计的基础，这一思想主要由 Hahn 等 (2001) 的连续性假定刻画。Hahn 等 (2001) 认为，如果关于 $x$ 的回归函数 $E [Y_{i} (1) | X_{i} = x]$ 和 $E [Y_{i} (0) | X_{i} = x]$ 在 $x = c$ 处连续，那么在清晰断点中，有：

E [Y_{i} (1) - Y_{i} (0) | X_{i} = c] = lim_{x \to c^{+}} E [Y_{i} | X_{i} = c] - lim_{x \to c^{-}} E [Y_{i} | X_{i} = c]

2.2 RD 效应的局部性质

与其他方法相比，精确断点回归仅仅在驱动变量的一个点上计算了平均差异，也就是局部的因果关系。换言之，其处理效应只包含有限的外部有效性，即精确断点回归设计的结果对远离断点的个体没有良好的代表性。

一般来讲，断点回归估计结果的外部有效性或代表性取决于特殊的应用。举例来说，在图 3 的左图中， $E [Y_{i} (1) | X_{i} = x]$ 和 $E [Y_{i} (0) | X_{i} = x]$ 之间的垂直距离在 $x = c$ 处明显为正，并且平均潜在结果的差异基本上处处为正。但右图在断点处的差异为 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)”。

为了更细致的定义这些箱体，我们假设驱动变量在 $[x_{l}, x_{u}]$ 之间取值。我们使用正号和负号下标表示处理组和控制组，分别构造不相交的区间，并使用 $J_{-}$ ， $J_{+}$ 表示左右箱体的总数。定义：

Control Bins : B_{-, j} = {\begin{aligned} [x_{l}, b_{-, j}) & j = 1 \\ [b_{-, j - 1}, b_{-, j}) & j = 2, \dots, J_{-} - 1 \\ [b_{-, j - 1}, c) & j = J_{-} \end{aligned}

Treated Bins : B_{+, j} = {\begin{aligned} [c, b_{+, j}) & j = 1 \\ [b_{+, j - 1}, b_{+, j}) & j = 2, \dots, J_{+} - 1 \\ [b_{+, j - 1}, x_{u}) & j = J_{+} \end{aligned}

上述的定义方式说明 $B_{-, 1} \cup B_{-, 2} \cup, \dots, B_{J -}, B_{+, 1} \cup B_{+, 2}, \dots, B_{J_{+}}$ 划分了驱动变量的范围 $[x_{l}, x_{u}]$ ，其中心为断点处。

在进一步定义区间之前，我们首先约定 $X_{-, i}$ 和 $X_{+, i}$ 表示控制组和处理组子区间内的样本个数， $⌊ \cdot ⌋$ 为取整函数。接下来我们可以分别定义两种类型的箱体。

Evenly-spaced（ES）箱体：区分驱动变量取值范围的不相交区间，所有区间长度相同；

b_{-, j} = x_{l} + \frac{j (c - x_{l})}{J_{-}}, b_{+, j} = c + \frac{j (x_{u} - c)}{J_{+}}

Quantile-spaced（QS）箱体：区分驱动变量取值范围的不相交区间，每个区间大致包含相同数量的观测个数。

b_{-, j} = X_{-, ⌊ j / J_{-}} ⌋, b_{+, j} = X_{+, ⌊ j / J_{+}} ⌋

两种箱体具有一定的差异性。尽管 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 的方式选定后，还需要确定断点两侧的箱体数量。接下来我们介绍两种在给定箱体类型后，以数据驱动，自动进行的选择 $J_{-}$ 和 $J_{+}$ 方法。

第一个方法通过最小化局部均值的积分均方误差 (Integrated mean-squared error, IMSE) 来选择断点两侧的箱体数量。积分均方误差定义为方差和平方误差的期望和。如果我们选择较大的箱体个数，偏误将会很小，因为区间范围更小并且局部常数拟合也会更好，这一方式的缺陷在于每个箱体内的观测值非常少，导致箱体内的可变性非常大。IMSE 方法的思路正在于平衡误差平方和变异。

选择 IMSE 最优的箱体将得到 “追踪” 潜在回归函数的均值，这意味着 IMSE 最优的箱体更有利于刻画回归函数的形状。然而，IMSE 方法通常在局部均值与全局多项式拟合几乎重合时会绘制出非常平滑的图片，在在断点附近难以捕捉数据的局部变异。

IMSE 最优的 $J_{-}$ 和 $J_{+}$ 分别可以做如下的定义：

J_{-}^{I M S E} = ⌈ L_{-}^{I M S E} n^{1 / 3} ⌉, J_{+}^{I M S E} = ⌈ L_{+}^{I M S E} n^{1 / 3} ⌉

这里的 $n$ 是观测值得总数， $⌈ \cdot ⌉$ 为取顶函数， $L_{-}^{I M S E}$ 和 $L_{+}^{I M S E}$ 取决于 ES 或者 QS 箱体类型。实践中，未知常数 $L_{-}^{I M S E}$ 和 $L_{+}^{I M S E}$ 使用原始的、客观的数据驱动方法估计。

在软件中，可以使用 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 箱体，并且断点两侧的区间长度有所不同，其原因在于算法需要保证每个箱体内拥有一样多的观测值。

第二种选择 $J_{-}$ 和 $J_{+}$ 的方法称为 Minmicking Variance 方法，该方法中箱体数量的选择依据是使区间均值的总体变化等于原始数据散点图的变化。MV 方法选择的箱体数量定义为：

J_{-}^{M V} = ⌈ L_{-}^{M V} \frac{n}{l o g (n)^{2}} ⌉, J_{+}^{M V} = ⌈ L_{+}^{M V} \frac{n}{l o g (n)^{2}} ⌉

公式中的 $n$ 为样本数量，常数的确定取决于 ES 或者 QS 箱体类型，并且其大小与 IMSE 最优的选择不同。一般而言， $J_{-}^{M V} > J_{-}^{E S}$ 并且有 $J_{+}^{M V} > J_{+}^{E S}$ ，在软件中使用如下方式使用 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. 基于连续性假定的断点回归分析

本小节中讨论的方法主要基于连续性假定，或者说：

E [Y_{i} (1) - Y_{i} (0) | X_{i} = c] = lim_{x \to c^{+}} E [Y_{i} | X_{i} = c] - lim_{x \to c^{-}} E [Y_{i} | X_{i} = c]

在这一框架下，将 $τ_{S R D}$ 定义为感兴趣的参数，估计方法是在断点两侧使用多项式方法近似回归函数。实践中，我们往往使用最小二乘法拟合结果变量关于驱动变量的多项式函数。当我们使用所有观测值时，这种多项式拟合是全局的，但是，如果我们仅适用断点附近的观测值拟合多项式时，使用的是局部方法或非参数方法。接下来的讨论中，我们主要使用的是局部方法。

4.1 局部多项式的点估计

局部多项式方法采用线性回归，但仅仅使用断点附近的观测值。具体而言，该方法关注的观测值处于 $c - h$ 和 $c + h$ 区间内，其中带宽 $h > 0$ 。在带宽内，通常使用加权方法保证临近断点的观测值相对于其他观测具有更高的权重，这一过程主要通过核函数实现 $K (\cdot)$ 。局部多项式估计由以下步骤组成：

选择多项式阶数 $q$ 和核函数 $K (\cdot)$ ；
选择带宽 $h$ ；
对断点右侧的观测值，结果变量 $Y_{i}$ 使用加权最小二乘法拟合。解释变量包括常数项和 $(X_{i} - c), (X_{i} - c)^{2}, \dots, (X_{i} - c)^{p}$ ，其中 $p$ 为选定的多项式阶数，权重为 $K (\frac{X_{i} - c}{h})$ 。截距的估计结果 ${\hat{μ}}_{+}$ 正是 $μ_{+} = E [Y_{i} (1) | X_{i} = c]$ 的估计:

{\hat{μ}}_{+} : {\hat{Y}}_{i} = {\hat{μ}}_{+} + {\hat{μ}}_{+, 1} (X_{i} - c) + {\hat{μ}}_{+, 2} (X_{i} - c)^{2} + \dots + {\hat{μ}}_{+, p} (X_{i} - c)^{p}

对断点左侧的观测值使用同样的方式，截距的估计结果 ${\hat{μ}}_{-}$ 正是 $μ_{-} = E [Y_{i} (0) | X_{i} = c]$ 的估计：

{\hat{μ}}_{-} : {\hat{Y}}_{i} = {\hat{μ}}_{-} + {\hat{μ}}_{-, 1} (X_{i} - c) + {\hat{μ}}_{-, 2} (X_{i} - c)^{2} + \dots + {\hat{μ}}_{-, p} (X_{i} - c)^{p}

计算精确断点的点估计： ${\hat{τ}}_{S R D} = {\hat{μ}}_{+} - {\hat{μ}}_{-}$ 。

选择核函数和多项式阶数

核函数 $K (\cdot)$ 为每一个观测值指派非负的权重，计算依据是观测值对应的驱动变量值与断点的距离 $\frac{X_{i} - c}{h}$ 。一般推荐使用三角核函数进行加权，其函数形式为：

K (u) = (1 - | u |) I (| u | \leq 1)

其原因是，在最优均方误差带宽条件下，三角核函数可以得到均方误差意义下的最优点估计。虽然三角核函数具有优良的性质，但有时我们仍然可能使用均匀核函数，其函数形式为 $K (u) = I (| u | \leq 1)$ ，带宽以外的权重仍为 0，但区间内的权重全部相同，此时的加权最小二乘法等同于线性回归。

除了核函数外，多项式阶数的选取也会影响估计结果。首先，如果选择 0 阶多项式，也就是对常数拟合，在边界处无法得到优良的理论性质。其次，对于给定带宽，增加多项式的阶数不但可以改变近似的程度，也可以增加处理效应的可变性。实践中，0 阶多项式和线性多项式拟合几乎相同，但后者拥有更小的渐进偏误。最后，高阶多项式可能产生过拟合导致结果不可靠。

带宽选择及其实现

带宽 $h$ 用于控制局部多项式回归使用的样本情况，带宽选择是断点回归分析的基础，因为其直接影响者估计和推断的性质，并且经验分析对带宽选择非常敏感。

选择一个更小的带宽可以减少模型误设的偏误，但同时有可能提高参数估计的方差，因为带宽越窄，使用的样本越少。反过来，使用更大的带宽可以提高估计的精度，但可能提高模型设定的误差。因此，带宽选择又称为 “偏误-方差权衡”。

由于带宽选择非常重要，所以我们需要使用一些数据驱动的方法来避免数据挖掘。多数带宽选择方法正是试图平衡偏误和方法。实践中使用最为广泛的方式即为最小化局部多项式估计的均方误差 (MSE)。断点回归处理效应渐进均方误差的通常形式为：

M S E ({\hat{τ}}_{S R D}) = B i a s^{2} ({\hat{τ}}_{S R D}) + V a r i a n c e ({\hat{τ}}_{S R D}) = B^{2} + V

进一步定义渐进偏误和估计方差为：

B = h^{2 (p + 1)} B ， V = \frac{1}{n h} V

标量 $B$ 和 $V$ 分别表示处理效应点估计的偏误和方差，其中使用了样本数量和带宽做调整。尽管省略了很多技术细节，我们展示了偏误和方差的一般形式以澄清局部多项式估计中包含最优均方误差选择的关键权衡。

偏误的一般形式由 $h^{2 (p + 1)}$ 决定，同时：

B = B_{+} - B_{-} ， B_{-} = μ_{-}^{p + 1} b_{-} ， B_{+} = μ_{+}^{p + 1} b_{+}

μ_{+}^{p + 1} = lim_{x \to c^{+}} \frac{d^{p + 1} E [Y_{i} (1) | X = x]}{d x^{p + 1}} ， μ_{-}^{p + 1} = lim_{x \to c^{-}} \frac{d^{p + 1} E [Y_{i} (0) | X = x]}{d x^{p + 1}}

以上的两个极限与未知函数的斜率有关，常数 $B_{+}$ 和 $B_{-}$ 与核函数和多项式阶数有关，该计算方法假设了共同的带宽 $h$ ，但该表达式也可以扩展到断点左右的不同带宽。

偏误项取决于回归函数的 $(p + 1)$ 阶导数。当我们近似 $E [Y_{i} (1) | X = x]$ 和 $E [Y_{i} (0) | X = x]$ 时，该近似必存在一个误差。例如我们使用局部线性模型逼近时，误差的主要部分取决于未知函数的二阶导数。

方差 $V$ 取决于样本大小和带宽，也可做如下的定义：

V = V_{-} + V_{+} ， V_{-} = \frac{σ_{-}^{2}}{f} v_{-} ， V_{+} = \frac{σ_{+}^{2}}{f} v_{+}

σ_{+}^{2} = lim_{x \to c^{+}} V [Y_{i} (1) | X_{i} = x] ， σ_{-}^{2} = lim_{x \to c^{-}} V [Y_{i} (0) | X_{i} = x]

该定义与结果变量在断点处的条件变动有关， $f$ 为断点处驱动变量的密度函数，常数 $v_{-}$ 和 $v_{+}$ 与和函数和多项式阶数有关。

为得到最优的 MSE 点估计 ${\hat{τ}}_{S R D}$ 我们选择最小化 MSE 的带宽：

m i n_{h > 0} [h^{2 (p + 1)} B \dots \dots 2 + \frac{1}{n h} V]

依据该准则，我们可以得到：