笔记

断点回归设计

断点回归分析

断点回归函数(Regression Discontinuity Design)

定义:\exists 一个处理变量DD,及其潜在结果Y0,Y1Y_0,Y_1。关于观测运行变量RR的间断函数,其表达式为:

Di=1{Ric}D_i = \mathbb{1}\{R_i \geq c\}

其中cc为断点。RR低于阈值cc时处理为0,否则为1。这意味着,当我们观测到RcR \geq c的条件时,我们观测到的是Y1Y_1,而我们观测到

R<cR < c的条件时,我们观测到的是Y0Y_0。也就是说虽然Y0(R)Y_0(R)Y1(R)Y_1(R)同时存在,但我们其实是没有办法同时观测的。这个cc就如一个阈值切换开关,我们只能观察到其中一个

1.00.80.60.40.20012345678910

Fig. 1. Assignment probabilities (SRD).

54Y132Y010012345678910

Fig. 2. Potential and observed outcome regression functions.

由于用的是/xymatrix库画的图表,因此很不美观,还请见谅

在图1中,我们能看到概率是从0直接到1,而没有其他取值,这种叫做精确断点回归设计

而图2展示的是潜在结果与观测到的结果回归函数。虚线表示在另一种处理状态下不可同时观测到的反事实趋势,实线表示断点两侧实际观测到的条件期望。断点处从Y0Y_0Y1Y_1的垂直跳跃,就是精确RDD在阈值处识别出来的局部处理效应。模糊断点回归设计不是由图2这种结果曲线定义的,而是由处理概率在断点处只发生部分跳跃定义的:断点右侧接受处理的概率高于左侧,但不一定从0跳到1。

精确断点回归设计

RDD的关键假设:

  • 连续性假设:潜在结果的条件期望在断点处连续
limrcE[YtiRi=r]=limrc+E[YtiRi=r]\lim_{r \rightarrow c^-} E[Y_{ti} | R_i = r] = \lim_{r \rightarrow c^+} E[Y_{ti} | R_i = r]
  • 局部随机化:在断点附近,个体无法精确操纵运行变量

  • 确定性的处理规则:在精确RDD中,处理状态由断点规则完全决定,即Di=1{Ric}D_i=\mathbb{1}\{R_i\ge c\}。如果处理概率只是在断点处增加、但没有从0跳到1,那就进入模糊RDD。

若以上条件成立,我们便可以估计阈值处的因果效应:

limrc+E[YtiRi=r]limrcE[YtiRi=r]=limrc+E[Y1iRi=r]limrcE[Y0iRi=r]=E[Y1iRi=r]E[Y0iRi=r]=E[Y1iY0iRi=r]\begin{aligned} \lim_{r \rightarrow c^+} E[Y_{ti} | R_i = r] - \lim_{r \rightarrow c^-} E[Y_{ti} | R_i = r] &= \lim_{r \rightarrow c^+} E[Y_{1i} | R_i = r] - \lim_{r \rightarrow c^-} E[Y_{0i} | R_i = r] \\ &= E[Y_{1i} | R_i = r]-E[Y_{0i} | R_i = r] \\ &= E[Y_{1i} - Y_{0i}| R_i = r] \end{aligned}

实际上,由于我们只估计了阈值处识别的效应,所以这估计的是断点处的局部平均处理效应

τSRD=E[Yi(1)Yi(0)Ri=c]\tau_{SRD}=E[Y_i(1)-Y_i(0)\mid R_i=c]

在这种条件下,RDD可以理解为一个局部的随机试验。对于正好处在阈值附近的人来说,是否接受处理像是随机的,因为有些人刚好低于阈值,有些人刚好高于阈值。虽然不像RCT那样的“金标准”,但已经非常接近了。

为了估计阈值处的处理效应,我们只需估计上述公式中的两个极限并进行比较。这时候又可以使用我们的老朋友:线性回归

我们设立一个线性方程,将“处理是否高于阈值”设置为一个虚拟变量来指示概率跳跃,然后与运行变量进行交互。为了让截距直接表示断点处的左右极限,通常先把运行变量中心化:

xi=Ric,Zi=1{Ric}x_i=R_i-c,\qquad Z_i=\mathbb{1}\{R_i\ge c\}

局部线性模型可以写成:

Yi=β0+β1xi+τZi+β2Zixi+εiY_i = \beta_0 + \beta_1 x_i + \tau Z_i + \beta_2 Z_i x_i+\varepsilon_i

这等同于在阈值之上和之下分别拟合一个线性回归。由于xi=0x_i=0正好对应断点cc,所以β0\beta_0代表断点左侧的极限,β0+τ\beta_0+\tau代表断点右侧的极限,τ\tau就是我们关心的跳跃。

更巧妙的有:

β0=limrcE[YiRi=r]β0+τ=limrc+E[YiRi=r]τ=limrc+E[YiRi=r]limrcE[YiRi=r]=E[Yi(1)Yi(0)Ri=c]\beta_0=\lim_{r \rightarrow c^-} E[Y_i | R_i = r] \\ \beta_0+\tau = \lim_{r \rightarrow c^+} E[Y_i | R_i = r] \\ \tau = \lim_{r \rightarrow c^+} E[Y_i | R_i = r] - \lim_{r \rightarrow c^-} E[Y_i | R_i = r] = E[Y_i(1)-Y_i(0)\mid R_i=c]

于是乎,我们就知道了我们关注的系数是τ\tau

那么我们在看断点前后的两个模型:断点前,斜率为β1\beta_1,而在断点后,斜率为β1+β2\beta_1+\beta_2。这样我们就得到了两个线性回归模型。

核加权

由于我们关注的其实是阈值前后的变化,这是我们所关注的,但,采用线性回归一个缺陷在于,它拟合其他数据点很好,但在阈值周围的点却拟合效果很差,那就得不到正确的结论了。那么,我们也就想到了,我们在阈值周围的点给予更高的权重,使得线性回归会更考虑阈值周围的点。

三角核函数

K(R,c,h)=1{Rch}(1Rch)K(R, c, h) = \mathbb{1}\{|R-c| \leq h\}\cdot \left(1-\frac{|R-c|}{h}\right)

我们详细看一下这个核函数:

第一部分:1{Rch}\mathbb{1}\{|R-c| \leq h\}: 为指示函数,这个hh为带宽,代表着我们接近阈值的程度。

第二部分:1Rch1-\frac{|R-c|}{h} 为加权函数,越远离阈值,即权重越小。带宽控制着权重减小的速度。

当我们画出函数图像的时候,是一个三角形状,其中最高点为我们关注的断点,也是权重最高的点,越远离两侧,权重越低。

于是乎,我们在模型输入里,就不能使用普通的最小二乘法。而是加权最小二乘法。

# OLS:
import statsmodels.formula.api as smf
model1 = smf.wls("y~x * threshold", data).fit() #等价于OLS

# WLS:
def kernel(R, c, h):
    indicator = (np.abs(R-c) <= h).astype(float)
    return indicator * (1 - np.abs(R-c)/h)
model2 = smf.wls("y~x * threshold", data,
                weights=kernel(data["x"], c=0, h=1)).fit()

# A*B = A + B + A : B

模糊断点回归设计

有些时候,并非所有个体都严格遵循断点规则,

π=limrc+P(Di=1Ri=r)limrcP(Di=1Ri=r),0<π<1\pi = \lim_{r \to c^+}P(D_i = 1 | R_i = r) - \lim_{r \to c^-}P(D_i = 1 | R_i = r), \qquad 0<\pi<1

其中 π\pi 是断点处处理概率的跳跃。

我们可以想象,越过阈值本应提高接受治疗的概率,但那些从不接受者并不会因此接受治疗;同样,低于阈值本应降低接受治疗的概率,但总是接受者仍然会接受治疗。在此情境下,我们定义潜在处理状态:Di(1)D_i(1)表示若个体位于阈值右侧时是否接受处理,Di(0)D_i(0)表示若个体位于阈值左侧时是否接受处理。

这时可以把阈值指示器Zi=1{Ric}Z_i=\mathbb{1}\{R_i\ge c\}看成一个局部工具变量。结果变量在断点处的跳跃是“简约式效应”(reduced-form effect),处理变量在断点处的跳跃是“第一阶段”(first stage)。模糊RDD用两者之比,把断点规则对结果的影响转换成对实际接受处理者的局部处理效应。

同样的,我们也假设潜在结果和潜在处理概率在断点处连续,并且需要一个处理状态的单调性假设:断点规则不能让某些个体反向改变处理状态。写成潜在处理状态就是

Di(1)Di(0),iD_i(1)\ge D_i(0),\qquad \forall i

也就是说,跨过阈值可以让一部分人从不处理变为处理,但不能让另一部分人从处理变为不处理。这里的单调性不是Y1i>Y0iY_{1i}>Y_{0i},因为结果是否更好是需要估计的问题,不能作为识别假设直接放进去。

一样的,我们依旧分两步:

第一阶段:

Di=γ0+γ1Zi+γ2(Ric)+γ3Zi(Ric)+νiD_i = \gamma_0 + \gamma_1 Z_i + \gamma_2 (R_i - c) + \gamma_3 Z_i \cdot (R_i - c) + \nu_i

其中,ZiZ_i为阈值指示器,即工具变量。

第二阶段:

Yi=α+τD^i+β1(Ric)+β2Zi(Ric)+ϵiY_i = \alpha + \tau \hat{D}_i + \beta_1 (R_i - c) + \beta_2 Z_i \cdot (R_i - c) + \epsilon_i

于是乎,我们可以得到局部平均处理效应的Wald估计量:

LATE^=limrc+E[YiRi=r]limrcE[YiRi=r]limrc+E[DiRi=r]limrcE[DiRi=r]=E[Yi(1)Yi(0)Di(1)>Di(0),Ri=c]\widehat{LATE} = \dfrac{\lim_{r \to c^+} E[Y_i|R_i=r] - \lim_{r \to c^-} E[Y_i|R_i=r]}{\lim_{r \to c^+} E[D_i|R_i=r] - \lim_{r \to c^-} E[D_i|R_i=r]} = E[Y_i(1) - Y_i(0) | D_i(1) > D_i(0), R_i=c]
def wald_rdd(data):
    weights = kernel(data["x"], c=0, h=1)
    denominator = smf.wls("treatment ~ x * threshold", data, weights=weights).fit()
    numerator = smf.wls("y ~ x * threshold", data, weights=weights).fit()
    return numerator.params["threshold"]/denominator.params["threshold"]

在现代RDD应用中使用局部多项式估计的异方差稳健标准误,并配合Calonico-Cattaneo-Titiunik提出的稳健偏差修正(robust bias correction, RBC)来构造置信区间。Bootstrap法常作为补充。

RDD检验

麦克雷利检验 (McCrary检验)

运行变量的密度检验

我们担心人们能够操纵自己在阈值附近的位置。比如治疗分配,需要病情恶化时才能接受的治疗,患者可以自报更严重的病情以获得新的治疗办法。在这些情况下,我们往往会观察到一种称为“聚集”的现象出现在运行变量的密度上。这意味着会有大量实体恰好位于阈值之上或之下。为了检验这一点,我们可以绘制运行变量的密度函数图,查看阈值附近是否存在尖峰

形式上检验:

H0: limrcfR(r)=limrc+fR(r)H_0:\ \lim_{r\to c^-}f_R(r)=\lim_{r\to c^+}f_R(r)

也可以写成对数密度跳跃:

θ=lnfR(c+)lnfR(c)\theta=\ln f_R(c^+)-\ln f_R(c^-)

若拒绝H0H_0,说明运行变量在断点处有明显不连续。这个结果不一定会直接推翻RDD,但需要解释排序机制,并且考虑甜甜圈RDD、限制样本或重新定义运行变量等的稳健性检验。

协变量平衡性检验

对预定的协变量使用RDD分析,检验是否存在显著性跳跃预定协变量不应该在断点处发生系统性跳跃。对每个协变量XkX_k重复RDD估计:

τ^Xk=limrc+E[XkiRi=r]limrcE[XkiRi=r]\hat{\tau}_{X_k}=\lim_{r\to c^+}E[X_{ki}\mid R_i=r]-\lim_{r\to c^-}E[X_{ki}\mid R_i=r]

如果多个政策实施前就已经确定的协变量在断点处显著跳跃,说明断点两侧个体可能不是局部可比的,连续性假设会受到怀疑。

这里应注意不是追求每一个协变量都显著,因为在多重检验下总会有偶然显著。更重要的是看跳跃是否系统性地出现,以及这些跳跃是否和结果变量有实质关系。预处理结果变量(pre-treatment outcome)尤其重要:如果实施政策发生前的结果变量已经在断点处跳跃,那么在实施政策后的跳跃就不一定能解释为处理效应。

参数选择

在RDD里,“参数选择”最核心的是选择断点附近多宽的局部窗口。这个窗口由带宽hh控制:

Rich|R_i-c|\le h

带宽太小,估计主要依赖离断点极近的样本,偏差小但方差大;带宽太大,样本更多,方差下降,但远离断点的函数形状会被带进来,偏差可能变大。因此RDD的参数选择本质上是在偏差和方差之间做权衡。

更一般地,令xi=Ricx_i=R_i-c,在断点两侧分别估计局部多项式:

minα+,β1+,,βp+i:xi0K(xih)(Yiα+j=1pβj+xij)2\min_{\alpha_+,\beta_{1+},\ldots,\beta_{p+}} \sum_{i:x_i\ge0}K\left(\frac{x_i}{h}\right) \left(Y_i-\alpha_+-\sum_{j=1}^p\beta_{j+}x_i^j\right)^2 minα,β1,,βpi:xi<0K(xih)(Yiαj=1pβjxij)2\min_{\alpha_-,\beta_{1-},\ldots,\beta_{p-}} \sum_{i:x_i<0}K\left(\frac{x_i}{h}\right) \left(Y_i-\alpha_- -\sum_{j=1}^p\beta_{j-}x_i^j\right)^2

MSE最优带宽:Imbens-Kalyanaraman

Imbens-Kalyanaraman带宽选择的目标是让点估计的**均方误差(MSE)**尽可能小。局部线性RDD的误差可以粗略分解为:

MSE(h)Bias(h)2+Var(h)MSE(h)\approx Bias(h)^2+Var(h)

其中局部线性估计的主导偏差与h2h^2同阶,方差与(nh)1(nh)^{-1}同阶:

Bias(h)2h4,Var(h)1nhBias(h)^2\propto h^4,\qquad Var(h)\propto \frac{1}{nh}

pp阶局部多项式,一般的MSE近似可以写成:

MSE{τ^p(h)}h2(p+1)Bp2+VpnhMSE\{\hat{\tau}_p(h)\} \approx h^{2(p+1)}B_p^2+\frac{V_p}{nh}

其中BpB_p对应主导偏差常数,VpV_p对应主导方差常数。最小化这个近似MSE,可以得到概念上的最优带宽:

hMSE=(Vp2(p+1)Bp2)12p+3n12p+3h_{MSE} = \left( \frac{V_p}{2(p+1)B_p^2} \right)^{\frac{1}{2p+3}} n^{-\frac{1}{2p+3}}

因此局部线性(p=1p=1)时,MSE最优带宽大致具有如下形式:

hMSE=CMSEn1/5h_{MSE} = C_{MSE}\cdot n^{-1/5}

其中CMSEC_{MSE}由断点附近的条件方差运行变量密度核函数常数以及断点两侧回归函数的二阶导数决定。更直观地说:

hMSE(断点附近的噪声大小断点附近的样本密度×曲率差异2)1/5n1/5h_{MSE} \propto \left( \frac{\text{断点附近的噪声大小}} {\text{断点附近的样本密度}\times \text{曲率差异}^2} \right)^{1/5} n^{-1/5}

如果断点附近结果变量噪声大,就需要更大的带宽来降低方差;如果断点两侧函数弯曲得很厉害,就需要更小的带宽来降低偏差。实现软件:R语言中的rdrobust/rdbwselect中的mserd

应该注意:MSE最优带宽主要用于点估计,但不一定能给出覆盖率最好的置信区间。因为当hh按MSE最优选择时,偏差往往还没有小到可以忽略,而置信区间可能低估不确定性。

稳健偏差修正:Calonico-Cattaneo-Titiunik

Calonico-Cattaneo-Titiunik(CCT)的核心贡献是:不要假装局部多项式估计没有偏差,而是显式估计并修正这个偏差,然后在标准误中把“估计偏差本身的不确定性”也纳入进去。

设主效应估计使用pp阶局部多项式和带宽hh,再用p+1p+1阶的局部多项式和另一个带宽bb估计主导偏差。可以把偏差修正后的估计量写成:

τ^bc=τ^p(h)Bias^p+1(h,b)\hat{\tau}_{bc} = \hat{\tau}_{p}(h)-\widehat{Bias}_{p+1}(h,b)

对应的稳健偏差校正(RBC, robust bias-corrected)置信区间为:

CIRBC=[τ^bc±z1α/2SE^RBC]CI_{RBC} = \left[ \hat{\tau}_{bc} \pm z_{1-\alpha/2}\cdot \widehat{SE}_{RBC} \right]

这里SE^RBC\widehat{SE}_{RBC}不仅包含原估计量的抽样误差,也包含偏差估计带来的额外不确定性。因此,RBC置信区间通常比“只对点估计做偏差修正、但仍用普通标准误”的做法更可靠。

CCT 的核心是RBC置信区间;在后续Calonico-Cattaneo-Farrell等工作和rdrobust实现中,覆盖误差率最优带宽(coverage error rate optimal bandwidth, CER-optimal bandwidth)被系统化。它的目标不是最小化点估计的MSE,而是改善置信区间的覆盖率。

交叉验证: Cross-validation

交叉验证的想法是:选择一个带宽hh,使得局部回归对观测结果的预测误差最小。形式上可以写成:

hCV=argminhi(Yim^i(Ri;h))2h_{CV} = \mathop{\arg\min}\limits_h \sum_i \left(Y_i-\hat{m}_{-i}(R_i;h)\right)^2

其中m^i\hat{m}_{-i}表示不使用第ii个观测值时得到的预测函数。

但在RDD中,交叉验证不是最常用的主方法。原因是RDD关心的是断点处的左右极限,而普通预测误差可能更重视整个窗口内的拟合效果。一个带宽如果让整体预测误差很小,并不代表它对断点处的跳跃估计最合适。因此,交叉验证更适合作为辅助参考。

其他参数选择

除了带宽,RDD还需要说明以下选择:

  1. 多项式阶数:主规格通常使用局部线性(p=1p=1)。局部二次(p=2p=2)可以作为稳健性检查,但一般而言,不使用高阶全局多项式,因为它们容易被远离断点的数据牵引,导致断点附近拟合不稳定。
  2. 核函数:三角核是常用主规格。矩形核或Epanechnikov核可作为稳健性检查。
  3. 左右带宽:可使用共同带宽(通常做法),也可以允许断点左右两侧有不同带宽。若两侧样本密度、方差或曲率明显不同,不同带宽可能更合理。
  4. 标准误:若数据存在聚类、重复观测或离散运行变量,应使用相应的聚类稳健标准误或离散运行变量调整。否则显著性可能过于乐观。

稳健性检验

带宽敏感性分析

主结果不应该只在一个特定带宽下成立。可以围绕主带宽h\*h^\*报告一组估计值:

h{0.5h\*,0.75h\*,h\*,1.25h\*,1.5h\*,2h\*}h\in\{0.5h^\*,\,0.75h^\*,\,h^\*,\,1.25h^\*,\,1.5h^\*,\,2h^\*\}

并分别计算

τ^(h)=α^+(h)α^(h)\hat{\tau}(h) = \hat{\alpha}_+(h)-\hat{\alpha}_-(h)

如果结论只在某一个非常特殊的带宽下显著,而稍微扩大或缩小带宽后就消失,那么RDD结果的可信度会下降。相反,如果估计值的方向和量级在合理带宽范围内都比较稳定,说明结果不依赖任意的窗口选择。

多项式阶数和核函数敏感性

主规格一般采用局部线性模型:

Yi=β0+β1xi+τZi+β2Zixi+εiY_i=\beta_0+\beta_1x_i+\tau Z_i+\beta_2Z_ix_i+\varepsilon_i

稳健性检验可以改用局部二次模型:

Yi=β0+β1xi+β2xi2+τZi+γ1Zixi+γ2Zixi2+εiY_i=\beta_0+\beta_1x_i+\beta_2x_i^2+\tau Z_i+\gamma_1Z_ix_i+\gamma_2Z_ix_i^2+\varepsilon_i

如果局部线性和局部二次给出非常不同的结果,需要检查是否存在带宽过大、函数形状过于复杂、或者断点附近样本不足的问题。

核函数也可以从三角核改为矩形核或Epanechnikov核。若更换核函数后结论完全改变,说明估计可能对权重设定比较敏感。

安慰剂断点检验

在真正断点cc以外选择一些“假断点”c0c_0,重复RDD估计:

τ^(c0)=limrc0+E[YiRi=r]limrc0E[YiRi=r],c0c\hat{\tau}(c_0) = \lim_{r\to c_0^+}E[Y_i\mid R_i=r] - \lim_{r\to c_0^-}E[Y_i\mid R_i=r], \qquad c_0\ne c

如果在很多非政策断点处也出现显著跳跃,说明结果可能来自函数形状、样本结构或其他制度变化,而不是断点规则本身。

安慰剂断点的选择通常是预先给出一组远离真实断点、且不会受到真实处理规则影响的位置,或者系统地在多个候选断点上画出估计结果。

甜甜圈RDD

甜甜圈RDD(donut RDD)会排除断点附近极近的一小段样本:

τ^donut(δ,h)=τ^(h; Ric>δ)\hat{\tau}_{donut}(\delta,h) = \hat{\tau}\left(h;\ |R_i-c|>\delta\right)

其中δ\delta是被挖掉的半径。

这个检验适合处理两类问题:第一,运行变量在断点附近可能被精确操纵,即McCrary检验有尖峰存在;第二,运行变量存在堆积、取整或测量误差,导致大量观测刚好卡在断点附近。如果排除这些极近样本后估计结果仍然相近,说明结论对断点处的异常堆积不敏感。若结果变化大,需要重新审视断点操纵、测量误差或排序问题。

模糊RDD的第一阶段稳健性

模糊RDD还必须检查处理概率是否真的在断点处跳跃:

π=limrc+E[DiRi=r]limrcE[DiRi=r]\pi = \lim_{r\to c^+}E[D_i\mid R_i=r] - \lim_{r\to c^-}E[D_i\mid R_i=r]

如果π\pi很小或不显著,那么局部Wald估计量

τ^FRD=τ^Yπ^\hat{\tau}_{FRD}=\frac{\hat{\tau}_Y}{\hat{\pi}}

会变得不稳定,类似弱工具变量问题。在使用模糊RDD时,应该注意报告结果变量跳跃处理变量跳跃、Wald比值和第一阶段图形。

稳健性检验总结

所有的稳健性检验只能说明:在我们能观察到的维度上,没有发现明显违背RDD假设的证据。