RDD：离散变量可以作为断点回归的分配变量吗？

发布时间：2021-03-03 阅读 3983

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

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

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

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

✨ 课程主页： https://gitee.com/arlionn/Course

Stata 系列推文：

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

作者：邱紫烨（中山大学）
邮箱：qiuzy@mail2.sysu.edu.cn

编者按：本文主要参考自下文，特此致谢！
Source: Kolesár M, Rothe C. Inference in regression discontinuity designs with a discrete running variable[J]. American Economic Review, 2018, 108(8): 2277-2304. -PDF- R-Data-Prog

1. 背景介绍

谢谦等 (2019) 对目前学术界断点回归 (regression-discontinuity designs) 的应用进行了综述，但他们侧重于经济学的五大顶级期刊上已使用的 RD 方法，对于未出现的多配置变量 RDD、分位数 RDD、拐点回归设计、多断点 ED、以及离散型配置变量 RDD 等方法并未作详细介绍。为补充相关内容，我们将在本文中简要介绍 RDD 的最新进展——离散型变量配置 RDD。

2. RDD 理论回顾以及存在问题

断点回归是一种随机试验，其接受处置的概率是一个或者几个变量的间断函数。断点回归可以分为两类：

一是拥有确定的临界值 (Sharp RDD)，即在临界值一侧的所有观测点都接受了处置，反之，在临界值另一侧的所有观测点都没有接受处置；
二是临界点是模糊的 (Fuzzy RDD)，即在临界值附近，接受处置的概率是单调变化的。

Hahn 等 (2001) 证明了在一定的假设下，无论是哪一类型的断点回归，都可以利用临界值附近样本的系统性变化来研究处置与否和经济变量变动的因果关系，并提出了相应的估计方法。由此，断点回归在经济学中的应用逐渐普及开来。

估计 ATE (Average Treatment Effect) 的经典方法是多项式回归，即选取一定的带宽，并对临界值两边窗口内的经济变量进行线性回归，根据临界值附近经济变量值的跳跃估计平均处置效应。然而，由于真实的条件期望函数往往不是线性的，因而估计出的 ATE 往往是有偏误的。如果选取的带宽足够小，这种估计偏误可以忽略不计，从而可以使用基于 EHW 异方差稳健标准误的置信区间进行因果推断。

反之，由于变量间的数值差距过大，研究者就需要选择一个大的带宽以获得足够多的观察值，此时估计偏误就不能忽略不计了。例如学生的分数、学校的入学率、公司的员工数量以及个人出生年份等。现有研究主要通过对变量进行聚类来解决这个问题，即将值相同的变量聚类，并通过构造类稳健的标准误来估计处置效应的方差。

但是，在离散型变量的情形中，这种聚类的方法可能无法解决 ATE 的估计偏误问题，具体来看：

根据 2003 年至 2005 年人口普查数据，以 40 岁为临界值，绘制年龄 (横轴) 与时薪取对数 (纵轴) 之间的关系。可以发现，在临界值两端的时薪取值并没有明显的跳跃，说明平均处置效应不存在。接着，对人口普查数据进行小样本随机抽样 (样本量分别取 100、500、2000、10000)，并且每个子样本都按下式进行 OLS 估计：

l o g (W A G E_{i}) = α_{h} + τ_{h} \cdot I {A G E \geq 40} + \sum_{j = 1}^{p} β_{h j}^{-} \cdot (A G E_{i} - 40)^{p} + \sum_{j = 1}^{p} β_{h j}^{+} \cdot I {A G E_{i} \geq 40} \cdot (A G E_{i} - 40)^{p} + U_{i}

通过估计上式，可以得出临界值附近的处置效应。又由上图可知，临界值两端的经济变量并不是完全的线性的，故 OLS 估计必然存在一定的偏误，且随着带宽的变动而发生改变。

为了进一步评估系数估计的准确性，原文还计算了每一个子样本的 CRV 和 EHW 标准误。由于模型存在估计偏误的问题，因此 EHW 标准误构造的置信区间只能在小于 95% 的概率下捕捉到这种零值效应，而当方差足够大时，这种估计偏误可以忽略不计。此外，如果将变量聚类是合适的修正偏误的办法，那么 CRV 标准误应当大于 EHW 标准误，并且由 CRV 标准误构造的置信区间应该在大约 95% 的概率水平上捕捉到零值效应。

下表给出了估计的结果，在一次项回归中，保持样本量不变，系数估计存在向下偏误，且 CRV 标准误会随着带宽的增加而逐渐扩大。同时，尽管 CRV 标准误均大于 EHW 标准误，其所构造的置信区间在大样本情况下效果不佳。而在二次项回归中的情形中，保持样本量不变，无论选取的带宽大小如何，CRV标准误均较小，由其构造的置信区间都在 95% 的显著水平之下。

当所估计的模型拥有较小的估计偏误时，标准误也往往较小，但会导致其构造的置信区间相应的缩小，而当聚类数量较小时，这种问题尤为严重。另一方面，当样本量的选取不确定时，往往存在比较严重的估计偏误现象，但此时的置信区间依然很狭窄。

3. 数理证明

在清晰断点回归 (Sharp RDD) 中，从大样本中随机抽取 $N$ 个观察值， $Y_{i} (1)$ 和 $Y_{i} (0)$ 分别为接受和没有接受处置的观察值因变量， $T_{i} \in {0, 1}$ 为观察值是否接受处置的指示变量，且只有当变量达到经过正态化处理的临界值才会被标记为受到处置，即

T_{i} = I {X_{i} \geq 0}

若用 $μ (X_{i}) = E [Y_{i} | X_{i}]$ 表示期望值，平均处置效应 (ATE) 的表达式如下：

μ (X_{i}) = E (Y_{i} (1) - Y_{i} (0) | X_{i} = 0) = lim_{x ↓ 0} μ (x) - lim_{x ↑ 0} μ (x)

用来估计 ATE 的方法主要是局部多项式回归。其步骤为先确定带宽 $h > 0$ 和项数 $p \geq 0$ ，其中， $p$ 一般取 1 或者 2。然后剔除带宽之外的所有样本，对 $[- h, h]$ 内的样本进行 $X_{i}$ 对 $Y_{i}$ 的OLS 回归，各个系数的数学表达式如下所示：

\hat{τ} = e_{1}^{'} \hat{θ}, θ = {\hat{Q}}^{- 1} \frac{1}{N_{h}} \sum_{i = 1}^{N_{h}} M_{i} M_{i}^{'}, M_{i} = m (X_{i})

m (x) = (I {x \geq 0}, I {x \geq 0} x, . . ., I {x \geq 0} x^{p}, 1, x, . . ., x^{p})^{'}

在系数估计的基础上，得出 $\hat{τ}$ 的异方差稳健标准误 ${\hat{σ}}_{E H W} / {\sqrt{N}}_{h}$ 。为了解决离散型断点回归中 EHW 标准误可能存在的问题，引入了 CRV 标准误，并将拥有相同结果的经济变量归为同一类，具体表达式为：

{\hat{σ}}_{C R V} / {\sqrt{N}}_{h}

基于 EHW 标准误和 CRV 标准误分别构建出置信区间：

\hat{τ} \pm z_{1 - α / 2} \times {\hat{σ}}_{E H W} / {\sqrt{N}}_{h}, a n d \hat{τ} \pm z_{1 - α / 2} \times {\hat{σ}}_{C R V} / {\sqrt{N}}_{h}

那么，在离散型变量的情况下，基于 EHW 标准误构造的置信区间会产生什么问题呢？在有限样本中， $\hat{τ}$ 是 $τ_{h}$ 的渐近无偏估计：

τ_{h} = e_{1}^{'} θ_{h}, θ_{h} = Q^{- 1} E [M_{i} Y_{i} | | X_{i} | \leq h], Q = E [M_{i} M_{i}^{'} | | X_{i} | \leq h]

然而， $\hat{τ}$ 往往是 $τ$ 的有偏估计，估计偏误 $τ - \hat{τ}$ 的大小取决于 $p$ 次多项式能够对于 $μ (x)$ 的拟合程度。在大样本情形下，EHW 服从标准正态分布，即：

\sqrt{N_{h}} (\hat{τ} - τ_{h}) / {\hat{σ}}_{E H W} ⟶^{d} N (0, 1)

因此，基于 EHW 标准误构造的 $t$ 统计量和置信区间能够用来估计 $τ_{h}$ ，如果估计偏误 $τ_{h} - τ$ 能够依样本趋近于 0，即：

\frac{τ_{h} - τ}{{\hat{σ}}_{E H W} / \sqrt{N_{h}}} \to^{P} 0

那么就能够用这个方法估计 $τ$ ，且无论在连续型还是离散型变量的情况下都适用。通过选取一个合适的带宽 $h$ 使得上式成立的过程就是欠光滑过程。但是，如果临界值两端的变量值差距过大，为了拥有足够的样本，研究者不得不选取较大的带宽，使得最后的估计值出现偏误。

为了解决 EHW 标准误在离散型 RDD 中的估计偏误问题，引入了 CRV 标准误，并在随后被广泛使用。首先，令 $δ (x) = μ (x) - m (x)^{'} θ_{h}$ 代表 $p$ 次多项式对条件期望的估计偏误，并且令 $δ_{i} = δ (X_{i})$ 。接着，对于估计窗口内的变量，有：

Y_{i} = M_{i}^{'} θ_{h} + u_{i}, u_{i} = δ_{i} + ε_{i}

其中， $ε_{i} = Y_{i} - μ (X_{i})$ 为真实值与期望值之差，且所有实现值相同的变量有相同的 $δ_{i}$ 。此时， $δ_{i}$ 就变成了一个随机效应，而不是根据误差而发生改变的具体值，即：

δ_{i} = \sum_{g = 1}^{G_{h}} D_{g} I {X_{i} = x_{g}},

其中， $D = (D_{1}, . . ., D_{h})$ 是一列独立同分布的零均值随机向量。

4. 置信区间的性质

通过讨论 ${\hat{σ}}_{E H W}^{2}$ 的性质，我们能够了解其对应的置信区间的性质。在第 3 部分讨论基础上，存在常数 $b$ 和 $d_{1}$ ，使得当样本逐渐增加时，对每一个在窗口内的变量 $x_{g}$ ，有：

\sqrt{N_{h}} (τ_{h} - τ) \to b a n d \sqrt{N_{h}} δ (x_{g}) \to d_{g}

在大样本情形下， $\sqrt{N_{h}} (\hat{τ} - τ)$ 服从均值为零的正态分布，渐近方差为 $σ_{τ}^{2} = \sum_{g = 1}^{G_{h}} σ_{g}^{2} ω_{g}$ 。其中， $ω_{g} = e_{1}^{'} Q^{- 1} Q_{g} Q^{- 1} e_{1}$ 。EHW 方差估计值是 $σ_{τ}^{2}$ 的一致估计。

正如第 3 部分中所讨论的，在 $b$ 不为 0 的情况下，利用 EHW 标准误构造的置信区间往往会低估 $τ$ ，而利用 CRV 标准误 ${\hat{σ}}_{C R V} / \sqrt{N_{h}}$ 构造的置信区间也有相同的点估计，因此，为了让CRV置信区间有正确的覆盖范围，理论上 CRV 标准误平均值应该比 EHW 标准误平均值更大一些。

然而，事实也并非如此。在正确识别的情况下，利用 CRV 标准误构造的置信区间往往会低估平均处置效应，而在错误识别的情况下，存在低估和高估两种可能。 ${\hat{σ}}_{C R V}^{2}$ 的渐近性质取决于估计矩内支持点数量是固定不变的还是会随着样本大小而发生改变。

原文检验了在样本增加情况下保持估计窗口不变 ${\hat{σ}}_{C R V}^{2}$ 的性质，并以此来模拟估计窗口 $[- h, h]$ 内样本数量较少的情形。实证结果表明， ${\hat{σ}}_{C R V}$ 并没有收敛到一个常数，而是收敛到一个非退化极限 ${\hat{σ}}_{C R V, \infty}^{2}$ ，这意味着即使在大样本情形下它依然是随机值。根据定理 1，我们能够将 ${\hat{σ}}_{C R V, \infty}^{2}$ 与 $σ_{τ}^{2}$ 之间的期望差异分解为如下形式：

E ({\hat{σ}}_{C R V, \infty}^{2}) - σ_{τ}^{2} = \sum_{g = 1}^{G_{h}} d_{g}^{2} π_{g} ω_{g} + \sum_{g = 1}^{G_{h}} m {(x_{g})}^{'} Q^{- 1} (\sum_{j = 1}^{G_{h}} σ_{j}^{2} Q_{j} - 2 σ_{g}^{2} Q) Q^{- 1} m (x_{g}) \cdot π_{g} ω_{g}

等式右边第一项是正值，大小取决于错误识别的程度。第二项与错误识别的程度无关，当将估计残差 ${\hat{μ}}_{i}$ 用真实残差值 $μ_{i}$ 替代时会使第二项变为 0。第二项的符号往往难以判断，因此，需要假设条件方差服从渐近同分布，即 $σ_{g}^{2} = σ^{2}$ ，从而使得 $σ_{τ}^{2} = σ^{2} \sum_{g = 1}^{G_{h}} ω_{g}$ ，将化简后的结果除以 $σ_{τ}^{2}$ ，可以得到期望差值的变化率：

\frac{E ({\hat{σ}}_{C R V, \infty}^{2}) - σ_{τ}^{2}}{σ_{τ}^{2}} = \frac{\sum_{g = 1}^{G_{h}} d_{g}^{2} π_{g} ω_{g}}{σ^{2} \sum_{g = 1}^{G_{h}} ω_{g}} - \frac{\sum_{g = 1}^{G_{h}} m {(x_{g})}^{'} Q^{- 1} m (x_{g}) π_{g} ω_{g}}{\sum_{g = 1}^{G_{h}} ω_{g}} \equiv T_{1} + T_{2}

由于 EHW 方差估计量是 $σ_{τ}^{2}$ 的一致估计，根据上式，可以推断出 ${\hat{σ}}_{C R V}^{2}$ 的值大致是 ${\hat{σ}}_{E H W}^{2}$ 的 $(1 + T_{1} + T_{2})$ 倍。为了使 CI 置信区间能够很好地反映错误识别的程度， $(T_{1} + T_{2})$ 应该为正数，但事实往往并非如此。

首先， $T_{1}$ 项确实是正值，且会随着估计偏误程度的增大而增加。而 $T_{2}$ 项在同方差假设下是负值，因为 $m {(x_{g})}^{'} Q^{- 1} m (x_{g})$ 和 $ω_{g}$ 都是正值，这一假设在轻微异方差的情况下同样成立。而 $T_{2}$ 的大小只与连续变量的边缘分布有关，因此，如果误差项的方差 $σ^{2}$ 非常大的话， $T_{2}$ 就会成为整个式子正负符号的决定项。

此外， $T_{2}$ 还与最高次项系数 $p$ 有关。根据下式可知， $T_{2}$ 项是对 $2 (p + 1)$ 的再次加权平均，次项数 $p$ 越高 $T_{2}$ 也会随之增加，这也就意味着 $p$ 的增加会增大 ${\hat{σ}}_{C R V, \infty}^{2}$ 与 $σ_{τ}^{2}$ 之间的期望差异。

\sum_{g = 1}^{G_{h}} m {(x_{g})}^{'} Q^{- 1} m (x_{g}) π_{g} = trace (\sum_{g = 1}^{G_{h}} Q_{g} Q^{- 1}) = 2 (p + 1)

最后， $T_{2}$ 与小样本估计下的衰减偏误有一致的趋势。尽管理论上这种估计偏误非常普遍，但它对于小样本估计有很大影响。

总之，理论结果显示，如果错误识别的程度和聚类个数都可以忽略不计的话，CRV 标准误的值会比 EHW 标准误的值更小些。因此，使用连续变量聚类会加剧而不是解决估计偏误的问题，这种情况在正确识别的情形下尤为严重。

接着，原文考察了随着带宽 $G_{h}$ 和样本量的增加， ${\hat{σ}}_{C R V}^{2}$ 的性质会发生什么样的变化，并以此来模拟估计窗口 $[- h, h]$ 内有大量样本点的情形。随着 $G_{h}$ 和 $N_{h}$ 的无限增大， ${\hat{σ}}_{C R V}^{2}$ 会逐渐趋 0，由此我们得出：

{\hat{σ}}_{C R V}^{2} - σ_{τ}^{2} = \sum_{g = 1}^{G_{h}} d_{g}^{2} π_{g} ω_{g} + o_{P} (1)

这意味着在样本量足够大的情况下，CRV 标准误确实比 EHW 标准误更大些，此时对连续变量进行聚类有更好的估计效果。但这并不意味着 CRV 置信区间能够对所有的期望误差都有正确的估计，只是表明在估计窗口 $[- h, h]$ 内存在这样一个期望函数集 $M_{C R V}$ 能够得出一致的估计。

但这一函数集包含的内容要远远不止多项式方程，所以这并不能成为使用 CRV 标准误的理由。事实上，任何使得 CRV 置信区间宽于 EHW 置信区间的函数形式 (包括加入一个随机的常数项) 都能带来更好的估计效果。所以，判断基于 CRV 标准误的估计结果是否稳健要看 $M_{C R V}$ 中是否有清晰易解释的条件期望函数，以帮助对 $μ (x)$ 做出合理的假设。

然而， $M_{C R V}$ 的范围受到 $ω_{g}$ 和 $π_{g}$ 的影响，而这取决于连续变量的分布情况。因此，当断点附近的样本量不同时，就算是相同的断点设计可能也会得到不同的结果。

基于这些原因，当样本点数量足够多并且存在估计偏误的可能时，一个科学的对策是选择更小的带宽。此外，原文还提出了利用 “真实置信区间(Honest Confidence Interval)” 推断的方法，以获得良好的估计效果。

5. 真实置信区间

如果 $μ (x)$ 的条件期望函数在临界值两端能够随意变动的话，任何一条拟合的直线都能给出一个 $τ$ 的预测值，此时难以对 ATE 进行有效且正确的推断。因此，需要对 $μ (x)$ 作出进一步限制。假设 $μ \in M$ ，其中 $M$ 是一个函数集，并且构造出满足如下条件的置信区间：

\underset{N \to \infty}{lim inf} inf_{μ \in M} P_{μ} (τ \in C^{1 - α}) \geq 1 - α

其中， $P_{μ}$ 表明置信区间的覆盖效果很大程度上受期望函数形式的影响，将其称作真实地反映了 $M$ 的形式，一般来说，置信区间最好能够反映清晰易解释、有意义的函数形式。RDHonest是 R 中的一个包，能够构造出真实置信区间，并通过加入两条约束条件证明 $μ (x)$ 能用 $p$ 次多项式进行良好的拟合。

首先对 $μ (x)$ 的二阶导形式进行限制，也就是对函数形式的平滑程度进行限制。假设临界值两端的 $μ (x)$ 二阶可微，且两端期望函数二阶导的差值为常数 $K$ 。通过假设 $K$ 趋近于 0，可以将函数形式限定较为光滑的线性函数，而大的 $K$ 值则代表着更不光滑的函数形式。只需要假设 $μ (x)$ 处处二阶可导，得到如下定义：

M_{H} (K) = {μ : | μ^{'} (a) - μ^{'} (b) | \leq K | a - b | for all a, b \in R_{-} and all a, b \in R_{+}}

基于连续变量的真实置信区间更好构造，用 ${\tilde{τ}}_{h} = E (\hat{τ} ∣ X_{1}, \dots, X_{N_{h}})$ 指代处置效应估计系数的条件期望， ${\hat{σ}}_{N N}^{2} / N_{h}$ 指代 $\hat{τ}$ 的条件方差 $V (\hat{τ} ∣ X_{1}, \dots, X_{N_{h}})$ 的最近邻估计，则估计系数可以表达为如下形式：

{\hat{σ}}_{N N}^{2} = e_{1}^{'} {\hat{Q}}^{- 1} {\hat{Ω}}_{N N} {\hat{Q}}^{- 1} e_{1}, {\hat{Ω}}_{N N} = \frac{1}{N_{h}} \sum_{g = 1}^{G_{h}} n_{g} {\hat{σ}}_{g}^{2} m (x_{g}) m {(x_{g})}^{'}

其中， $n_{g}$ 是满足 $X_{i} = x_{g}$ 的观察值的个数， ${\hat{σ}}_{g}^{2} = \sum_{i : X_{i} = x_{g}} {(Y_{i} - {\bar{Y}}_{g})}^{2} / (n_{g} - 1)$ 是条件方差 $σ_{g}^{2} = V (Y_{i} ∣ X_{i} = x_{g})$ 的无偏估计，且 ${\bar{Y}}_{g} = n_{g}^{- 1} \sum_{i : X_{i} = x_{g}} Y_{i}$ 。由此， $t$ 统计量可以由下式得出：

\frac{\hat{τ} - τ}{{\hat{σ}}_{N N} / \sqrt{N_{h}}} = \frac{\hat{τ} - {\tilde{τ}}_{h}}{{\hat{σ}}_{N N} / \sqrt{N_{h}}} + \frac{{\tilde{τ}}_{h} - τ}{{\hat{σ}}_{N N} / \sqrt{N_{h}}}

中心极限定理使得上式第一部分在大样本情形下服从标准正态分布，而第二项拥有一个上界值 $r_{s u p} = \frac{sup_{μ \in M_{H} (K)} | {\tilde{τ}}_{h} - τ |}{{\hat{σ}}_{N N} / \sqrt{N_{h}}}$ 。由于上界值是由 $μ (x)$ 决定的，当 $x \geq 0$ 时为 $- K x^{2}$ ，当 $x < 0$ 时为 $K x^{2}$ ，可以将第二项的上界值表示为如下形式：

r_{sup} = - \frac{K}{2} \frac{\sum_{i = 1}^{N_{h}} w (X_{i}) X_{i}^{2} sign (X_{i})}{{\hat{σ}}_{N N} / \sqrt{N_{h}}}, w (X_{i}) = \frac{1}{N_{h}} \cdot e_{1}^{'} {\hat{Q}}^{- 1} M_{i}

这种对二阶导设定上界值方式可以构造出如下置信区间：

Proposition 1

令 $C V_{1 - α} (r)$ 指代 $∣ N (r, 1) ∣$ 分布的 $1 - α$ 分位数，那么所构造的置信区间，即下式，是反映 $M_{H} (K)$ 的真实置信区间。这一置信区间适用于离散型和连续型两种情形，并且由于考虑进了有限样本情形下估计偏误问题，因此能够适用于不同的带宽。

C_{B S D}^{1 - α} = (\hat{τ} \pm {c v}_{1 - α} (r_{sup}) \times {\hat{σ}}_{N N} / \sqrt{N_{h}})

第二是对临界值的设定误差进行限制，即所选择的临界值在给定带宽内能够得出拟合效果最好的模型，并且临界值的左极限 $lim_{x ↑ 0} δ (x)$ 不能大于所有临界值以下的变量的预测误差，临界值的右极限同理。

Proposition 2

利用