第 24 章 工具变量方法的应用：模糊断点回归
第 20 章的断点回归，和第 21–23 章的 IV，都是自然实验的重要例子。研究设计不如第二部分的随机化实验那么理想，可它们带着与随机化实验相近的特征。所以才叫自然实验。
把断点回归与 IV 叠在一起，就得到模糊断点回归（fuzzy regression discontinuity），又一种重要的自然实验。我先举例子，再给数学表述。
24.1 动机例子
第 20 章介绍的是断点回归。下面三个例子略有不同：接受的处理并不是驱动变量的确定性函数。驱动变量做的事，是在阈值处把接受处理的概率一下子改掉。
例 24.1 2000 年，印度政府启动总理乡村道路计划（Prime Minister’s Village Road Program）；到 2015 年，这项计划已资助将近 20 万个村庄修了全天候道路。Asher and Novosad (2020) 用村庄一级的数据，以断点回归去估新建支线公路对若干经济变量的效应。国家指南按 2001 年人口普查划出阈值，优先照顾更大的村庄，阈值本身带点任意。处理变量等于 1，表示结果被测量的那一年之前，村庄已经接到了新路。村庄人口与阈值之差并不决定处理，可它在阈值 \(0\) 处不连续地改了处理的概率。
例 24.2 Li et al. (2015) 用意大利两所大学 2004 至 2006 年入学的一年级学生数据，评估大学助学金对辍学率的因果效应。标准化家庭收入低于 15,000 欧元的学生有资格申请这笔助学金。为简单起见，我们把驱动变量定义成 15,000 减去标准化家庭收入。要真正拿到助学金，学生必须先申请。因此，资格与申请状态一起决定最终的处理状态。驱动变量单独并不决定处理，可它在阈值 \(0\) 处改了处理概率。
例 24.3 Amarante et al. (2016) 估计孕期暴露于一项社会救助计划，对儿童出生结局的影响。他们用的是乌拉圭国家社会紧急关注计划（Plan de Atención Nacional a la Emergencia Social）带来的断点。这是一项临时社会救助，瞄准最穷的 10% 家庭，实施于 2005 年 4 月至 2007 年 12 月。预测低收入得分低于事先定好的阈值的家庭被分配进计划。预测收入得分并不决定母亲是否在怀孕期间至少收到一次计划转账，可它改了最终接受处理的概率。出生结局包括出生体重、孕周等。
相对第 20 章那种（清晰）断点回归（sharp regression discontinuity），上面这些例子叫做模糊断点回归。例 24.1 与 24.2 的数据，我放在第 24.3 节来分析。
24.2 数学表述
令 \(X_i\) 为驱动变量，它决定
\[ Z_i=I(X_i\ge x_0), \]
其中 \(x_0\) 是阈值。接受的处理 \(D_i\) 未必等于 \(Z_i\)，可 \(\operatorname{pr}(D_i=1\mid X_i=x)\) 在 \(x_0\) 处有一个跳跃。图 24.1 比较清晰断点回归与模糊断点回归里接受处理的概率。它画的是模糊断点回归的一个特例：\(\operatorname{pr}(D=1\mid X<x_0)=0\)，与例 24.2 相符。
清晰： 模糊（左端为 0 的特例）：
pr(D=1 | X=x) pr(D=1 | X=x)
1 | ________ 1 | ........
| | | /
| | | ___/
0 |_____| 0 |_____|
x0 X x0 X原书图 24.1：清晰断点回归（左）与模糊断点回归（右）的处理分配。左边是从 \(0\) 到 \(1\) 的台阶；右边在 \(x_0\) 左侧为 \(0\)，阈值处跳起来，再往右慢慢升，虚线标着 \(1\)。
令 \(Y_i\) 为关心的结果。把 \(Z_i\) 看成分配的处理，就可以定义潜在结果 \(\{D_i(1),D_i(0),Y_i(1),Y_i(0)\}\)。\(Z\) 上的清晰断点回归，允许我们按定理 20.2 识别
\[ \begin{aligned} \tau_D(x_0) &= \mathrm{E}\{D(1)-D(0)\mid X=x_0\} \\ &= \lim_{\varepsilon\to 0^{+}}\mathrm{E}(D\mid Z=1,X=x_0+\varepsilon) - \lim_{\varepsilon\to 0^{+}}\mathrm{E}(D\mid Z=0,X=x_0-\varepsilon) \end{aligned} \]
以及
\[ \begin{aligned} \tau_Y(x_0) &= \mathrm{E}\{Y(1)-Y(0)\mid X=x_0\} \\ &= \lim_{\varepsilon\to 0^{+}}\mathrm{E}(Y\mid Z=1,X=x_0+\varepsilon) - \lim_{\varepsilon\to 0^{+}}\mathrm{E}(Y\mid Z=0,X=x_0-\varepsilon). \end{aligned} \]
把 \(Z\) 当作 \(D\) 的 IV，并在 \(X=x_0\) 处加上 IV 假定，再用定理 21.1，就可以识别局部依从者平均因果效应。
定理 24.1 假定处理由 \(Z=I(X\ge x_0)\) 决定，其中 \(x_0\) 是事先定好的阈值。假定在 \(x_0\) 的无穷小邻域里，单调性
\[ D_i(1)\ge D_i(0) \]
与排除限制
\[ D_i(1)=D_i(0) \implies Y_i(1)=Y_i(0) \]
成立。局部依从者平均因果效应（local complier average causal effect）定义为
\[ \tau_{\mathrm{c}}(x_0) = \mathrm{E}\{Y(1)-Y(0)\mid D(1)>D(0),X=x_0\}, \]
可由
\[ \tau_{\mathrm{c}}(x_0) = \frac{\mathrm{E}\{Y(1)-Y(0)\mid X=x_0\}}{\mathrm{E}\{D(1)-D(0)\mid X=x_0\}} \]
识别。再假定 \(\mathrm{E}\{D(1)\mid X=x\}\) 与 \(\mathrm{E}\{Y(1)\mid X=x\}\) 在 \(x=x_0\) 从右边连续，\(\mathrm{E}\{D(0)\mid X=x\}\) 与 \(\mathrm{E}\{Y(0)\mid X=x\}\) 在 \(x=x_0\) 从左边连续。若 \(\mathrm{E}(D\mid X=x)\) 在 \(x=x_0\) 处有非零跳跃，则局部依从者平均因果效应可由
\[ \tau_{\mathrm{c}}(x_0) = \frac{ \lim_{\varepsilon\to 0^{+}}\mathrm{E}(Y\mid Z=1,X=x_0+\varepsilon) - \lim_{\varepsilon\to 0^{+}}\mathrm{E}(Y\mid Z=0,X=x_0-\varepsilon) }{ \lim_{\varepsilon\to 0^{+}}\mathrm{E}(D\mid Z=1,X=x_0+\varepsilon) - \lim_{\varepsilon\to 0^{+}}\mathrm{E}(D\mid Z=0,X=x_0-\varepsilon) } \]
识别。
定理 24.1 是定理 20.2 与定理 21.1 的叠加。证明留给习题 24.1。
清晰与模糊断点回归里，关键都是指定阈值附近的邻域。实践里，邻域小则偏倚小、方差大，邻域大则偏倚大、方差小。也就是偏倚–方差权衡。也有一些按统计准则自动选邻域的手续，可它们依赖相当强的条件。更稳妥的做法，似乎是在一串邻域选择上做敏感性分析。
假定我们已经用带宽 \(h\) 指定了 \(x_0\) 的邻域。对 \(X_i\in[x_0-h,x_0+h]\) 的数据，可以用
\[ \hat\tau_D(x_0) = \text{\(D_i\) 对 \(\{1,Z_i,R_i,L_i\}\) 做 OLS 时 \(Z_i\) 的系数} \]
去估 \(\tau_D(x_0)\)，用
\[ \hat\tau_Y(x_0) = \text{\(Y_i\) 对 \(\{1,Z_i,R_i,L_i\}\) 做 OLS 时 \(Z_i\) 的系数} \]
去估 \(\tau_Y(x_0)\)，其中 \(R_i=\max(X_i-x_0,0)\)、\(L_i=\min(X_i-x_0,0)\)，与第 20 章相同。于是局部依从者平均因果效应可以估成
\[ \hat\tau_{\mathrm{c}}(x_0)=\hat\tau_Y(x_0)/\hat\tau_D(x_0). \]
这是一个间接最小二乘估计量。由定理 23.1，它与下面这件事数值上相同：
把 \(Y_i\) 对 \(\{1,D_i,R_i,L_i\}\) 做 TSLS，\(D_i\) 用 \(Z_i\) 当工具，得到的 \(D_i\) 的系数。
一句话：指定了 \(h\) 之后，估 \(\tau_{\mathrm{c}}(x_0)\) 就化成：用阈值附近的局部数据做一次 TSLS。
24.3 应用
24.3.1 重新分析 Asher and Novosad (2020) 的数据
回到例 24.1。我们可以按一串 \(h\)，基于结果 occupation_index_andrsn 计算点估计与标准误。
library("car")
road_dat = read.csv("indianroad.csv")
table(road_dat$t, road_dat$r2012)
road_dat$runv = road_dat$left + road_dat$right
## sensitivity analysis
seq.h = seq(10, 80, 1)
frd_sa = lapply(seq.h, function(h){
road_sub = subset(road_dat, abs(runv) <= h)
road_sub$r2012hat = lm(r2012 ~ t + left + right,
data = road_sub)$fitted.values
tslsreg = lm(occupation_index_andrsn ~ r2012hat + left + right,
data = road_sub)
res = with(road_sub, {
occupation_index_andrsn -
cbind(1, r2012, left, right) %*% coef(tslsreg)
})
tslsreg$residuals = as.vector(res)
c(coef(tslsreg)[2],
sqrt(hccm(tslsreg, type = "hc2")[2, 2]),
length(res))
})
frd_sa = do.call(rbind, frd_sa)原书图 24.2 给出结果。除非 \(h\) 很大，处理效应都不显著。
原书图 24.2：重新分析 Asher and Novosad (2020) 的数据，点估计与标准误来自 TSLS。上方面板：横轴是带宽 \(h\)，纵轴是估计；黑线是点估计，灰带大约是正负一个标准误，虚线在 \(0\)。点估计多在 \(-0.5\) 附近晃，\(h\) 较小时灰带盖住 \(0\)，\(h\) 大了才大体落到 \(0\) 下面。下方面板：样本量随 \(h\) 近似直线增加。
包 rdrobust 会自动选带宽。结果提示：接到新路，并没有显著改这个结果。
> library("rdrobust")
> frd_road = with(road_dat, {
+ rdrobust(y = occupation_index_andrsn,
+ x = runv,
+ c = 0,
+ fuzzy = r2012)
+ })
> res = cbind(frd_road$coef, frd_road$se)
> round(res, 3)| Coeff | Std. Err. | |
|---|---|---|
| Conventional | \(-0.253\) | \(0.301\) |
| Bias-Corrected | \(-0.283\) | \(0.301\) |
| Robust | \(-0.283\) | \(0.359\) |
24.3.2 重新分析 Li et al. (2015) 的数据
回到例 24.2。驱动变量是 15,000 减去标准化收入。分析时，我先把数据限制到驱动变量落在 \([-5000,5000]\) 的子集，再把驱动变量除以 \(5000\)，于是它夹在 \([-1,1]\) 里，阈值在零。我们可以按一串 \(h\) 计算点估计与标准误。
library("car")
italy = read.csv("italy.csv")
italy$left = pmin(italy$rv0, 0)
italy$right = pmax(italy$rv0, 0)
## sensitivity analysis
seq.h = seq(0.1, 1, 0.01)
frd_sa = lapply(seq.h, function(h){
italy_sub = subset(italy, abs(rv0) <= h)
italy_sub$Dhat = lm(D ~ Z + left + right,
data = italy_sub)$fitted.values
tslsreg = lm(outcome ~ Dhat + left + right,
data = italy_sub)
res = with(italy_sub, {
outcome -
cbind(1, D, left, right) %*% coef(tslsreg)
})
tslsreg$residuals = as.vector(res)
c(coef(tslsreg)[2],
sqrt(hccm(tslsreg, type = "hc2")[2, 2]),
length(res))
})
frd_sa = do.call(rbind, frd_sa)原书图 24.3 给出结果，提示大学助学金并没有显著改辍学率。
原书图 24.3：重新分析 Li et al. (2015) 的数据，点估计与标准误来自 TSLS。上方面板：点估计多在 \(-0.1\) 到 \(-0.2\) 之间，灰带在几乎所有 \(h\) 上都盖住 \(0\)。下方面板：样本量随 \(h\) 近似直线增加。
基于包 rdrobust 的结果，结论相同。
> library("rdrobust")
> frd_italy = with(italy, {
+ rdrobust(y = outcome,
+ x = rv0,
+ c = 0,
+ fuzzy = D)
+ })
> res = cbind(frd_italy$coef, frd_italy$se)
> round(res, 3)| Coeff | Std. Err. | |
|---|---|---|
| Conventional | \(-0.149\) | \(0.101\) |
| Bias-Corrected | \(-0.155\) | \(0.101\) |
| Robust | \(-0.155\) | \(0.121\) |
24.4 讨论
第 20 章与这一章，都把断点回归立在潜在结果对驱动变量的条件期望的连续性上。这个视角数学上更简单，可它识别的只是驱动变量阈值处那一点的局部效应。Hahn et al. (2001) 开了这条文献。
另一条不那么占主流的视角，立在局部随机化（local randomization）上（Cattaneo et al., 2015; Li et al., 2015）。若把驱动变量看成某个底层真相的带噪声度量，而阈值又带点任意，那么阈值附近的单元并不系统性地不同。这提示：在阈值的一个小邻域里，单元接受处理或对照，就像随机化实验那样带着随机。与第一种视角里选 \(h\) 的麻烦相似，这里也得决定：断点回归下的这场「随机化实验」该有多局部。把这套直觉写成数学并不容易；在第二种视角里，按一串 \(h\) 做敏感性分析，看起来同样合理。
更偏概念的讨论，见 Sekhon and Titiunik (2017)。
24.5 习题
24.1 定理 24.1 的证明
证明定理 24.1。
24.2 数据分析
第 24.3.1 节报告的是对 occupation_index_andrsn 的效应估计。另外四个结果变量是 transport_index_andrsn、firms_index_andrsn、consumption_index_andrsn 与 agriculture_index_andrsn，含义见原文。请估对这些结果的效应。
24.3 对 Li et al. (2015) 数据分析的一点反思
在 Li et al. (2015) 里，决定处理状态的一个关键变量是二值申请状态 \(A\)，它有潜在结果 \(A(1)\) 与 \(A(0)\)，分别对应处理 \(Z=1\) 与对照 \(Z=0\)。按定义，
\[ D(1)=A(1), \qquad D(0)=0, \]
于是依从者 \(\{D(1),D(0)\}=(1,0)\) 等价于 \(A(1)=1\)。所以
\[ \tau_{\mathrm{c}}(x_0) = \mathrm{E}\{Y(1)-Y(0)\mid A(1)=1,X=x_0\}. \]
第 24.3.2 节用整套数据去估 \(\tau_{\mathrm{c}}(x_0)\)。
另一种分析是只用 \(A=1\) 的那些单元。这时处理状态由 \(X\) 决定。可这种分析可能有问题，因为
\[\begin{align} &\lim_{\varepsilon\to 0^{+}}\mathrm{E}\{Y\mid A=1,X=x_0+\varepsilon\} - \lim_{\varepsilon\to 0^{+}}\mathrm{E}\{Y\mid A=1,X=x_0-\varepsilon\} \nonumber\\ &\qquad = \mathrm{E}\{Y(1)\mid A(1)=1,X=x_0\} - \mathrm{E}\{Y(0)\mid A(0)=1,X=x_0\}. \tag{24.1} \end{align}\]
证明 (24.1)，并解释这种分析为什么可能有问题。
注： (24.1) 左边是给定 \(A=1\) 时，\(X=x_0\) 处局部平均处理效应的识别公式。右边是两组单元上潜在结果均值之差：一组是 \((A(1)=1,X=x_0)\)，一组是 \((A(0)=1,X=x_0)\)。相关讨论见后面第 26 章。
24.4 推荐阅读
Imbens and Lemieux (2008) 基于潜在结果框架，给过断点回归的实践指南。Lee and Lemieux (2010) 综述过断点回归及其在经济学里的应用。