第20章断点回归

20.1 它是如何起作用的？

圣地亚哥是美国的一座大城市，地理覆盖范围相当广，超过 \(300\) 平方英里。它也相当富裕，截至 \(2019\) 年，平均家庭年收入超过 \(85000\) 美元，比美国全国平均水平高出约 \(50\%\)。当你从该地区向南行进时，景象会略有变化。城市的一些南部地区没那么富裕了。当你到达圣伊西德罗，这座城市紧邻墨西哥边境的行政区时，家庭年收入已经下降到更接近 \(50000\) 美元或 \(55000\) 美元。你往南走得越远，预期收入就越低。

但我们在这座城市的旅程，与我们跨越边境进入墨西哥蒂华纳时发生的情况相比，根本不算什么。我们从市中心开车 \(16\) 英里到圣伊西德罗，目睹收入下降 \(25\%\)。但你可以在一个家庭年收入为 \(50000\) 美元的地区下车，然后只需走几步越过边境进入蒂华纳。家庭年收入突然大幅下降到 \(20000\) 美元。

当然，圣地亚哥南部与市中心地区相比，可能在人口、机会或地理方面存在差异，这些差异可以解释收入上的部分不同。但仅仅是往南走，无法解释当我们跨出那一步越过边境时看到的收入暴跌情况。从市中心到圣伊西德罗，变化是持续的。但在边境处，变化是跳跃式的。有些东西从根本上就不同了。我们遇到了一个临界点。

每当某种处理是不连续分配的 —— 刚好在一条线一侧的人能得到它，而刚好在另一侧的人得不到它时，我们就可以相当确定，我们所看到的差异不只是因为他们往南走了多远而存在的不同。如果没有那条线，他们可能会有一点不同，但不会有那么大的不同。我们可以把这种巨大的变化归因于分隔墨西哥和美国的那条线

所以如果他们存在差异，我们大概可以将这种差异归因于在那个临界点上发生的无论是什么情况¹。

断点回归关注的是在临界点处分配的处理。刚好在一侧，没有处理；刚好在另一侧，有处理²。

在这里有必要明确一些术语：

运行变量\(（Running \space variables）\)：运行变量，也被称为强制变量，是决定你是否接受处理的变量。例如，如果一位医生测量你的血压，并且如果你的收缩压高于 \(135\)，就会给你开降压药，那么你的血压就是运行变量。
临界点\(（Cutoff）\)：临界点是运行变量的取值，它决定你是否接受处理。用同样的血压例子，在这种情况下，临界点就是收缩压为 \(135\)。如果你的收缩压高于 \(135\)，你就会得到药物；如果低于 \(135\)，就得不到³。
带宽\(（Bandwidth）\)：认为刚好处于临界点两侧的人，除了临界点的因素外，本质上是相同的，这是合理的。但离临界点更远的人（比如圣地亚哥市中心与墨西哥更深处的人）可能会因为临界点之外的原因而存在差异。带宽指的是你愿意将临界点周围多大范围的区域视为具有可比性。是美国 / 墨西哥边境两侧各十英尺？一千英尺？八十英里？

断点回归研究设计的核心在于：\(（a）\)解释运行变量通常如何影响结果；\(（b）\) 聚焦于临界点附近、带宽范围内的观察对象；\(（c）\)对比刚好接受处理的人群和刚好未接受处理的人群，以得出处理的效应。

世界上有很多这样的临界点。相当多的政策在设计时就明确考虑了临界点。收入刚好略高一点？你可能就没有资格参加某些经经济状况调查的项目。或者考试分数刚好略低一点？你可能就没有资格参加资优项目。

这些临界点可能是地理上的 —— 我们已经有了圣地亚哥 / 蒂华纳的例子。但同样，就住在时区边界的某一侧或另一侧附近？你可能要早起一个小时。或者，如果你就住在某个警察辖区边界的一侧附近，你可能会经历不同的警务政策。

政治是另一个存在硬性临界点的领域。在选举中，在选举中排名前两位的候选人里，你获得了 \(49.999\%\) 的选票？你得不到那个职位。但如果获得 \(50.001\%\) 的选票呢？好好享受你的职位吧！现在我们就能看到你的当选效应了。

在这些案例中的每一个，我们都能相当合理地想象，刚好处于临界点两侧附近的情况是具有可比性的，而它们之间的任何差异实际上都是处理因素导致的。换句话说，我们已经识别出了处理的效应。

这种方法的目标是仔细挑选变异情况，让我们可以忽略各种后门干扰因素，而实际上无需对它们进行控制。

我们什么时候可以应用断点回归？我们要找的是某种基于临界点分配的处理。存在决定处理的运行变量，并且存在某个临界值。如果你刚好在临界值的一侧，你不会得到处理；如果你在另一侧，你会得到处理⁴。另外，由于我们的策略是假设接近临界点的人实际上是被随机分配的⁵，不应该有任何明显的因素阻碍这种随机性 —— 人们不应该能够操纵运行变量来选择自己的处理，而且那些决定临界值的人也不应该能够根据发现谁具有哪些运行变量值来做出设定临界值的选择。

这给了我们一个如 \(图20.1\) 所示的因果图。在这里，我们有一个运行变量\(（Running \space Variables）\)，它有一条通过 \(Z\) 的后门路径（或许我们无法对其进行控制）。运行变量对结果\(（Outcome）\)也有直接影响。

现在，我们实际上并不关心运行变量的影响。我们关心的是处理变量\(（Treatment）\)的影响。处理变量有一条通过运行变量的后门路径。然而，如果我们控制了运行变量中除了是否高于临界值之外的所有变异，那么我们就可以关闭这条后门路径，识别出处理对结果的影响⁶。

基于这个图表，断点回归似乎有点像是控制一个变量（运行变量）然后声称我们已经识别出了效应。但实际上它比这更有意思。我们实际上是在隔离一条前门路径，就像第\(9\)章（寻找前门）或第\(19\)章（工具变量）中那样，而不是关闭后门路径。通过只观察临界值附近的情况，我们消除了所有不位于” \(AboveCutoff \rightarrow Treatment \rightarrow Outcome\) “这条路径上的变异。一旦我们把范围限定在那个临界值附近，其他东西实际上都不应该有变化了，所以所有其他变量的影响都消失了！

所以，当然，在\(图20.1\)中，我们可以通过对运行变量进行常规的统计调整来识别效应。但如果存在其他后门路径呢？比如，甚至存在处理本身的后门路径，例如如果图中有 \(Z \rightarrow Treatment\) 的路径？控制运行变量并不能解决这个问题，但断点回归可以，它通过隔离”\(AboveCutoff \rightarrow Treatment \rightarrow Outcome\)“这条路径，让我们可以忽略任何其他指向处理变量的箭头⁷。

那么，断点回归到底是如何做到这些超棒的事情的呢？我们只需做出几个基本选择，就能进行断点回归分析。一旦做出这些选择，过程就相当简单了。

选择一种方法来预测临界值两侧的结果。
选择一个带宽（可选）。

就这么简单⁸！让我们把这些选择用在一些模拟数据上。你可以在\(图20.2\)中看到整个过程是如何展开的。

在部分\(（a）\)中，我们看到了一个典型的断点回归设置：\(X\) 轴上是一个运行变量。即便在我们不在临界值附近时，运行变量似乎也与结果相关，这是没问题的。这里有一个明显的向上的斜率。我们还有一个施加了处理的临界值。并且我们看到在那个临界值处结果出现了一个跳跃！这正是我们想要看到的。

接着看部分\(（b）\)，我们需要一种方法来预测刚好在临界值两侧的结果会是怎样的。我们可以使用很多不同的预测方法 —— 回归、局部均值等等⁹。在这种情况下，我们使用带有\(LOESS\)的局部均值，它只是一条基于回归预测的平滑线，而该回归只使用每个 \(X\) 轴值一定范围内的点。

在部分\(（c）\)中，我们可以看到我们限制了带宽（这会修改我们在\(（b）\)中使用的方法所做出的预测）。我们真正只对临界值附近的值感兴趣。远离临界值的值仍然有价值 —— 我们可以用它们来找出趋势和线条，让我们能更好地预测临界值处的值。但如果离得太远，你就会引入大量我们原本通过断点回归想要避免的后门变异。不过，也许我们计划通过调整结果与运行变量之间的关系以及回归来解决这个问题，那样我们就可以跳过带宽这一步。我们之后会讲到这方面的细节¹⁰，但现在我们选择使用带宽。在这里，我们只使用非常接近临界值的数据来计算我们对临界值附近的预测。其他所有数据都被忽略。

最后，基于我们的带宽和预测方法，我们得到了刚好在临界值两侧的结果预测值。然后，我们的断点回归估计值就是预测值在临界值处的跳跃幅度。我们可以在部分\(（d）\)中看到这一点，其中估计值是从刚好在左侧的预测值到刚好在右侧的预测值的垂直线。

到目前为止，我展示的例子中，临界值是处理的唯一决定因素。但通常情况并非如此。在很多断点回归的应用中，处于临界值的这一侧或那一侧只会改变接受处理的概率。在这些情况下，我们会有所谓的模糊断点回归，而不是处理率从\(0\%\) 跃升至\(100\%\)的精确断点回归。

模糊断点回归：一种断点回归设计，其中处理并非完全由处于断点某一侧所决定。

\(图20.3\) 展示了在精确\(（sharp）\)和模糊\(（fuzzy）\)设计中，处理概率可能会如何变化¹¹。虽然两者都显示在临界值处处理概率出现了跳跃，但精确设计表明，除了临界值处，其他任何地方都完全没有变化。

模糊断点回归被大量应用的一个主题例子是研究退休的影响。在许多职业和国家，存在特定的年龄，或者你在工作岗位上的年数，在这些节点会有新政策生效 —— 养老金收入、退休基金的获取途径等 —— 并且退休率在这些节点会显著跃升。不过，并不是每个人一有资格就退休，有些人还会提前退休，这就使得这成为一种模糊的不连续性（模糊断点回归的情形）。

\(Battistin \space et \space al.（2009）\)的论文是众多将模糊断点回归应用于退休研究的论文之一。在他们的论文中，他们研究了消费情况以及消费在退休时是如何变化的。具体来说，他们想知道退休是否会导致消费立即下降¹²。如果确实下降了，又是为什么？¹³

Battistin, Erich, Agar Brugiavini, Enrico Rettore, and Guglielmo Weber. 2009. “The Retirement Consumption Puzzle: Evidence from a Regression Discontinuity Approach.” American Economic Review 99 (5): 2209–26.

他们研究了\(1993-2004\)年意大利的数据，在这些数据中，他们掌握了人们何时有资格领取养老金的信息。在刚好达到领取资格的年龄时，男性退休的比例有超过 \(30\) 个百分点的跃升。这不是从 \(0\%\) 到 \(100\%\) 的跃升（因此属于模糊断点的情况），但 \(30\) 个百分点的跃升是我们可以利用的！我们可以在\(图20.4\)中看到他们所研究的其中一年里的这种跃升情况。很明显，在那个临界值处确实发生了一些变化。

他们想要研究在那个临界值处消费是如何变化的，就像你通常用非模糊的断点回归所做的那样。但如果他们只是把它当作普通情况来处理，他们的结果会大错特错。他们会表现得好像临界值让所有人都退休了，而实际上它只使退休率提高了\(30\)个百分点。退休率\(30\)个百分点的变化对消费的改变，会比\(100\)个百分点的变化带来的改变要小，所以我们估计的效应自然会令人失望！但这种失望只是因为我们的估计预期会有\(100\)个百分点的变化。

作者们想要恰当地调整估计量的预期。所以他们用临界值的效应来进行缩放。大致来说，他们通过将观察到的效应除以\(30\)个百分点的跃升来缩放这些效应，以了解如果每个人都接受了处理（都退休了），变化会有多大。

他们发现了什么呢？他们确实发现，在人们退休的那一刻，消费出现了下降。但很方便的是，他们的断点回归设计不仅识别出了退休对整体消费的影响，还识别出了退休对几乎所有方面的影响，在这种情况下，就是对特定类型消费以及家中孩子数量的影响。他们发现，退休导致的消费下降，很大程度上是由于诸如利用额外的闲暇时间做饭而不是从餐馆点餐、不再需要漂亮的工作服装，以及成年子女往往会搬出家这些原因。这并不像消费下降看起来的那么糟糕。

肯定有什么陷阱。嗯，某种程度上是这样的。和任何事情一样，我们需要问问自己，要让这一切都成立需要做出什么假设，利弊是什么，以及我们最终到底估计的是什么。

首先，和任何其他我们只隔离能识别效应的那部分变异的方法一样，我们需要假设在那部分变异中没有可疑的情况发生。我们当然是在估计处理效应，但实际上我们是在估计临界值的效应，并将其（或者如果是模糊断点的话，将其部分）归因于处理。所以，如果在临界值处还有其他任何东西发生了变化，我们就会遇到麻烦。我们需要做出的假设是，结果在临界值处是平滑的。也就是说，如果处理状态在临界值处实际上没有变化（如果临界值附近没有人被处理，或者所有人都被处理，或者所有人被处理的概率都相同），那么就不会有什么跳跃或不连续性可言¹⁴。

假设我们想了解一个资优学校项目的效应，该项目要求数学测试成绩达到\(100\)分中的\(95\)分才能入学。数学测试得到\(96\)分或更高，是既能让你进入资优项目，又能让你进入科学博物馆实地考察项目吗？如果是这样的话，你就无法用断点回归只得到资优项目的效应。没办法把它和那些实地考察项目的效应区分开。

即使在临界值处没有另一个明确的处理发生，这也是个问题。也许资优项目是唯一一个以\(95\)分成绩为临界值分配的项目，但如果碰巧得\(94\)分的人和得\(96\)分的人里高个子的比例差异极大。即使得\(94\)分和得\(96\)分是完全随机的，身高差异也是偶然的，也很难将资优项目的效应与身高的效应区分开来。

其次，这件事的关键在于，如果我们放大到足够近的范围，人们基本上是被随机分配到临界值两侧的¹⁵。，¹⁶所以我们需要假设，运行变量在临界值附近是随机分配的，并且没有被操纵。如果老师们想让更多学生进入资优项目，于是把得\(94\)分的孩子的试卷偷偷重新评分，让他们变成得\(96\)分，那么突然之间，你得\(94\)分还是 \(96\)分就不是真的随机了。我们就不能再说\(94\)分和\(96\)分的人之间的比较能识别出效应了！

此外，由于这种方法基于能够足够放大样本，使得在断点两侧的分配是随机的，我们需要一个测量足够精确的运行变量，以便能够放大到必要的水平。当然，或许\(94\)分和\(96\)分的考生是可比较的。但如果测试数据仅以\(5\)分为区间报告，因此我们必须将\(90-94\)分区间的考生与\(95-100\)分区间的考生进行比较呢？在这种情况下，我们可能对随机分配假设的信任度会大打折扣。

除了关键假设外，我们的第三个担忧是关于统计功效的。断点回归设计必然聚焦于数据中一个极小的片段——即恰好位于断点附近的个体。实际上，\(94\)分、\(95\)分或\(96\)分的考生数量并不多。或者，从\(图20.2（c）\)中可以看出我们有多少数据被灰化并丢弃。因此，仅比较断点两侧的个体时，我们严重限制了样本量，从而使估计结果更加嘈杂。我们可以通过纳入远离断点的更多数据来解决这一问题，但这会重新引入我们试图消除的偏差¹⁷。

最后，如果我们完成了所有这些步骤，我们估计的是哪种类型的处理效应（第\(10\)章）？在这一点上，断点回归设计相对容易理解。我们仅利用断点附近的变化。因此，我们得到的是恰好位于断点附近个体的处理效应。这是一种局部平均处理效应，即针对那些恰好处于处理分配边际上的个体的加权平均处理效应。

在断点回归设计中，获取局部平均处理效应有其利弊。利在于，如果我们利用研究结果将处理扩展到更多个体，我们可能会通过稍微移动断点以涵盖更多人来实现。新接受处理的个体与我们刚刚估计的断点处个体非常相似。毕竟，了解天才项目对一个数学考试只得\(5\)分（满分\(100\)）的人的影响并不太有用——他们不太可能被纳入该项目。了解局部效应可能比了解整体平均处理效应更有实际意义。

当然，在其他情况下，我们可能更希望获得平均处理效应。在许多案例中，我们研究的是本可适用于所有人的处理效果，只是恰巧能够通过基于断点的分配方式识别其效应。例如，假设一项救济计划向大部分民众发放\(1600\)美元的救济支票，并以\(75000\)美元作为收入门槛。如果我们使用断点回归分析这些支票的效果，我们只能得到收入非常接近\(75000\)美元人群的效应。但我们永远无法知道年收入\(10\)万美元人群的效应如何。更令人担忧的是，我们无法了解年收入接近\(2\)万美元人群的效应——而这些人才是最需要帮助的群体。尽管年收入\(2\)万美元的人群确实收到了支票，但由于他们远离断点，我们无法通过断点回归来分析支票对他们的影响。

20.2 它是怎么工作的

20.2.1 使用普通最小二乘法的断点回归

至此，我已从”断点”角度阐述了断点回归，但尚未涉及”回归”部分。这似乎不妥，因此让我们完善它。预测断点两侧结果的方法有很多，但我们首先从一个基于交互项的简单方法开始：

\[Y = \beta_0 + \beta_1 (Running - Cutoff) + \beta_2 Treated + \beta_3 (Running - Cutoff) \times Treated + \varepsilon \quad (20.1)\]

这是断点回归的简单线性方法，其中\(Running\)是运行变量，我们通过使用\(（Running-Cutoff）\)将其以断点为中心进行标准化。\(（Running-Cutoff）\)在断点左侧取负值，在断点处为零，在右侧取正值。这里讨论的是清晰断点回归，因此 \(Treated\) 既是接受处理的指示变量，也是高于断点的指示变量——这两者是相同的¹⁸，您完全可以改用 \(AboveCutoff\) 来代替 \(Treated\)。该模型通常使用异方差稳健标准误进行估计，因为在我们拟合的直线中，断点及整体形状在大多数情况下都可能存在异方差性¹⁹。我们暂时忽略带宽问题，后续再讨论——无论使用全部数据还是限制在断点附近的带宽范围内，这种回归方法都适用。

注意\(方程（20.1）\)中未包含控制变量。在大多数其他章节中，我这样做是为了帮助您专注于研究设计。但此处是刻意为之。断点回归的核心思想是在断点两侧存在近乎随机的分配。您不需要控制变量，因为设计本身应该关闭所有后门。没有开放的后门？就不需要控制变量。添加控制变量意味着您不相信断点回归方法所需的假设，这会使整个分析在读者眼中显得可疑²⁰。

添加控制变量并非错误，但如果您认为处理完全由断点分配，则应仔细考虑是否真的需要它们——在这种情况下，您要关闭的是哪条后门？您必须认为在断点附近存在某些后门（这并非不可能，但并非既定事实）。在模糊断点回归中，控制变量可能更为必要，因为除了断点之外，肯定还存在其他决定处理的因素，因此可能存在需要担忧的后门。也许您有理由相信，只有在控制某些变量的条件下才存在随机分配，特别是当为了增加有效样本量而不得不接受更多远离断点的观测值时。一如既往，您仍然可以添加控制变量以尝试减小残差规模，从而缩小标准误。

我们最终得到两条直线——一条在断点未处理侧，截距为 \(\beta_0\)，斜率为 \(\beta_1\)；另一条在断点处理侧，截距为 \(\beta_0 + \beta_2\)，斜率为 \(\beta_1 + \beta_3\)。将此回归应用于某些模拟数据，我们得到\(图20.5\)。

啊，但别忘了——由于我们将运行变量以断点为中心调整为\(0\)，这两条直线的截距就是我们在断点处的预测值²¹。因此，我们对断点处两条直线跳跃幅度的估计就是截距的变化量，即 \((\beta_0 + \beta_2) - (\beta_0) = \beta_2\)。\(\beta_2\) 给出了我们的断点回归效应估计值。

这种方法的优缺点是什么？一个明显的优点是操作简单，但这可能不是选择估计方法的好理由。

真正的优缺点其实是相互关联的，这与我们使用拟合形状来预测断点处结果的事实有关。这样做有一些明显的好处，尤其是在统计精度方面。使用形状让我们能够整合远离断点的信息，以揭示结果随着越来越接近断点而呈现的变化轨迹。看到 \(图20.5\) 中断点右侧明显的正向趋势了吗？在趋势变得明显之前，数据中存在一段平坦甚至略微下降的区域——如果不整合远离断点的数据，就很难清晰地捕捉到该趋势，这可能导致我们高估来自右侧的断点处的值，如\(图20.6\)所示。

然而，应用拟合形状仅在我们使用正确的拟合形状时才有效。如果我们选择了错误的形状（如\(图20.5\)中断点左侧明显所示），从形状得到的预测将是错误的。

当然，当我们使用普通最小二乘法和拟合形状时，这总是个问题。但这里问题尤其严重，因为拟合形状在可用数据边缘往往误差最大。当达到数据的最边缘时，预测可能突然转向奇怪的方向，因为您选择的形状强制其继续沿该方向延伸。我们再次在 \(图20.6\) 左侧看到这一点，比较虚线（突然转向）与显示断点附近实际斜率的实线。

一个明显的直觉是尝试更灵活的形状。我们可以采用与\(方程（20.1）\)完全相同的思路，但对
\[(Running - Cutoff)\]使用更灵活的函数：

\[Y = \beta_0 + f(Running - Cutoff, Treated) + \beta_2Treated + \epsilon \quad (20.2)\]

其中 \(f()\) 是 \(Running - Cutoff\) 的某个非线性函数，针对 \(Treated = 0\) 和 \(Treated = 1\) 分别有不同版本。例如，这可以是二阶多项式：

\[Y = \beta_0 + \beta_1(Running - Cutoff) + \beta_2(Running - Cutoff)^2 + \beta_3Treated + \beta_4(Running - Cutoff) \times Treated + \beta_5(Running - Cutoff)^2 \times Treated + \epsilon \quad (20.3)\]

\(方程（20.3）\)再次产生两条线——断点左侧的抛物线和右侧的抛物线。与之前类似，\(Treated\) 的系数（此处为 \(\beta_3\)）显示了预测在断点处的跳跃幅度，给出了断点回归估计值。

然而，我们需要对更灵活的\(OLS\)线保持谨慎。通过使形状更加灵活来修正\(RDD\)中的”错误形状”问题是很困难的，因为更灵活的形状在数据边缘附近会做出比简单形状更奇怪的预测和急剧转向无穷大的行为。\(图20.7\)显示了使用\(（a）\)方程\(（20.3）\)中的\(2\)阶多项式、\(（b）\)\(4\)阶多项式以及 \((c)\) \(25\)阶多项式对相同数据估计的断点回归结果。添加所有这些多项式后，我们得到了越来越奇怪的结果，在断点附近出现巨大波动，甚至似乎不能很好地跟踪数据。

由于数据是模拟的，我可以告诉您真实效应为\(0.15\)，真实模型为\(2\)阶多项式。线性模型给出的效应为\(0.29\)。\(2\)阶、\(6\)阶和\(25\)阶多项式模型分别给出\(0.17\)、\(0.26\)和…即使\(25\)阶多项式碰巧接近真实效应，也仅仅是运气使然_。看看那条曲线在断点附近多么疯狂！灵活性并不总是好事。

由于这个问题，在使用普通最小二乘法进行断点回归时，不建议使用超过\(2\)阶的项\(（Gelman \space \& \space lmbens, 2019）\)。如果需要拟合复杂形状，可以尝试通过带宽限制数据范围并使用更简单的形状。当您通过带宽放大观察时，所有这些扭曲和曲线往往会消失，一条直线可能就足够了。这可以在\(图20.7（d）\)中看到，其给出的效应为\(0.19\)，而真实效应为\(0.15\)。它的表现与\(2\)阶多项式差不多，但不需要知道真实底层模型是\(2\)阶的。

Gelman, Andrew, and Guido W. Imbens. 2019. “Why High-Order Polynomials Should Not Be Used in Regression Discontinuity Designs.” Journal of Business & Economic Statistics 37 (3): 447–56.

我们可以将这种”放大并简化”的方法视为局部回归。局部回归在估计某个预测变量\(X\)与某个结果\(Y\)之间的关系时，允许这种关系在\(X\)的范围内自由变化²²。大多数研究人员选择使用某种局部回归来实现他们的断点回归设计，至少在他们拥有大样本时如此²³。

局部回归。一种针对预测变量不同值采用不同拟合方式的回归方法，根据观测值与该预测变量值的接近程度进行加权处理。

局部回归的具体做法是：对于每个\(X\)值，比如\(X=5\)，它通过运行自己特殊的\(X=5\)回归来获得\(Y\)的估计值。这个\(X=5\)回归会拟合某种特定于\(X=5\)的回归形状。如何特定？在估计回归时，它会给靠近\(X=5\)的观测值赋予更大的权重。这种特定的加权形式称为核函数。有时核函数是一个截断——靠近时权重为\(1\)，太远时权重为\(0\)。其他时候则是渐变的——距离越远权重越低，直到缓慢地降为\(0\)。断点回归中一个非常常用的核函数是三角核，如\(图20.8\)所示，正如其名——在您查看的值处权重最大，随着距离增加权重线性下降，直到在某点降为\(0\)。这个点就是”带宽”——即您能接受的最大距离。这形成了一个整洁的小三角形。

核函数：描述数值如何加权的函数。在局部回归的背景下，它描述了如何根据数值与关注点的接近程度进行加权。

从计算角度而言，局部回归需要大量计算资源。但概念上这非常简单：只需沿着\(X\)轴滚动，在途中基于当前最近的观测值创建一系列新回归，然后从中预测值，这就是您的曲线。

但应该在每个点运行哪种回归呢？

我们可能只取平均值而不运行回归，这种情况下我们得到的是核回归。核回归在断点回归中可能带来问题，特别是在结果变量在接近断点处存在任何趋势时。由于它只是取平均值，往往会将趋势扁平化，导致在断点处的预测将实际存在的部分趋势错误归因于断点处的跳跃。\(图20.9\)显示了一个这种问题可能发生的夸张示例。请注意断点左侧的预测值远低于所有实际数据，因为它忽略了上升趋势。

更常见的是，在断点回归应用中，我们会进行局部回归，其中最常见的形式是\(LOESS\)（我们在第\(4\)章”描述关系”和本章前面的\(图20.2\)中已经讨论过）。\(LOESS\)不是取局部均值，而是在每个点估计一个线性（“局部线性回归”）或二阶多项式（“局部多项式回归”或”局部二次回归”）回归。

局部回归方法允许\(X\)和\(Y\)之间的关系不是完全平坦的。根据其弯曲程度，我们选择的线性或二阶多项式函数可能不足以完全捕捉这种关系。然而，由于我们已经放大了很多，在如此狭窄的 \(x\) 轴范围内，任何关系都不太可能具有高度非线性。

当然，考虑到我们使用局部回归的计划，实际上不需要做任何特殊处理。虽然在整个数据范围内应用局部回归估计量可以让我们很好地了解运行变量与结果之间的整体关系，但如果我们只对估计断点回归感兴趣…我们实际上只对两个局部点感兴趣——刚好在断点左侧和刚好在右侧。因此，如果我们估计一个常规的交互线性回归，如\(方程（20.1）\)（用于局部线性回归）或\(方程（20.3）\)（用于局部多项式回归），但使用带宽（或者可能是缓慢衰减核的加权），我们就完成了。

让我们编写一些断点回归代码！我将用两种方式实现。首先，我将运行一个普通最小二乘模型，应用带宽和核权重，并使用异方差稳健标准误。然而，本章我们将讨论许多其他调整，将这个任务交给一个知道如何处理的包可能是个好主意。因此我们也将使用这些专门的命令，稍后我会讨论它们为我们做的一些额外工作。

在这个例子中，我们将使用\(Manacorda、Miguel \space \& \space Vigorito（2011）\)的《政府转移支付与政治支持》中的数据。该论文研究了乌拉圭的一项大型扶贫计划，该计划向大部分人口发放了可观的支票²⁴。他们感兴趣的是，收到这些资金是否使人们更有可能支持发放这些资金的新上台的中左翼政府。

Manacorda, Marco, Edward Miguel, and Andrea Vigorito. 2011. “Government Transfers and Political Support.” American Economic Journal: Applied Economics 3 (3): 1–28.

谁获得了支付？您的收入必须足够低。但他们并没有完全使用收入作为运行变量；那可能太容易操纵。相反，政府使用了一系列因素——住房、工作、报告收入、学校教育——并从中预测您的收入。然后，预测收入成为运行变量，处理分配基于是否低于某个断点。约 \(14\%\) 的人口最终获得了支付。

研究人员在收入断点附近对一群人进行了调查，以检查他们事后对政府的支持度。收入刚好低于断点的人是否比刚好高于断点的人更支持政府？

我们拥有的政府转移支付数据集中的预测收入变量已经预先中心化，因此断点位于\(0\)。然后，对政府的支持度取三个值：您认为他们比上一届政府更好\((1)\)、相同\((1/2)\)或更差\((0)\)。数据仅包括断点附近的个体——数据中的中心化收入变量仅从\(-0.02\)到\(0.02\)²⁵。

在估计我们的模型之前，让我们做一个几乎是强制性的图形断点回归检查，以确认在我们预期的位置确实存在某种断点²⁶。这是一个”分箱均值图”。有一些预编程的方法可以做到这一点，比如 \(rdrobust\) 包中的 \(rdplot\)（所有三种语言）或 \(Stata\) 中 \(binscatter\) 包的 \(binscatter\)，但这很简单，我们不妨手动完成并展示我们的绘图能力。

对于此类图表，通常需要：\(（1）\)将运行变量分割成若干区间（确保断点是两个区间的边界），\(（2）\)计算每个区间内结果变量的均值，\(（3）\)绘制结果图表，并在断点处添加垂直线以便识别断点位置。然后通常会用处理变量替代结果变量重复此过程，生成类似 \(图20.3\) 的图表。

library(tidyverse)

## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.5
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ ggplot2   3.5.1     ✔ tibble    3.2.1
## ✔ lubridate 1.9.4     ✔ tidyr     1.3.1
## ✔ purrr     1.0.4     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

gt <- causaldata::gov_transfers

# Use cut() to create bins, using breaks to make sure it breaks at 0
# (-15:15)*.02/15 gives 15 breaks from -.02 to .02
binned <- gt %>%
    mutate(Inc_Bins = cut(Income_Centered,breaks = (-15:15)*(.02/15))) %>%
    group_by(Inc_Bins) %>%
    summarize(Support = mean(Support),Income = mean(Income_Centered))

# Taking the mean of Income lets us plot data roughly at the bin midpoints
ggplot(binned, aes(x = Income, y = Support)) + 
    geom_line() + 
    # Add a cutoff line
    geom_vline(aes(xintercept = 0), linetype = 'dashed')

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn.objects as so
from causaldata import gov_transfers

d = gov_transfers.load_pandas().data
# cut at 0, and 15 places on either side
edges = np.linspace(-.02,.02,31)
d['Bins'] = pd.cut(d['Income_Centered'], bins = edges)
# Mean within bins
binned = d.groupby(['Bins']).agg('mean')

## <string>:2: FutureWarning: The default of observed=False is deprecated and will be changed to True in a future version of pandas. Pass observed=False to retain current behavior or observed=True to adopt the future default and silence this warning.

# And plot
fig = plt.figure()
(so.Plot(binned, x = 'Income_Centered', y = 'Support').on(fig).add(so.Line()).plot())

## <seaborn._core.plot.Plotter object at 0x000001C558234440>

fig.axes[0].axvline(0)

library(RStata)

dotext = "

  causaldata gov_transfers.dta, use clear download
  
  * Create bins. 15 bins on either side of the cutoff from -.02 to .02, plus 0, means we want \"steps\" of...
  local step = (.02)/15
  egen bins = cut(income_centered), at(-.02(`step').02)
  
  * Means within bins
  collapse (mean) support, by(bins)
  
  * And graph with a cutoff line
  twoway (line support bins) || (function y = 0, horiz range(support)), xti(\"Centered Income\") yti(\"Support\")

"
stata(dotext, data.in = NULL, data.out = TRUE)

## . 
## . 
## .   causaldata gov_transfers.dta, use clear download
## .   
## .   * Create bins. 15 bins on either side of the cutoff from -.02 to .02, plus 
## > 0, means we want "steps" of...
## .   local step = (.02)/15
## .   egen bins = cut(income_centered), at(-.02(`step').02)
## .   
## .   * Means within bins
## .   collapse (mean) support, by(bins)
## .   
## .   * And graph with a cutoff line
## .   twoway (line support bins) || (function y = 0, horiz range(support)), xti("
## > Centered Income") yti("Support")
## .

添加一些额外的主题设置后，我们得到\(图20.10\)²⁷。断点处明显存在一个跳跃，并且在断点左侧（由于需要低收入才能获得支付，因此处理组在左侧）出现了较高的支持度，这与根据断点右侧趋势预期的结果不一致。

现在我们已经有了图表的概念，实际上可以估计我们的模型了。我们将首先使用\(OLS\)和二阶多项式进行估计，然后使用三角核权重线性模型，将断点两侧的带宽限制在\(0.01\)范围内。在这种情况下，带宽并不是特别必要——数据已经限制在断点两侧\(0.02\)范围内——但这只是一个演示。

需要注意的是，使用\(OLS\)进行断点回归的第一步是将运行变量以断点为中心进行中心化。这很简单，只需创建变量 \(RunningVariable-Cutoff\)，并转换为您使用的任何语言。第二步是创建一个处理变量 \(RunningVariable < Cutoff\)（因为在这种情况下，处理应用于断点以下）。但在这个案例中，运行变量已经预先中心化，且低于断点的处理变量已经存在于数据中，因此我将省略这些部分。

library(tidyverse); library(modelsummary)

gt <- causaldata::gov_transfers

# Linear term and a squared term with "treated" interactions
m <- lm(Support ~ Income_Centered*Participation + I(Income_Centered^2)*Participation, data = gt)

# Add a triangular kernel weight
kweight <- function(x) {
    # To start at a weight of 0 at x = 0, and impose a bandwidth of .01, 
    # we need a "slope" of -1/.01 = 100, and to go in either direction use the absolute value
    w <- 1 - 100*abs(x)
    # if further away than .01, the weight is 0, not negative
    w <- ifelse(w < 0, 0, w)
    return(w)
}

# Run a linear version of the but with the weight
mw <- lm(Support ~ Income_Centered*Participation, data = gt,weights = kweight(Income_Centered))

# See the results with heteroskedasticity-robust SEs
msummary(list('Quadratic' = m, 'Linear with Kernel Weight' = mw), stars = c('*' = .1, '**' = .05, '***' = .01), vcov = 'robust')

## Warning in residuals^2/(1 - diaghat)^2: longer object length is not a multiple
## of shorter object length
## Warning in residuals^2/(1 - diaghat)^2: longer object length is not a multiple
## of shorter object length

	Quadratic	Linear with Kernel Weight
* p < 0.1, p < 0.05, * p < 0.01
(Intercept)	0.769***	0.819***
	(0.034)	(0.015)
Income_Centered	-11.567	-23.697***
	(8.101)	(3.219)
Participation	0.093**	0.033
	(0.044)	(0.021)
I(Income_Centered^2)	562.247
	(401.982)
Income_Centered × Participation	19.300*	26.594***
	(10.322)	(4.433)
Participation × I(Income_Centered^2)	-101.103
	(502.789)
Num.Obs.	1948	937
R2	0.036	0.041
R2 Adj.	0.034	0.038
AIC	1007.5	750.9
BIC	1046.6	775.2
Log.Lik.	-496.764	-370.470
F	14.128	55.653
RMSE	0.31	0.34
Std.Errors	HC3	HC3

import numpy as np
import statsmodels.formula.api as smf
from causaldata import gov_transfers

d = gov_transfers.load_pandas().data

# Run the polynomial model
m1 = smf.ols('''Support~Income_Centered*Participation + I(Income_Centered**2)*Participation''', d).fit()

# Create the kernel function
def kernel(x):
    # To start at a weight of 1 at x = 0,
    # and impose a bandwidth of .01, we need a "slope" of -1/.01 = 100 and to go in either direction use the absolute value
    w = 1 - 100*np.abs(x)
    # if further away than .01, the weight is 0, not negative
    w = np.maximum(0,w)
    return w

# Run the linear model with weights using wls
m2 = smf.wls('Support~Income_Centered*Participation', d,weights = kernel(d['Income_Centered'])).fit()

m1.summary()

OLS Regression Results
Dep. Variable:	Support	R-squared:	0.036
Model:	OLS	Adj. R-squared:	0.034
Method:	Least Squares	F-statistic:	14.63
Date:	周日, 14 9月 2025	Prob (F-statistic):	4.24e-14
Time:	02:04:22	Log-Likelihood:	-496.76
No. Observations:	1948	AIC:	1006.
Df Residuals:	1942	BIC:	1039.
Df Model:	5
Covariance Type:	nonrobust

	coef	std err	t	P>\|t\|	[0.025	0.975]
Intercept	0.7690	0.034	22.321	0.000	0.701	0.837
Income_Centered	-11.5666	7.777	-1.487	0.137	-26.819	3.686
Participation	0.0929	0.046	2.025	0.043	0.003	0.183
Income_Centered:Participation	19.3000	10.445	1.848	0.065	-1.185	39.785
I(Income_Centered ** 2)	562.2473	372.182	1.511	0.131	-167.672	1292.166
I(Income_Centered ** 2):Participation	-101.1025	500.196	-0.202	0.840	-1082.079	879.874

Omnibus:	335.257	Durbin-Watson:	1.997
Prob(Omnibus):	0.000	Jarque-Bera (JB):	526.668
Skew:	-1.246	Prob(JB):	4.32e-115
Kurtosis:	3.532	Cond. No.	9.81e+04

Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
[2] The condition number is large, 9.81e+04. This might indicate that there are
strong multicollinearity or other numerical problems.

m2.summary()

WLS Regression Results
Dep. Variable:	Support	R-squared:	0.041
Model:	WLS	Adj. R-squared:	0.040
Method:	Least Squares	F-statistic:	28.04
Date:	周日, 14 9月 2025	Prob (F-statistic):	9.52e-18
Time:	02:04:22	Log-Likelihood:	-inf
No. Observations:	1948	AIC:	inf
Df Residuals:	1944	BIC:	inf
Df Model:	3
Covariance Type:	nonrobust

	coef	std err	t	P>\|t\|	[0.025	0.975]
Intercept	0.8194	0.019	43.280	0.000	0.782	0.857
Income_Centered	-23.6967	4.455	-5.320	0.000	-32.433	-14.960
Participation	0.0335	0.026	1.292	0.196	-0.017	0.084
Income_Centered:Participation	26.5937	6.177	4.305	0.000	14.480	38.708

Omnibus:	854.936	Durbin-Watson:	2.003
Prob(Omnibus):	0.000	Jarque-Bera (JB):	4673.816
Skew:	-2.025	Prob(JB):	0.00
Kurtosis:	9.417	Cond. No.	1.21e+03

Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
[2] The condition number is large, 1.21e+03. This might indicate that there are
strong multicollinearity or other numerical problems.

dotxt <- "

  causaldata gov_transfers.dta, use clear download
  
  * Include the running variable, its square (by interaction with itself)
  * and interactions of both with treatment, and robust standard errors
  reg support i.participation##c.income_centered##c.income_centered, robust
  
  * Create the triangular kernel weight. To start at a weight of 0 at x = 0,
  * and impose a bandwidth of .01, we need a \"slope\" of -1/.01 = 100 and to go in either direction use the absolute value
  
  g w = 1 - 100*abs(income_centered)
  
  * if further away than .01, the weight is 0, not negative
  
  replace w = 0 if w < 0
  
  reg support i.participation##c.income_centered [aw = w], robust

"
stata(dotxt, data.in = NULL, data.out = TRUE)

## . 
## . 
## .   causaldata gov_transfers.dta, use clear download
## .   
## .   * Include the running variable, its square (by interaction with itself)
## .   * and interactions of both with treatment, and robust standard errors
## .   reg support i.participation##c.income_centered##c.income_centered, robust
## 
## Linear regression                               Number of obs     =      1,948
##                                                 F(5, 1942)        =      14.18
##                                                 Prob > F          =     0.0000
##                                                 R-squared         =     0.0363
##                                                 Root MSE          =     .31274
## 
## ------------------------------------------------------------------------------
##              |               Robust
##      support |      Coef.   Std. Err.      t    P>|t|     [95% Conf. Interval]
## -------------+----------------------------------------------------------------
## participat~n |
## Got a tra..  |   .0928553   .0436336     2.13   0.033     .0072818    .1784289
## income_cen~d |  -11.56663   8.070016    -1.43   0.152    -27.39344    4.260173
##              |
## participat~n#|
##           c. |
## income_cen~d |
## Got a tra..  |   19.29999   10.28725     1.88   0.061    -.8752289     39.4752
##              |
##           c. |
## income_cen~d#|
##           c. |
## income_cen~d |   562.2473   400.4888     1.40   0.161    -223.1858     1347.68
##              |
## participat~n#|
##           c. |
## income_cen~d#|
##           c. |
## income_cen~d |
## Got a tra..  |  -101.1025   501.1407    -0.20   0.840    -1083.933    881.7277
##              |
##        _cons |   .7689957   .0336602    22.85   0.000     .7029817    .8350096
## ------------------------------------------------------------------------------
## .   
## .   * Create the triangular kernel weight. To start at a weight of 0 at x = 0,
## .   * and impose a bandwidth of .01, we need a "slope" of -1/.01 = 100 and to g
## > o in either direction use the absolute value
## .   
## .   g w = 1 - 100*abs(income_centered)
## .   
## .   * if further away than .01, the weight is 0, not negative
## .   
## .   replace w = 0 if w < 0
## (1,011 real changes made)
## .   
## .   reg support i.participation##c.income_centered [aw = w], robust
## (sum of wgt is 459.9937024153769)
## 
## Linear regression                               Number of obs     =        937
##                                                 F(3, 933)         =      12.92
##                                                 Prob > F          =     0.0000
##                                                 R-squared         =     0.0415
##                                                 Root MSE          =     .30909
## 
## ------------------------------------------------------------------------------
##              |               Robust
##      support |      Coef.   Std. Err.      t    P>|t|     [95% Conf. Interval]
## -------------+----------------------------------------------------------------
## participat~n |
## Got a tra..  |   .0334818   .0441959     0.76   0.449    -.0532531    .1202166
## income_cen~d |  -23.69669   6.680465    -3.55   0.000    -36.80717   -10.58621
##              |
## participat~n#|
##           c. |
## income_cen~d |
## Got a tra..  |   26.59367   9.198351     2.89   0.004     8.541818    44.64553
##              |
##        _cons |   .8194073   .0321537    25.48   0.000     .7563053    .8825094
## ------------------------------------------------------------------------------
## .

由此，我们得到\(表20.1\)中的结果。在二阶多项式模型中，我们的估计是接收支付使对政府的支持度增加了\(9.3\)个百分点——不错！²⁸使用核权重的版本显示了一个统计不显著的\(3.3\)个百分点效应。第一个模型在我看来更合理——由于数据本身只包含断点附近的观测值，带宽已经相当窄，因此进一步限制数据没有必要。您还可以看到样本量从\(1948\)减少到\(937\)——这是因为权重为\(0\)的个体不被计入样本。减少偏差的潜力（通过减少远离断点的个体）是以样本量减少和精度降低为代价的。

现在我们完成了手动版本，可以尝试使用一些预打包命令，包括用于估计模型效应和制作断点回归风格图表的工具。对于所有三种语言，都有\(rdrobust\)包，它属于”\(RD \space packages\)“生态系统的一部分，包含许多有用的工具，我将在本章中陆续介绍²⁹。

\(rdrobust\)包中的\(rdrobust\)函数为我们做了几件事。首先，我们不需要选择带宽。它会为我们执行最优带宽选择程序。更多关于其工作原理的细节请参见”专业人士如何操作”部分。

其次，它将使用您选择的多项式阶数（默认为线性）和三角核来正确实现局部回归。

第三，如其名称所示，它将应用异方差稳健标准误。具体来说，它使用一个考虑了断点回归问题结构的稳健标准误估计量。它会寻找每个观测值在运行变量上的最近邻，并利用该信息来判断异方差程度并进行校正。

第四，它应用了偏差校正。最优带宽选择方法倾向于选择较宽的带宽，以帮助增加有效样本量并提高精度。然而，标准误的估计依赖于我们对误差项分布的理解。而这些宽带宽会扰乱这些分布的形状。因此，这里有一种方法来校正这种偏差，改进标准误的估计。更多关于此内容请参见”专业人士如何操作”部分。

library(tidyverse); library(rdrobust)

gt <- causaldata::gov_transfers

# Estimate regression discontinuity and plot it
m <- rdrobust(gt$Support, gt$Income_Centered, c = 0, masspoints = "adjust")

## Warning in rdrobust(gt$Support, gt$Income_Centered, c = 0, masspoints =
## "adjust"): Mass points detected in the running variable.

summary(m)

## Sharp RD estimates using local polynomial regression.
## 
## Number of Obs.                 1948
## BW type                       mserd
## Kernel                   Triangular
## VCE method                       NN
## 
## Number of Obs.                 1127          821
## Eff. Number of Obs.             291          194
## Order est. (p)                    1            1
## Order bias  (q)                   2            2
## BW est. (h)                   0.005        0.005
## BW bias (b)                   0.010        0.010
## rho (h/b)                     0.509        0.509
## Unique Obs.                     841          639
## 
## =============================================================================
##         Method     Coef. Std. Err.         z     P>|z|      [ 95% C.I. ]       
## =============================================================================
##   Conventional     0.025     0.062     0.396     0.692    [-0.098 , 0.147]     
##         Robust         -         -     0.624     0.533    [-0.097 , 0.188]     
## =============================================================================

# Note, by default, rdrobust and rdplot use different numbers of polynomial terms. You can set the p option to standardize them.
rdplot(gt$Support, gt$Income_Centered)

## [1] "Mass points detected in the running variable."

import pandas as pd
from rdrobust import rdrobust, rdplot
from causaldata import gov_transfers

gt = gov_transfers.load_pandas().data

# Estimate regression discontinuity and plot it
m = rdrobust(gt['Support'], gt['Income_Centered'], c = 0)

## Mass points detected in the running variable.
## Mass points detected in the running variable.

print(m)

## Call: rdrobust
## Number of Observations:                  1948
## Polynomial Order Est. (p):                  1
## Polynomial Order Bias (q):                  2
## Kernel:                            Triangular
## Bandwidth Selection:                    mserd
## Var-Cov Estimator:                         NN
## 
##                                 Left      Right
## ------------------------------------------------
## Number of Observations          1127        821
## Number of Unique Obs.            841        639
## Number of Effective Obs.         291        194
## Bandwidth Estimation           0.005      0.005
## Bandwidth Bias                  0.01       0.01
## rho (h/b)                      0.509      0.509
## 
## Method             Coef.     S.E.   t-stat    P>|t|       95% CI      
## -------------------------------------------------------------------------
## Conventional       0.025    0.062    0.396   6.920e-01    [-0.098, 0.147]
## Robust                 -        -    0.624   5.328e-01    [-0.097, 0.188]

# Note, by default, rdrobust and rdplot use different numbers
# of polynomial terms. You can set the p option to standardize them.
rdplot(gt['Support'], gt['Income_Centered'])

## Mass points detected in the running variable.
## Call: rdplot
## Number of Observations:                  1948
## Kernel:                               Uniform
## Polynomial Order Est. (p):                  4
## 
##                                 Left      Right
## ------------------------------------------------
## Number of Observations          1127        821
## Number of Effective Obs         1127        821
## Bandwith poly. fit (h)          0.02       0.02
## Number of bins scale               1          1
## Bins Selected                     35         35
## Average Bin Length             0.001      0.001
## Median Bin Length              0.001      0.001
## IMSE-optimal bins                3.0        7.0
## Mimicking Variance bins         35.0       35.0
## 
## Relative to IMSE-optimal:
## Implied scale                 11.667        5.0
## WIMSE variance weight          0.001      0.008
## WIMSE bias weight              0.999      0.992

dotxt <- "
  * STATA CODE
  * If necessary: ssc install rdrobust
  causaldata gov_transfers.dta, use clear download
  
  * Run the RDD model and plot it. Note, by default, 
  * rdrobust and rdplot use different
  * numbersof polynomial terms. 
  * You can set the p() option to standardize them.
  rdrobust support income_centered, c(0)
  rdplot support income_centered, c(0)
"
stata(dotxt, data.in = NULL, data.out = TRUE)

## . 
## .   * STATA CODE
## .   * If necessary: ssc install rdrobust
## .   causaldata gov_transfers.dta, use clear download
## .   
## .   * Run the RDD model and plot it. Note, by default, 
## .   * rdrobust and rdplot use different
## .   * numbersof polynomial terms. 
## .   * You can set the p() option to standardize them.
## .   rdrobust support income_centered, c(0)
## Mass points detected in the running variable.
## 
## Sharp RD estimates using local polynomial regression.
## 
##       Cutoff c = 0 | Left of c  Right of c            Number of obs =       194
## > 8
## -------------------+----------------------            BW type       =      mser
## > d
##      Number of obs |      1127         821            Kernel        = Triangula
## > r
## Eff. Number of obs |       291         194            VCE method    =         N
## > N
##     Order est. (p) |         1           1
##     Order bias (q) |         2           2
##        BW est. (h) |     0.005       0.005
##        BW bias (b) |     0.010       0.010
##          rho (h/b) |     0.509       0.509
##         Unique obs |       841         639
## 
## Outcome: support. Running variable: income_centered.
## -------------------------------------------------------------------------------
## > -
##             Method |   Coef.    Std. Err.    z     P>|z|    [95% Conf. Interval
## > ]
## -------------------+-----------------------------------------------------------
## > -
##       Conventional |   .0247     .06236   0.3961   0.692    -.09752      .14692
## > 4
##             Robust |     -          -     0.6238   0.533   -.097391      .18832
## > 5
## -------------------------------------------------------------------------------
## > -
## Estimates adjusted for mass points in the running variable.
## .   rdplot support income_centered, c(0)
## Mass points detected in the running variable.
## 
## RD Plot with evenly spaced mimicking variance number of bins using polynomial r
## > egression.
## 
##          Cutoff c = 0 | Left of c  Right of c        Number of obs  =       194
## > 8
## ----------------------+----------------------        Kernel         =    Unifor
## > m
##         Number of obs |      1127         821
##    Eff. Number of obs |      1126         820
##   Order poly. fit (p) |         4           4
##      BW poly. fit (h) |     0.020       0.020
##  Number of bins scale |     1.000       1.000
## 
## Outcome: support. Running variable: income_centered.
## ---------------------------------------------
##                       | Left of c  Right of c
## ----------------------+----------------------
##         Bins selected |        35          35
##    Average bin length |     0.001       0.001
##     Median bin length |     0.001       0.001
## ----------------------+----------------------
##     IMSE-optimal bins |         3           7
##   Mimicking Var. bins |        35          35
## ----------------------+----------------------
## Rel. to IMSE-optimal: | 
##         Implied scale |    11.667       5.000
##     WIMSE var. weight |     0.001       0.008
##     WIMSE bias weight |     0.999       0.992
## ---------------------------------------------

如果您运行此代码，会发现\(rdrobust\)的输出与我们手动操作时的结果有较大差异。结果不显著且与原始论文方向相反³⁰。这归结于\(rdrobust\)选择了一个较窄的带宽——如果强制将带宽设置为整个范围（考虑到原本范围已经很窄，这可能是合理的）——您会得到与\(表20.1\)中二次模型几乎相同的结果。做出这些决定很困难！仅仅因为我们有预制包并不意味着我们可以让它做出所有决定。动动脑筋，通过阅读文档来探索函数中的选项。

20.2.2 模糊断点回归

到底有多少种处理干预是完全由一个运行变量的临界值所决定的呢？当然是有一些的。但我们生活在现实世界里。即使在政策本应基于临界值来分配的情况下，肯定也可能会有少数人，虽然处于临界值之上或之下，却设法得到了错误的处理。抛开这一点不谈，还有很多场景，政策未必是要基于临界值来分配的，但出于某种原因，在临界值处，接受处理的概率会显著跃升。或许你想评估拥有高中学位对某个结果的影响。你拥有高中学位的概率在你\(18\)岁生日之后的那个夏天会显著上升。但有相当多的人在 \(17\) 岁或 \(19\) 岁毕业，或者不是在夏天毕业，又或者根本就没毕业。在那个夏天的临界值处，处理概率（拥有高中学位的概率）会有一个跃升，但并非从 \(0\%\) 跃升到 \(100\%\)。我们可以在 \(图20.11\) 中看到，处理概率在运行变量上可能呈现出的样子。

我已经在”它是如何工作的”部分涵盖了基本概念。但我们实际上该如何实施模糊断点回归设计呢？

我们可以通过修改\(图20.1\)中的因果图来弄清楚这一点。我们需要考虑的唯一变化是，除了临界值之外，还有其他一些方式能让运行变量和处理（干预）产生关联，这就是为什么我们在\(图20.11\)的非临界值点能看到处理概率的变化。我们也可能会有一些非运行变量的处理干预决定因素。

我们在\(图20.12\)中得到了修改后的因果图³¹。从这个图中我们能得到什么呢？

嗯，需要注意的一点是，识别效应的思路基本上还是一样的。如果我们控制了运行变量，那么从高于临界值到结果就没有后门路径。通过排除由运行变量在临界值之外任何地方所驱动的变化，我们就只剩下没有后门路径的变化，而这些变化能识别处理（干预）对结果的效应。

唯一真正的变化是：我们不能再简单地将数据限定在临界值附近的区域，然后就认为在控制运行变量方面大功告成了。那样做会导致我们低估效应。如果我们看到的处理率跃升只有，比如说\(15\)个百分点，那么我们应该只预期看到\(15\%\)的效应。我们必须将其调整回完整的规模。

值得庆幸的是，这可以通过应用工具变量来实现，就像在第\(19\)章中那样，使用与精确断点回归设计基本相同的断点回归方程，比如\(方程（20.1）\)或\(方程（20.3）\)。只不过这些方程现在变成了我们的第二阶段方程。我们的第一阶段使用高于临界值作为接受处理（干预）的工具变量（并且使用与高于临界值的交互项作为与接受处理（干预）的交互项的工具变量）。

这会完全按照我们想要的方式对估计值进行调整。在其基本形式中，工具变量法是用工具变量对结果的效应除以工具变量对内生/处理变量的效应³²。换句话说，我们是在调整高于临界值对结果的效应（也就是我们从典型的精确断点回归模型中得到的效应），以考虑到高于临界值仅导致处理概率部分上升这一事实。

就像精确断点回归一样，我们可以用代码实现模糊断点回归，要么使用我们的标准工具（和第\(19\)章中使用的相同的工具变量估计器），要么使用专门的断点回归专用工具。我们在这里两种方法都试试。

为此，我们将使用\(Fetter（2013）\) 的数据。\(Fetter\)的研究也让我们有机会看到在一个稍微没那么干净的环境中进行的断点回归分析。\(Fetter\)主要感兴趣的问题是，\(20\)世纪中期美国住房拥有率的上升，有多少是源于政府提供的抵押贷款补贴。

Fetter, Daniel K. 2013. “How Do Mortgage Subsidies Affect Home Ownership? Evidence from the Mid-Century GI Bills.” American Economic Journal: Economic Policy 5 (2): 111–47.

他是如何研究这个问题的呢？他关注的是那些年龄刚好符合参加像第二次世界大战或朝鲜战争等重大战争的退伍军人条件的人。任何参加过这些战争的退伍军人都能获得特殊的抵押贷款补贴。

这和断点回归有什么关系呢？参军是有年龄要求的。如果你出生得太晚，晚了一年，没能参军去参加朝鲜战争，那么你就无法获得这些抵押贷款补贴（或者说，至少符合条件的退伍军人要少得多）。要是出生得刚刚好呢？你就可以参军并获得补贴！所以，这里存在一个基于出生年份的断点回归。是出生得刚刚好，还是稍微晚了一点呢？

当然，并非每一个在特定年份出生的人都会参军³³。所以，有资格获得抵押贷款补贴这一处理干预），只会适用于那些在恰当时间出生的一部分人。处理率从\(0\%\)跃升到…… 高于\(0\%\)，但不是\(100\%\)。这就是模糊的情况！

那我为什么说这比其他一些例子”没那么干净”呢？因为这不是一个明确设计的临界值。在谁能得到处理干预方面，存在更多的选择，而且，哎呀，有些人可能在未达法定年龄的情况下就参军了。在这里，我们必须非常仔细地思考，是否相信使用断点回归所需的那些假设是成立的。我认为这里的设计相当合理，但我对它的看法，不像对那些处理干预分配涉及更少可变因素的情况的看法那样绝对可靠³⁴。

作者还得进行一些检验和测试，而对于那些研究故事更简单的人来说，可能就不用做这些。比如，虽然我不会在下面的代码中展示，但他进行了一项安慰剂检验，即使用其他时代的数据重复分析 —— 在那些时代，退伍军人身份不会对获得抵押贷款补贴产生影响，结果发现在那些时候没有效果。你们即将看到的分析包含了出生季节、种族以及他们出生的美国州份这些控制变量 —— 对于精确断点回归来说，控制变量通常不是必需的，但在这里纳入它们是有意义的。如果设计本身不能为你完成所有工作，那么天哪，我猜你就得更努力一点了。

让我们把分析用代码写出来。首先，我们会用第\(19\)章用过的工具变量代码自己来做。我们使用的运行变量是出生季度\(（qob）\)，这个变量已经以有资格因参加朝鲜战争而获得抵押贷款补贴所需的出生季度为中心进行了调整\(（qob\_minus\_kw）\)。这决定了你是否是朝鲜战争或第二次世界大战的退伍军人\(（vet\_wwko）\)。所有这些都是为了看看退伍军人身份（以及随之而来的抵押贷款补贴）是否会影响你的住房拥有状况\(（home\_ownership）\)—— 记住，这里的目标是探寻那些抵押贷款补贴的影响。还有对是否为白人或非白人，以及出生州的控制变量。

为了简化问题，在进行估计之前，我们会应用一个带宽，跳过核加权，只使用线性模型。不过，如果你想添加这些内容，我们之前用于整合这些内容的代码在这里仍然有效。

library(tidyverse); library(fixest); library(modelsummary)

vet <- causaldata::mortgages
# Create an "above-cutoff" variable as the instrument
vet <- vet %>% mutate(above = qob_minus_kw > 0)
# Impose a bandwidth of 12 quarters on either side
vet <- vet %>% filter(abs(qob_minus_kw) < 12)
# Control for race + fixed effect controls + Instrument our standard RDD + with being above the cutoff
m <- feols(home_ownership ~ nonwhite | bpl + qob | qob_minus_kw*vet_wwko ~ qob_minus_kw*above, se = 'hetero',data = vet) 
# And look at the results
msummary(m, stars = c('*' = .1, '**' = .05, '***' = .01))

	(1)
* p < 0.1, p < 0.05, * p < 0.01
fit_qob_minus_kw	-0.007***
	(0.002)
fit_vet_wwko	0.170***
	(0.046)
fit_qob_minus_kw × vet_wwko	-0.003
	(0.003)
nonwhite	-0.190***
	(0.007)
Num.Obs.	56901
R2	0.053
R2 Adj.	0.052
R2 Within	0.037
R2 Within Adj.	0.037
AIC	68659.5
BIC	69187.5
RMSE	0.44
Std.Errors	Heteroskedasticity-robust
FE: bpl	X
FE: qob	X

import pandas as pd
from linearmodels.iv import IV2SLS
from causaldata import mortgages

d = mortgages.load_pandas().data
# Create an above variable as an instrument
# Apply a bandwidth of 12 quarters on either side
d = (d.assign(above = d['qob_minus_kw'] > 0).query('abs(qob_minus_kw) < 12'))
# Now we estimate!
# Add nonwhite, bpl, qob controls and instrumnet for the vet_wwko interaction with the above-cutoff interaction
m = IV2SLS.from_formula('''home_ownership ~ nonwhite + C(bpl) + C(qob) + [qob_minus_kw*vet_wwko ~ qob_minus_kw*above]''', data = d)

# With robust standard errors
m.fit(cov_type = 'robust')

IV-2SLS Estimation Summary
Dep. Variable:	home_ownership	R-squared:	0.0532
Estimator:	IV-2SLS	Adj. R-squared:	0.0522
No. Observations:	56901	F-statistic:	2.536e+04
Date:	周日, 9月 14 2025	P-value (F-stat)	0.0000
Time:	02:04:27	Distribution:	chi2(59)
Cov. Estimator:	robust

Parameter Estimates
	Parameter	Std. Err.	T-stat	P-value	Lower CI	Upper CI
nonwhite	-0.1904	0.0069	-27.646	0.0000	-0.2039	-0.1769
C(bpl)[Alabama]	0.2775	0.0208	13.343	0.0000	0.2367	0.3182
C(bpl)[Alaska]	0.3695	0.0811	4.5546	0.0000	0.2105	0.5286
C(bpl)[Arizona]	0.3001	0.0368	8.1525	0.0000	0.2280	0.3723
C(bpl)[Arkansas]	0.2758	0.0210	13.144	0.0000	0.2347	0.3170
C(bpl)[California]	0.2474	0.0246	10.069	0.0000	0.1992	0.2955
C(bpl)[Colorado]	0.2805	0.0292	9.6234	0.0000	0.2234	0.3377
C(bpl)[Connecticut]	0.1439	0.0288	4.9906	0.0000	0.0874	0.2004
C(bpl)[Delaware]	0.3632	0.0489	7.4263	0.0000	0.2674	0.4591
C(bpl)[District of Columbia]	0.1997	0.0349	5.7172	0.0000	0.1312	0.2681
C(bpl)[Florida]	0.3115	0.0263	11.823	0.0000	0.2598	0.3631
C(bpl)[Georgia]	0.2736	0.0199	13.766	0.0000	0.2346	0.3126
C(bpl)[Hawaii]	0.2191	0.0333	6.5771	0.0000	0.1538	0.2843
C(bpl)[Idaho]	0.2891	0.0367	7.8749	0.0000	0.2171	0.3610
C(bpl)[Illinois]	0.2243	0.0239	9.3871	0.0000	0.1775	0.2712
C(bpl)[Indiana]	0.3226	0.0247	13.076	0.0000	0.2742	0.3709
C(bpl)[Iowa]	0.2475	0.0252	9.8094	0.0000	0.1981	0.2970
C(bpl)[Kansas]	0.2727	0.0267	10.214	0.0000	0.2204	0.3251
C(bpl)[Kentucky]	0.2476	0.0218	11.354	0.0000	0.2049	0.2903
C(bpl)[Louisiana]	0.2888	0.0224	12.903	0.0000	0.2450	0.3327
C(bpl)[Maine]	0.2640	0.0307	8.6073	0.0000	0.2039	0.3241
C(bpl)[Maryland]	0.2603	0.0259	10.044	0.0000	0.2095	0.3111
C(bpl)[Massachusetts]	0.1243	0.0262	4.7524	0.0000	0.0730	0.1755
C(bpl)[Michigan]	0.3537	0.0240	14.735	0.0000	0.3067	0.4008
C(bpl)[Minnesota]	0.2886	0.0254	11.350	0.0000	0.2388	0.3385
C(bpl)[Mississippi]	0.2635	0.0207	12.705	0.0000	0.2228	0.3041
C(bpl)[Missouri]	0.2677	0.0240	11.141	0.0000	0.2206	0.3148
C(bpl)[Montana]	0.2636	0.0357	7.3917	0.0000	0.1937	0.3335
C(bpl)[Nebraska]	0.2536	0.0279	9.0796	0.0000	0.1989	0.3084
C(bpl)[Nevada]	0.3520	0.0779	4.5176	0.0000	0.1993	0.5047
C(bpl)[New Hampshire]	0.1739	0.0416	4.1774	0.0000	0.0923	0.2555
C(bpl)[New Jersey]	0.1564	0.0246	6.3593	0.0000	0.1082	0.2046
C(bpl)[New Mexico]	0.2974	0.0320	9.3025	0.0000	0.2348	0.3601
C(bpl)[New York]	0.1375	0.0225	6.1048	0.0000	0.0933	0.1816
C(bpl)[North Carolina]	0.2429	0.0199	12.210	0.0000	0.2040	0.2819
C(bpl)[North Dakota]	0.2543	0.0323	7.8845	0.0000	0.1911	0.3176
C(bpl)[Ohio]	0.2978	0.0227	13.095	0.0000	0.2532	0.3423
C(bpl)[Oklahoma]	0.2696	0.0240	11.246	0.0000	0.2227	0.3166
C(bpl)[Oregon]	0.2612	0.0341	7.6482	0.0000	0.1943	0.3281
C(bpl)[Pennsylvania]	0.2216	0.0230	9.6273	0.0000	0.1765	0.2667
C(bpl)[Rhode Island]	0.1618	0.0354	4.5743	0.0000	0.0925	0.2312
C(bpl)[South Carolina]	0.2713	0.0216	12.569	0.0000	0.2290	0.3136
C(bpl)[South Dakota]	0.2729	0.0321	8.4904	0.0000	0.2099	0.3359
C(bpl)[Tennessee]	0.2477	0.0215	11.534	0.0000	0.2056	0.2897
C(bpl)[Texas]	0.2922	0.0198	14.728	0.0000	0.2533	0.3311
C(bpl)[United States, ns]	0.1570	0.0458	3.4305	0.0006	0.0673	0.2468
C(bpl)[Utah]	0.3116	0.0333	9.3667	0.0000	0.2464	0.3768
C(bpl)[Vermont]	0.2164	0.0438	4.9416	0.0000	0.1305	0.3022
C(bpl)[Virginia]	0.2287	0.0218	10.474	0.0000	0.1859	0.2715
C(bpl)[Washington]	0.2595	0.0287	9.0523	0.0000	0.2033	0.3157
C(bpl)[West Virginia]	0.1951	0.0243	8.0277	0.0000	0.1475	0.2427
C(bpl)[Wisconsin]	0.2145	0.0237	9.0352	0.0000	0.1680	0.2610
C(bpl)[Wyoming]	0.1914	0.0435	4.4034	0.0000	0.1062	0.2766
C(qob)[T.2]	-0.0080	0.0053	-1.5095	0.1312	-0.0183	0.0024
C(qob)[T.3]	-0.0110	0.0052	-2.0936	0.0363	-0.0212	-0.0007
C(qob)[T.4]	-0.0022	0.0054	-0.4097	0.6820	-0.0128	0.0083
qob_minus_kw	-0.0072	0.0018	-4.0342	0.0001	-0.0106	-0.0037
vet_wwko	0.1702	0.0459	3.7067	0.0002	0.0802	0.2602
qob_minus_kw:vet_wwko	-0.0029	0.0026	-1.0822	0.2792	-0.0080	0.0023

Endogenous: qob_minus_kw, vet_wwko, qob_minus_kw:vet_wwko
Instruments: qob_minus_kw, above, qob_minus_kw:above
Robust Covariance (Heteroskedastic)
Debiased: False
id: 0x1c5584d56a0

library(RStata)

dotext <- "

* We will be using ivreghdfe, 
* which must be installed with ssc install ivreghdfe
* along with ftools, ranktest, reghdfe, and ivreg2

causaldata mortgages.dta, use clear download

* Create an above variable as an instrument
g above = qob_minus_kw > 0

* Impose a bandwidth of 12 quarters on either side
keep if abs(qob_minus_kw) < 12

* Regress, using above as an instrument for veteran status.  
* Note that qob_minus_kw  by itself doesn't need instrumentation so we separate it. Usually I advise against it,  but in this case this is * * easier if we just make our  own interaction with g (generate)

g interaction_vet = vet_wwko*qob_minus_kw
g interaction_above = above*qob_minus_kw

ivreghdfe home_ownership nonwhite qob_minus_kw (vet_wwko interaction_vet = above interaction_above), absorb(bpl qob) robust

* We can also use regular ol' ivregress but it will be slower
encode bpl, g(bpl_num)
ivregress 2sls home_ownership nonwhite qob_minus_kw i.bpl_num i.qob  (vet_wwko interaction_vet = above interaction_above), robust

"
stata(dotext, data.in = NULL, data.out = TRUE)

## . 
## . 
## . * We will be using ivreghdfe, 
## . * which must be installed with ssc install ivreghdfe
## . * along with ftools, ranktest, reghdfe, and ivreg2
## . 
## . causaldata mortgages.dta, use clear download
## . 
## . * Create an above variable as an instrument
## . g above = qob_minus_kw > 0
## . 
## . * Impose a bandwidth of 12 quarters on either side
## . keep if abs(qob_minus_kw) < 12
## (157,243 observations deleted)
## . 
## . * Regress, using above as an instrument for veteran status.  
## . * Note that qob_minus_kw  by itself doesn't need instrumentation so we separa
## > te it. Usually I advise against it,  but in this case this is * * easier if w
## > e just make our  own interaction with g (generate)
## . 
## . g interaction_vet = vet_wwko*qob_minus_kw
## . g interaction_above = above*qob_minus_kw
## . 
## . ivreghdfe home_ownership nonwhite qob_minus_kw (vet_wwko interaction_vet = ab
## > ove interaction_above), absorb(bpl qob) robust
## last estimates not found
## r(301);
## . 
## . * We can also use regular ol' ivregress but it will be slower
## . encode bpl, g(bpl_num)
## . ivregress 2sls home_ownership nonwhite qob_minus_kw i.bpl_num i.qob  (vet_wwk
## > o interaction_vet = above interaction_above), robust
## .

这些结果表明，在这个边际上，退伍军人身份使住房拥有率提高了\(17\)个百分点。还不错。我们也可以应用上次完全相同的断点回归绘图代码（没有专门的 “模糊” 绘图，不过现在我们会把典型的断点回归图称为”第二阶段”图），并且也应用该绘图代码，将处理变量作为结果（即”第一阶段”），以确保分配变量如预期那样发挥作用。这就得到了\(图20.13\)。³⁵看起来不错！我们确实看到，在临界值附近，符合条件的退伍军人比例有一个急剧的变化，并且在同一点上，住房拥有率也存在一个断点。模糊断点回归估计值是\(图20.13（b）\)中临界值处的跃升除以\(图20.13（a）\)中临界值处的跃升。

接下来，我们可以像之前一样，使用 \(rdrobust\) 函数来进行局部多项式回归、核加权、便捷绘图等所有这些好用的操作。我们只需要将 \(vet\_wwko\) 用作我们的”处理状态指示器”，并把它提供给 \(fuzzy\) 选项。

library(tidyverse); library(rdrobust)

vet <- causaldata::mortgages
# It will apply a bandwidth anyway, but having it check the whole bandwidth space will be slow. So let's
# pre-limit it to a reasonable range of 12 quarters
vet <- vet %>%
    filter(abs(qob_minus_kw) <= 12)

# Create our matrix of controls
controls <- vet %>%
    select(nonwhite, bpl, qob) %>%
    mutate(qob = factor(qob))
# and make it a matrix with dummies
conmatrix <- model.matrix(~., data = controls)

# This is fairly slow due to the controls, beware!
m <- rdrobust(vet$home_ownership, vet$qob_minus_kw, fuzzy = vet$vet_wwko, c = 0,covs = conmatrix)

## Warning in rdrobust(vet$home_ownership, vet$qob_minus_kw, fuzzy = vet$vet_wwko,
## : Multicollinearity issue detected in covs. Redundant covariates dropped.

## Warning in rdrobust(vet$home_ownership, vet$qob_minus_kw, fuzzy = vet$vet_wwko,
## : Mass points detected in the running variable.

summary(m)

## Covariate-adjusted Fuzzy RD estimates using local polynomial regression.
## 
## Number of Obs.                56901
## BW type                       mserd
## Kernel                   Triangular
## VCE method                       NN
## 
## Number of Obs.                28776        28125
## Eff. Number of Obs.            6911         6756
## Order est. (p)                    1            1
## Order bias  (q)                   2            2
## BW est. (h)                   3.387        3.387
## BW bias (b)                   5.354        5.354
## rho (h/b)                     0.633        0.633
## Unique Obs.                      12           12
## 
## First-stage estimates.
## 
## =============================================================================
##         Method     Coef. Std. Err.         z     P>|z|      [ 95% C.I. ]       
## =============================================================================
##   Conventional    -0.248     0.018   -13.637     0.000    [-0.284 , -0.212]    
##         Robust         -         -   -12.968     0.000    [-0.334 , -0.247]    
## =============================================================================
## 
## Treatment effect estimates.
## 
## =============================================================================
##         Method     Coef. Std. Err.         z     P>|z|      [ 95% C.I. ]       
## =============================================================================
##   Conventional    -1.514     0.133   -11.378     0.000    [-1.775 , -1.253]    
##         Robust         -         -   -10.479     0.000    [-2.039 , -1.396]    
## =============================================================================

import pandas as pd
from rdrobust import rdrobust
from causaldata import mortgages

vet = mortgages.load_pandas().data
# It will apply a bandwidth anyway, but having it check the whole bandwidth space will be slow. So let's
# pre-limit it to a reasonable range of 12 quarters
vet = vet.query('abs(qob_minus_kw) <= 12')
# Create our matrix of controls with dummies
controls = pd.concat([d[['nonwhite']], pd.get_dummies(d[['bpl']])], axis = 1)
                  
d['qob'] = pd.Categorical(d['qob'])
# Drop one since we already have full rank with bpl (we'd also drop_first with bpl if we did add_constant)
controls = pd.concat([controls, pd.get_dummies(d[['qob']], drop_first = True)], axis = 1)

# As of this writing, actually adding those controls is unworkably slow. 
# But maybe by the time you read this  you can add them with covs = controls
m = rdrobust(vet['home_ownership'], vet['qob_minus_kw'],fuzzy = vet['vet_wwko'], c = 0)

## Mass points detected in the running variable.
## Mass points detected in the running variable.

print(m)

## Call: rdrobust
## Number of Observations:                 56901
## Polynomial Order Est. (p):                  1
## Polynomial Order Bias (q):                  2
## Kernel:                            Triangular
## Bandwidth Selection:                    mserd
## Var-Cov Estimator:                         NN
## 
##                                 Left      Right
## ------------------------------------------------
## Number of Observations         28776      28125
## Number of Unique Obs.             12         12
## Number of Effective Obs.        6911       6756
## Bandwidth Estimation           2.797      2.797
## Bandwidth Bias                 5.225      5.225
## rho (h/b)                      0.535      0.535
## 
## Method             Coef.     S.E.   t-stat    P>|t|       95% CI      
## -------------------------------------------------------------------------
## Conventional       1.879     3.35    0.561   5.750e-01    [-4.688, 8.445]
## Robust                 -        -     1.26   2.076e-01   [-2.817, 12.963]

library(RStata)

dotext = "

* This uses rdrobust, which must be installed with ssc install rdrobust
causaldata mortgages.dta, use clear download

* It will apply a bandwidth anyway, but having it check the whole bandwidth space will be slow. 
* So let's pre-limit it to a reasonable range of 12 quarters

keep if abs(qob_minus_kw) <= 12

* And run rdrobust
rdrobust home_ownership qob_minus_kw, c(0) fuzzy(vet_wwko) covs(nonwhite bpl qob)

"
stata(dotext, data.in = NULL, data.out = TRUE)

## . 
## . 
## . * This uses rdrobust, which must be installed with ssc install rdrobust
## . causaldata mortgages.dta, use clear download
## . 
## . * It will apply a bandwidth anyway, but having it check the whole bandwidth s
## > pace will be slow. 
## . * So let's pre-limit it to a reasonable range of 12 quarters
## . 
## . keep if abs(qob_minus_kw) <= 12
## (157,243 observations deleted)
## . 
## . * And run rdrobust
## . rdrobust home_ownership qob_minus_kw, c(0) fuzzy(vet_wwko) covs(nonwhite bpl 
## > qob)
## no observations
## r(2000);
## .

它没有包含在代码中，但我们可以很容易地应用 \(rdplot\)（也许也可以将处理变量作为因变量来应用它，只是为了确保它看起来符合我们的预期）。

\(rdrobust\) 估计值与我们之前做的线性工具变量估计值有很大不同，退伍军人身份使住房拥有率提高了 \(51.8\) 个百分点，而不是 \(17\) 个！这差别相当大。了解这些估计方法何时会出现这样的差异是有用的。我们应该相信哪一个呢？好吧，让我们仔细想想它们为何不同：\(rdrobust\) 使用二阶多项式而不是线性回归，它使用局部回归，并且使用更窄的带宽（它在临界值两侧各选取 \(3.39\) 个出生季度，而不是 \(12\) 个）。

在思考应该相信哪一个时，我们的第一步必须是询问其中一个模型是否做出了更有可能准确的假设。这些都是你在描述自己的研究设计时必须去辩护的内容。例如，如果我们认为线性是合理的，那么线性模型可能比二阶多项式更少噪声，是更精确的估计。\(图20.13\) 看起来是线性的吗？如果我们认为更宽的带宽可能包含太多不可比的人，那么我们可能会认为更窄的带宽更好。我们是否认为临界值两侧各 \(12\) 个月的范围太远而不可比？

如今，进行断点回归的研究者通常会观察估计值对不同带宽的敏感性，多次运行模型并报告效应如何变化，让研究的读者去思考应该将范围限制在离临界值多近的地方。³⁶我在”回归弯折”部分讨论的\(Bana、Bedard \space \& \space Rossin-Slater（2020）\)的论文就做了这样的事。

20.2.3 安慰剂检验

断点回归令人惊叹的一点是，它关闭了所有的后门路径，甚至是那些通过我们无法测量的变量的后门路径。这就是核心思想 —— 因为我们在运行变量的一个极窄窗口内隔离出了变异，所以声称在临界值处唯一发生变化的就是处理干预，进而还有处理干预所影响的任何事物（比如结果），是相当可信的。

更方便的是，断点回归带来的这个极其美妙的结果，为我们提供了一种相当直接的方法来检验它是否成立。只需对处理干预本不应该影响的某个事物运行你的断点回归模型即可。我们通常可能用作控制变量的任何事物都能用于这个目的³⁷。如果我们确实发现了效应，那我们最初的断点回归设计可能就不太对。或许出于某种原因，在临界值处的情况并不像我们希望的那样具有随机性。

这里没太多要讲的。方程、代码，所有东西对于这个检验和常规的断点回归来说完全一样。我们只是把实际的结果变量换成处理干预本不应该影响的某个其他变量。在很多情况下，我们可以用一长串可能用作控制变量的变量来进行尝试。

\(图20.14\)用《政府转移支付与政治支持》论文中的两个变量展示了这个检验。这里没有效应！正是我们所期望的。

关于这些安慰剂检验，有一点很重要需要记住：由于你可以对一长串潜在的安慰剂结果进行检验，所以仅靠随机偶然，你就有可能发现一些非零效应。这种情况会发生。如果你检验了一长串变量并发现了一些差异，这对你的设计来说并非致命问题。在这些情况下，你可能想要把那些安慰剂检验未通过的变量作为控制变量添加到模型中。如果你在使用局部回归，添加控制变量可能会比较棘手，但有一些方法可以做到这一点\(（Calonico \space et \space al.,2019）\)，而且预先打包好的断点回归命令通常也包含了这些方法。

Calonico, Sebastian, Matias D. Cattaneo, Max H. Farrell, and Rocio Titiunik. 2019. “Regression Discontinuity Designs Using Covariates.” Review of Economics and Statistics 101 (3): 442–51.

20.2.4 密度断点检验

要让断点回归奏效，我们需要处于这样一种情况：刚好在临界值两侧的人是被有效随机分配到处理干预组的。如果有人在暗中操纵，这种情况就不会发生！

操纵可能以两种方式发生。首先，负责设定临界值的人（或事物）可能在完全清楚这会导致谁得到处理（干预）的情况下设定临界值。想象一下田径队选拔时的教练，在他儿子以 \(5\) 分 \(35\) 秒的成绩达标后，就”恰巧”决定队里的任何人都必须在 \(5\) 分 \(37\) 秒或更短时间内跑完一英里。

其次，个体自身可能对他们的运行变量有一定的控制。有时候他们有直接的控制 —— 如果我当地的超市在我购买 \(12\) 罐或更多金枪鱼时给我折扣价，那么我可以精确地选择我的运行变量（我购买的金枪鱼罐数）是多少，而且说买 \(12\) 罐的人和买 \(11\) 罐的人具有可比性，从而能识别出折扣的效应，这其实是不太可信的。

更常见的是，他们有间接的控制 —— 如果我要参加那支田径队的选拔，我可以决定放松跑，或者拼命跑以跑得更快。又或者，计时的人真的很喜欢我，决定在给我记录的时间上减掉几秒。如果我们有间接控制，同时又知道临界值是多少，这就成了问题。如果我知道临界值是 \(5:37\)，而我通常跑 \(5:40\)，我知道我必须格外努力才有机会进入队伍，所以会非常拼命地跑，而通常跑 \(4:30\) 的人就会放松跑。又或者，计时的人可能看到我实际跑了 \(5:38\)，觉得我离达标太近了，稍微作弊一下，把时间写成 \(5:36\)，这样我就能进入队伍。不过，跑 \(6:30\) 的人不太可能被砍掉整整一分钟；那样要求太多了。

在所有这些情况中，我们都希望确保自己了解运行变量是如何被确定的。我们是处于”金枪鱼罐”的情况吗？是”暗中选择临界值”的情况吗？还是处于”我们知道临界值，可能会为了它而操纵自己行为”的情况？对于任何这些情况，我们可能会担心，对临界值或临界值附近运行变量的操纵，会使临界值两侧的人变得不可比。如果我友好的计时员为我减掉了几秒，让我进入了队伍，但却没为同样跑了 \(5\) 分 \(38\) 秒但计时员不喜欢的丹尼³⁸这么做，那么处理干预分配就不是随机的，而是基于计时员喜欢谁，而这存在各种各样的后门路径。

在间接控制的情况下，我们确实有一个可以进行的检验，用来检查在临界值处是否似乎发生了操纵。从概念上讲，它非常简单。我们只需查看运行变量在临界值附近的分布。它应该是平滑的，就像我们预期的那样 —— 如果运行变量是在不考虑临界值的情况下被随机分布的。

但如果我们看到的分布在临界值的一侧似乎有一个下降，而那些观测值刚好落在另一侧，那就很像存在操纵。如果你查看田径选拔的跑步时间，发现只有丹尼的时间在 \(5:38\) 到 \(5:45\) 之间，而有十个人的时间刚好比临界值少一两秒，那就看起来非常可疑！那个计时员似乎没安好心。

将这种直觉付诸实践的想法来自\(McCrary（2008）\)，不过从那以后，对这个基本想法有了改进，包括更好的密度估计和检验效能。而且，就像这本书中几乎所有其他内容一样，这些检验为了在某些情况下有更好的表现仍在不断发展（例如，参见 \(Bugni \space \& Canay（2021）\)的研究），下面代码中的操作说明是针对\(Cattaneo、Michael \space \& \space Ma（2020）\)的方法的。

McCrary, Justin. 2008. “Manipulation of the Running Variable in the Regression Discontinuity Design: A Density Test.” Journal of Econometrics 142 (2): 698–714. Bugni, Federico A., and Ivan A. Canay. 2021. “Testing Continuity of a Density via g-Order Statistics in the Regression Discontinuity Design.” Journal of Econometrics 221 (1): 138–59. Cattaneo, Matias D., Michael Jansson, and Xinwei Ma. 2020. “Simple Local Polynomial Density Estimators.” Journal of the American Statistical Association 115 (531): 1449–55.

\(McCrary\) 的方法是这样的：

计算处理变量的密度分布。
允许该密度在临界值处存在断点。
查看在临界值处是否存在显著的断点。
同时，绘制密度分布并观察，看是否看起来存在断点。

一个大的断点不是好迹象。它意味着存在操纵。

我们可以使用预先打包好的命令来编写密度检验的代码。代码仍将使用《政府转移支付与政治支持》的示例数据。

在运行我们的代码之前，我们可以思考一下这里存在操纵的可能性。政府并非仅基于一个特征来分配处理干预，而是基于一大堆事物的总和。政府在确切知道谁符合条件之前就确定了这个过程（这使得政府难以操纵）。整个过程并非完全透明（这使得人们难以操纵），而且事先很难确切知道自己离临界值有多近（这使得人们难以知道是否应该费心去尝试操纵）。

所以，这项特定的研究似乎不太可能显示出操纵的证据。但谁也说不准！让我们看看。这个检验使用了略微扩展的数据版本，包括那些没有在合适的时间接受调查从而被纳入研究的人，以及那些在估计所用带宽之外的人（所以我们自己要限制在带宽范围内）。

library(tidyverse); library(rddensity)

# Limit to the bandwidth ourselves
gt <- causaldata::gov_transfers_density  %>%
   filter(abs(Income_Centered) < .02)
# Estimate the discontinuity
gt %>%
    pull(Income_Centered) %>%
    rddensity(c = 0) %>%
    summary()

## 
## Manipulation testing using local polynomial density estimation.
## 
## Number of obs =       20463
## Model =               unrestricted
## Kernel =              triangular
## BW method =           estimated
## VCE method =          jackknife
## 
## c = 0                 Left of c           Right of c          
## Number of obs         9077                11386               
## Eff. Number of obs    1927                2309                
## Order est. (p)        2                   2                   
## Order bias (q)        3                   3                   
## BW est. (h)           0.004               0.005               
## 
## Method                T                   P > |T|             
## Robust                -0.9238             0.3556

## Warning in summary.CJMrddensity(.): There are repeated observations. Point
## estimates and standard errors have been adjusted. Use option massPoints=FALSE
## to suppress this feature.

## 
## P-values of binomial tests (H0: p=0.5).
## 
## Window Length / 2          <c     >=c    P>|T|
## 0.000                      24      20    0.6516
## 0.000                      50     104    0.0000
## 0.000                      70     123    0.0002
## 0.000                      79     157    0.0000
## 0.000                     109     176    0.0001
## 0.000                     141     195    0.0038
## 0.000                     160     216    0.0045
## 0.000                     184     237    0.0112
## 0.000                     199     254    0.0111
## 0.000                     226     270    0.0534

import pandas as pd
from rddensity import rddensity
from causaldata import gov_transfers

gt = gov_transfers.load_pandas().data
# Limit to the bandwidth ourselves
gt = gt.query('abs(Income_Centered) < .02')
# Estimate the discontinuity
rdd = rddensity(gt['Income_Centered'], c = 0)
print(repr(rdd))

## Manipulation testing using local polynomial density estimation
## Number of obs:                              1948
## Model:                              unrestricted
## Kernel:                               triangular
## BW method:                             estimated
## VCE:                                   jackknife
## 
## c = 0                               Left of c             Right of c
## Number of obs:                           1127                    821
## Eff. number of obs:                       384                    275
## Order est. (p):                             2                      2
## Order bias. (q):                            3                      3
## BW est.                                0.0068                 0.0071
## 
## Method:                                     T                P > |T|
## Robust                                -0.8251                 0.4093
## 
## P-values of binomial tests (H0: p =  [0.5] ).
## 
## Window Length/2                           < c                     >= c                       P>|T|
## 
## 0.0196                                   1108                      813                        0.0

dotext = "
  * If necessary, findit rddensity and install the rddensity package
  causaldata gov_transfers_density.dta, use clear download
  
  * Limit to the bandwidth ourselves
  keep if abs(income_centered) < .02
  * Run the discontinuity check
  rddensity income_centered, c(0)
"
stata(dotext, data.in = NULL, data.out = TRUE)

## . 
## .   * If necessary, findit rddensity and install the rddensity package
## .   causaldata gov_transfers_density.dta, use clear download
## .   
## .   * Limit to the bandwidth ourselves
## .   keep if abs(income_centered) < .02
## (32,086 observations deleted)
## .   * Run the discontinuity check
## .   rddensity income_centered, c(0)
## Computing data-driven bandwidth selectors.
## 
## Point estimates and standard errors have been adjusted for repeated observation
## > s.
## (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 =        2046
## > 3
## -------------------+----------------------          Model         = unrestricte
## > d
##      Number of obs |      9077       11386          BW method     =         com
## > b
## Eff. Number of obs |      1927        2309          Kernel        =   triangula
## > r
##     Order est. (p) |         2           2          VCE method    =    jackknif
## > e
##     Order bias (q) |         3           3
##        BW est. (h) |     0.004       0.005
## 
## Running variable: income_centered.
## ------------------------------------------
##             Method |      T          P>|T|
## -------------------+----------------------
##             Robust |   -0.9238      0.3556
## ------------------------------------------
## 
## P-values of binomial tests. (H0: prob = .5)
## -----------------------------------------------------
##  Window Length / 2 |       <c         >=c |     P>|T|
## -------------------+----------------------+----------
##              0.000 |       11           9 |    0.8238
##              0.000 |       23          17 |    0.4296
##              0.000 |       31          26 |    0.5966
##              0.000 |       48         102 |    0.0000
##              0.000 |       56         118 |    0.0000
##              0.000 |       65         122 |    0.0000
##              0.000 |       76         142 |    0.0000
##              0.000 |       78         156 |    0.0000
##              0.000 |       82         160 |    0.0000
##              0.000 |      107         173 |    0.0001
## -----------------------------------------------------

正如预期的那样，我们发现在临界值处，收入分布没有统计上显著的断点。好极了！

当然，相比看能否在统计上显著地拒绝”无操纵”的假设，我们更感兴趣的是看”无操纵”这一假设是否合理。它能通过直觉检验吗？我们可以通过绘制分箱的运行变量的分布来检验这一点，如 \(图20.15\) 所示。收入分布在临界值处看起来当然没有任何类型的断点。你可以用 \(R\) 或 \(Python\) 中 \(rddensity\) 包的 \(rdplotdnsity\) 来绘制一个与之非常相似的图。在 \(Stata\) 中，它是常规 \(rddensity\) 命令的一个绘图选项。

20.3 专业人士是怎么做的

20.3.1 回归弯折

断点回归关注的是临界值处结果均值的变化。为什么不关注其他方面的变化呢？确实，为什么不呢？从理论上讲，你几乎可以关注任何事物的变化。为什么不关注结果的标准差的变化呢？或者中位数的变化？又或者结果对一大堆预测变量回归的\(R^2\)的变化？³⁹

不过，最常见的”其他事物”是关注结果与运行变量之间关系的斜率变化。在临界值处实施的处理干预并不会使结果\(Y\)本身发生变化/跃升 —— 相反，处理干预改变了结果变量与()运行变量之间关系的强度。由此得到的图形可能会有点像\(图20.16\)—— 在临界值处没有跃升，但肯定有事情在发生。

为什么处理干预会改变这种关系？它可能以两种主要方式出现。一种可能是，处理干预确实改变了这种关系。例如，假设我们进行一项断点回归，以时间为运行变量⁴⁰，以科技公司的就业数量为结果。然后，在某个特定时间（我们的临界值），一项新政策出台，鼓励对科技公司进行更多投资。我们不会立即看到科技就业人数上升（没有跃升），但我们会预期科技就业人数的增长速度会比以前更快（斜率增加）。

回归弯折更常见的应用场景是，处理干预本身存在弯折，就像 \(图20.16\) 那样，但纵轴是”处理水平/比率”而非”结果”。很多政府政策都是这样起作用的。以失业保险为例，它是回归弯折的常见应用场景。我会引用 \(Card \space et \space al.（2015）\)的研究作为众多例子中的一个，选择它主要是因为我下面还要因为另一个原因引用它。他们研究了奥地利的数据，在奥地利，失业保险福利是正常收入的 \(55\%\)，直到达到最大慷慨额度，之后就不再增加。所以，你的正常收入（运行变量）对你能获得的失业保险金额（处理）有积极影响，直到临界值，在这一点上，正常收入对失业保险金的影响变为\(0\)。

Card, David, David S. Lee, Zhuan Pei, and Andrea Weber. 2015. “Inference on Causal Effects in a Generalized Regression Kink Design.” Econometrica 83 (6): 2453–83.) as one example of many, picking it largely because I was already going to cite it for another reason

所以，比如，如果你认为更慷慨的失业救济金会让人们花更长时间找新工作，那么你应该预期，在临界值之前，你的正常收入（运行变量）和找到新工作的时间（结果）之间存在正相关关系，而在临界值之后，这种关系应该是平稳的。这正是他们的发现。

“处理存在弯折”的回归弯折设计的一个例子是在\(Bana、Bedard \space \& \space Rossin-Slater（2020）\)的研究中。他们研究了加利福尼亚州带薪家庭假的影响。带薪家庭假是一项政策，怀孕的母亲可以休一段延长的带薪假期，这样她们就可以，你懂的，生孩子并照顾新生儿，而不用回复电子邮件。在加利福尼亚州，该州支付正常收入的 \(55\%\)，最高不超过最高福利金额。所以，带薪家庭假的补贴（处理）随着正常收入（运行变量）的增加而增加，直到你的收入水平（临界值）达到能获得最高福利金额的程度。在那之后，额外的正常收入不会增加你的家庭假补贴。这种设定与我之前讨论的\(Card \space et \space al.（2015）\)的失业保险案例非常相似，巧合的是，报销比例也是 \(55\%\)。

Bana, Sarah H., Kelly Bedard, and Maya Rossin-Slater. 2020. “The Impacts of Paid Family Leave Benefits: Regression Kink Evidence from California Administrative Data.” Journal of Policy Analysis and Management 39 (4): 888–929. Card, David, David S. Lee, Zhuan Pei, and Andrea Weber. 2015. “Inference on Causal Effects in a Generalized Regression Kink Design.” Econometrica 83 (6): 2453–83.

首先要检查的是，他们是否确实在临界值处看到了处理的斜率变化，而他们肯定看到了。\(图20.17\)显示了样本中基于计算福利所依据的基期收入而获得的平均带薪家庭假补贴。即使没有在上面加上回归线，在我看来，这也像是斜率发生了变化。

接下来，我们可以将纵轴上的处理干预替换为一些结果，以查看这些结果是否也会发生斜率变化。如果它们发生了变化，那就是存在效应的证据⁴¹。他们在论文中研究了多个结果，但我特别挑选了两个：母亲休家庭假的时长，以及在重返工作岗位的前提下，未来三年内是否再次使用家庭假。我挑选这两个是因为你可以看到其中一个有效应，而另一个没有效应，如\(图20.18\)所示。

在 \(图20.18\) 中我们能看到什么呢？左边我们基本上什么也没看到。在临界值处斜率完全没有变化！这告诉我们，额外的家庭假补贴并没有影响家庭假最终的时长 —— 如果有影响的话，我们会预期当补贴停止增加时，斜率会发生变化。相反，在右边，临界值处的斜率有很大的变化。这告诉我们，更高的休假补贴确实增加了再次申请家庭假的可能性（也就是，又生了一个孩子），因为在福利增加的时候，这样做的可能性在增加，但一旦福利停止增加，这种情况就停止了。

不过，右边的图不只是斜率变化了，它还跃升了！我们该如何理解这一点呢？哦，不，这可能是一个需要关注的问题。处理干预并没有跃升，它只是斜率变化了。所以如果结果跃升了，而不只是斜率变化，这意味着某个假设可能是错误的。

或者，至少如果原始论文确实发现了这样的情况的话会是如此。在原始论文中，根本没有跃升，只有斜率的良好变化。我只是在两侧画了普通的直线，但他们做了一些更精巧的处理。⁴²无论如何，这个跃升的幅度还不足以令人信服地说明它是非零的。我自己把它加进来，是为了有理由谈论这个问题，并强调对结果进行直觉检验的重要性。

一旦你掌握了断点回归的概念，实际进行回归弯折就不是太难了。你只是在那基础上增加了一层复杂性。我已经讲过的大部分内容在这里仍然适用 —— 核函数、带宽、多项式、局部回归等等。同样的安慰剂检验在很大程度上也适用 —— 例如，如果你用控制变量作为结果来估计回归弯折模型，斜率不应该有变化，而且你要确保运行变量在临界值处不是不连续的。代码变化不大 —— 只需在你的 \(rdrobust\) 代码中把 \(deriv\)（你想要回归方程的几阶导数）选项设为 \(1\) 即可。

甚至基础的回归方程也没有太大变化。核心思想是，你仍然只是在临界值的左右两侧估计一条线，就像回到 \(方程（20.1）\) 时那样：

\(\begin{multline} Y = \beta_0 + \beta_1 (\text{Running-Cutoff}) + \beta_2 \text{Treated} + \beta_3 (\text{Running-Cutoff}) \times \text{Treated} + \varepsilon \end{multline}\)

方程有时会有差异：一些研究者会允许斜率变化，但不允许有任何类型的跃升，也就是说，从模型中去掉\(\text{Treated}\)本身，设定\(\beta_2 = 0\)。去掉\(\text{Treated}\)有助于减少临界值附近噪声的影响（而回归弯折可以利用这种降噪的帮助），不过如果那里实际上应该有一个跃升，你就不会注意到它。⁴³另一个差异是，现在，你不再关注代表跃升的系数\(\beta_2\)（如果它还在模型中的话），而是关注代表斜率变化的系数\(\beta_3\)。

当然，如果你使用的设定中，因为处理水平存在弯折而出现了弯折，你会像之前在模糊断点回归中那样，将临界值用作处理的工具变量⁴⁴。

那么还有其他不同的地方吗？当然有一些。首先，最优带宽选择程序是不同的。另外，你也更有可能使用均匀核而不是三角核（这很方便地让计算变得容易得多 —— 只需去掉带宽之外的观测值然后进行分析即可！）\((Card \space et \space al.，2017)\)。查看 \(rdrobust\) 中 \(kernel\) 的 \(uniform\) 选项。而且，正如你可能预期的那样，由于我们现在关注的是斜率而不仅仅是结果的跃升，未能恰当地对非线性关系进行建模，在这里比在断点回归中更成问题\(（Ando, 2017）\)。

Card, David, David S. 2017. “Regression Kink Design: Theory and Practice.” In Advances in Econometrics, 38:341–82. Bingley: Emerald Publishing Limited. Ando, Michihito. 2017. “How Much Should We Trust Regression-Kink-Design Estimates?” Empirical Economics 53 (3): 1287–1322.

需要注意的一个大差异是，统计功效问题 —— 断点回归中本来就存在的这个问题 —— 变得更加紧迫。精确估计关系的变化（斜率）比估计均值的跃升需要多得多的观测值。你试图仅使用临界值附近的观测值来估计那个斜率。你最好希望在那个临界值附近有大量的观测值。

20.3.2 临界值的划分

明确的临界值、模s糊的临界值、改变斜率的临界值。谁能应对所有这些临界值呢？我希望是你，因为你马上就要被临界值淹没了。当然，存在多临界值设计这种情况。

要是这还不够的话，实际上还有几种不同类型的多临界值设计。多重多临界值设计。这真让人疲惫。至少让我们简要地讨论一下它们吧。

有时候，临界值是一个移动的目标。给定的处理（干预）可能会有不同的临界值，这取决于它所处的位置、针对的对象…… 你能想到的情况都有可能！这相当常见。例如，大学招生办公室可能对大多数学生有一个入学考试的临界分数，而对申请校队职位的学生有一个不同的临界分数。或者，就考试而言，在美国，如果你想要一个 \(GED\)（高中文凭的替代证书），你必须参加考试并达到一定的最低分数。但每个州的临界分数是不同的。又或者，你可能想用断点回归来分析第三方势力盛行的选举。在一些选举中，你可能需要获得 \(50.1\%\) 的选票才能获胜，但在另一些选举中，你可能只需获得比如 \(42.7\%\) 的选票就能获胜。

又或者，可能只有一个临界值，但它会随着时间而变化。上一节中关于带薪家庭假的\(Bana、Bedard \space \& \space Rossin-Slater（2020）\)的论文就是这种情况的一个例子。你能拿到最高家庭假补贴的季度收入随着时间的推移而上升，从\(2005\)年略低于\(20000\)美元到\(2014\)年略高于\(25000\)美元。

所以，临界值会按群体、跨地区或跨时间发生变化，但想必我们希望把所有这些数据整合在一起 —— 毕竟，在统计功效低的断点回归中，观测值是很宝贵的。我们能做些什么呢？

最常见的方法就是简单地对数据进行中心化处理，并假装没有什么不同。就像你通常做的那样，使用\(运行变量-临界值\)来估计你的模型，只是允许每个人的\(运行变量-临界值\)根据他们的情况而有所不同。这是个简单的解决办法。\(Bana、Bedard \space \& \space Rossin-Slater\)正是这么做的。

但这是我们唯一能做的吗？不是的。对数据进行中心化处理是有效的，但它实际作用的对象略有不同，尤其是当你考虑处理效应的平均值时\(（Cattanel \space et \space al.,2016）\)。常规的断点回归给我们的是局部平均处理效应 —— 也就是刚好在临界值附近的那些人的平均处理效应。那么，如果我们有多个临界值，会得到什么呢？它并非恰好是那些局部效应的平均值。相反，你会得到那些效应的加权平均值，权重是接近各个临界值的观测值数量。这比仅仅 “临界值处那些的平均值” 要复杂一些，这可能会影响你的结果对于政策而言实际有多大的用处。

Cattaneo, Matias D., Rocío Titiunik, Gonzalo Vazquez-Bare, and Luke Keele. 2016. “Interpreting Regression Discontinuity Designs with Multiple Cutoffs.” The Journal of Politics 78 (4): 1229–48.

此外，把所有东西混在一起，你就丢弃了一些有用的信息。既然你可以利用不同的临界值来尝试获取远离临界值的处理效应，为什么要满足于某个奇怪的平均值呢？或许你有临界值低的 \(A\) 组和临界值高的 \(B\) 组。在对处理的外推方式做一些假设的情况下，我们可以对 \(A\) 组中远离 \(A\) 组临界值、但靠近 \(B\) 组临界值的人估计处理效应\(（Cattaneo \space et al.,2020）\)。关于同一想法的另一种方法，参见 \(Bertranha（2020）\)。⁴⁵

Cattaneo, Matias D., Luke Keele, Rocío Titiunik, and Gonzalo Vazquez-Bare. 2020. “Extrapolating Treatment Effects in Multi-Cutoff Regression Discontinuity Designs.” Journal of the American Statistical Association, 1–12. Bertanha, Marinho. 2020. “Regression Discontinuity Design with Many Thresholds.” Journal of Econometrics 218 (1): 216–41.

现在，多个临界值不再给我们一个像土豆泥一样的处理效应平均值，而是让我们得到比简单直接的单临界值版本更具信息量的估计值。用这个类比的话，这是炸薯条吗？我思路断了。

方便的是，这些针对多临界值的方法与 \(rdrobust\)、\(rddensity\)、\(rdpower\)一起，被纳入了断点回归软件包家族，在三种语言中都以 \(rdmulti\) 的形式存在。\(rdmulti\)包包含\(rdmc\)函数，可用于评估多临界值设置⁴⁶。

\(rdmulti\)包还提供了\(rdms\)函数，可用于另一种存在多临界值的情况 —— 当所有人都受相同的临界值约束，但同一运行变量的不同临界值会带来不同种类或水平的处理（干预）。例如，考虑血压药物治疗，如果你血压上升到”潜在担忧”范围，会给你低剂量；如果血压变得危险地高，会给你高剂量。这里的直觉与在每个临界值处分别进行一次断点回归的想法并没有太大不同。

当然，那都是在你的多个临界值仅来自一个运行变量的情况下。多没劲啊！我们这些”酷小孩”可是有两个运行变量的。在相当多的情况下，处理（干预）的分配不是基于一个运行变量的临界值，而是基于两个（或更多）运行变量各自的临界值。

这种情况在地理断点回归中很常见。在基于地理的断点回归中，处理（干预）发生在边界的某一侧（就像本章开头我们举的圣地亚哥/蒂华纳的例子）。在地理断点回归的许多应用中，要从一个未处理区域进入一个处理区域，你必须同时跨越某一纬度和某一经度。这就有两个运行变量！⁴⁷

虽然多个运行变量的应用随处可见，但在教育领域这种情况似乎相当常见\(（Reardon \space \& \space Robinon, 2012）\)。例如，或许要升入下一年级，学生需要在数学考试和英语考试中都达到最低分数。或者，要获得奖金，教师需要他们所教的两个班级的学生都达到某一标准。又或者，一所学校只有在多项表现阈值上都低于标准时，才会被认定为未达到州标准。

Reardon, Sean F., and Joseph P. Robinson. 2012. “Regression Discontinuity Designs with Multiple Rating-Score Variables.” Journal of Research on Educational Effectiveness 5 (1): 83–104.

这是个问题，因为现在实际上我们有了无限多的临界值。或者至少是多维的临界值。假设你需要数学和英语都考 \(60\) 分才能升入下一年级。这不仅仅是数学有个 \(60\) 分的临界值、英语有个 \(60\) 分的临界值。当英语是 \(61\) 分时，数学有个 \(60\) 分的临界值；当英语是 \(61.1\) 分时，数学有个 \(60\) 分的临界值；当英语是 \(73.7\) 分时，数学有个 \(60\) 分的临界值，还有……

\(图20.19\) 说明了这个问题。对于单一运行变量，很容易弄清楚你在比较谁 —— 那些在临界值附近带宽内的人。但对于两个运行变量，你仍然想要比较临界值附近的人，但你要把他们和谁比较呢？临界值是\(L\)形的！你可能不想把\(L\)形左上角的人与右下角的人进行比较 —— 他们可能并不那么相似。而且，如果\(L\)形的某些部分比其他部分人口更多，直接进行比较可能无法得到你想要的结果。所以你必须想办法进行恰当的比较。

不仅这些群体可能不具有可比性，而且处理（干预）的效应对于\(L\)形的所有部分可能也不一样。留级对数学不及格的人和英语不及格的人可能有非常不同的影响。所以你必须决定要估计\(L\)形的哪部分（或哪些部分）的效应。

在这些情况下，有几种估计方法，\(Reardon \space \& \space Robinon（2012）\)介绍了其中五种。也有一些特定版本，比如 \(Papay、Willett \space \& \space Murnane（2011）\)中讨论的那种，专门针对你所获得的处理（干预）因跨越 \(L\) 形的不同部分而不同的情况 —— 也许不是留级，数学不及格的孩子参加数学补习，而英语不及格的孩子参加英语补习。

Reardon, Sean F., and Joseph P. Robinson. 2012. “Regression Discontinuity Designs with Multiple Rating-Score Variables.” Journal of Research on Educational Effectiveness 5 (1): 83–104. Papay, John P., John B. Willett, and Richard J. Murnane. 2011. “Extending the Regression-Discontinuity Approach to Multiple Assignment Variables.” Journal of Econometrics 161 (2): 203–7.

一种可能的方法，也是\(rdmulti\)包中\(rdms\)函数采用的方法，是在\(L\)形上选取特定的点，然后使用另一个运行变量来估计相当常规的断点回归。在\(L\)形上选取足够多的点进行这样的操作，你就能了解效应是什么，以及它如何根据跨越边界的位置而变化。

20.3.3 当运行变量不安分时

断点回归很不错。它就像一场度假。在临界值附近的区域，因果图中指向处理变量的箭头很少甚至没有。事情很简单，结果也令人信服。我们悠闲地躺着，用椰子喝着伏特加，漫不经心地想着是用线性模型还是二阶多项式。

然后雨云来了。假期毁了。运行变量不”安分”了。

在断点回归中，运行变量除了要负责分配处理（干预）外，最主要的是要良好、连续且平滑，就像海滩上轻轻退去的潮水。

我们已经讨论过运行变量被操纵可能带来的问题。但还有另外两个潜在问题：运行变量的粒度太粗—— 测量的层级很粗略，或者运行变量出现堆点现象—— 有些值可疑地比其他值常见得多。这两种情况都可能引发问题。

粒度指的是变量测量的层级很粗略。例如，如果我们要测量年收入，我们可以说你去年挣了 \(40231.36\) 美元，或者说你挣了 \(40231\) 美元，或者说你挣了 \(40000\) 美元。或者我们甚至可以说你挣了”\(4\)万到\(5\)万美元”。或者，嘿，我们可以说”不到\(10\)万美元”。这些是对收入的测量，粒度依次降低。

你可以想象这对断点回归来说为什么可能是个问题。断点回归的核心思想是，临界值两侧的人应该很容易进行比较，处理（干预）是他们之间唯一真正的差异。如果临界值是\(40000\)美元，我们把挣 \(40231.36\) 美元的你和挣 \(39990.32\) 美元的人进行比较，这可能相当可信。但如果我们只把收入测度为 \(4\) 万到 \(5\) 万美元呢？那我们不仅把你和挣 \(39990.32\) 美元的人比较，还会把挣 \(49999\) 美元的人和挣 \(30001\) 美元的人进行比较。要说处理（干预）是区分他们的唯一因素，就不那么容易了。

多粗的粒度才算足够细呢？当然，这是主观的，但只要你能足够精确地测量变量，从而区分出你认为具有可比性的人和不具有可比性的人，那就可以了。

话虽如此，即使你能足够好地区分从而做到这一点，但如果变量的测量仍然不够精细，你可能无法恰当地对结果与运行变量在接近临界值处的关系形状进行建模。如果结果随着收入从 \(30000\) 美元到 \(40000\) 美元而增加，但你只有 \(30000\) 到 \(40000\) 美元这个区间，那么你会错误地预测临界值处的值，这样你的估计也会出错。至少，你无法对那种趋势进行建模会增加你的抽样变异性。

这里真正的解决办法是找到一个粒度足够细的运行变量。但没有变量的粒度是无限细的。所以，如果你担心粒度问题，你可能想要选择一个会考虑到该粒度的估计量。

有几篇源于 \(Kolesar \space \& \space Rothe（2018）\)的论文做了这件事。\(Kolesar \space \& \space Rothe\) 通过对结果与运行变量在我们无法深入探究的那个大而粗的\(30000\)到\(40000\)美元区间内的关系做出一些假设，从而更好地估计断点回归估计量的抽样变异性。他们的其中一个估计量假设关系不是太非线性，另一个估计量假设你对形状的判断不会比在临界值处的判断更错。这些估计量在 \(R\) 语言的 \(RDHonest\) 包中可以找到。

\(Kolesar \space \& \space Rothe（2018）\)指出的一点是，你不应该试图通过对运行变量的不同值进行标准误聚类来解决这个问题。长期以来，这是一种标准做法，但事实证明，在很多情况下，这比完全不解决这个问题还要糟糕。哎呀！别这么做。

Kolesár, Michal, and Christoph Rothe. 2018. “Inference in Regression Discontinuity Designs with a Discrete Running Variable.” American Economic Review 108 (8): 2277–2304.

非随机堆点是指运行变量更有可能取某些值而非其他值的情况。这种情况通常以四舍五入的形式出现。问人们多大年纪，绝大多数人会给你一个整数，比如36 岁。不过，有些人会更精确，说”\(36\) 岁 \(8\) 个月零 \(2\) 天”。典型的例子来自 \(Almond、Doyle \space \& \space Williams（2010）\)的研究，他们发现婴儿出生体重极有可能被报告为四舍五入到最近的 \(100\) 克。

这些过程中的任何一个都可能导致运行变量的分布看起来像\(图20.20（a）\)中的模拟数据。你可以看到，每当运行变量达到某个特定值（这里是一个整数）时，取该值的观测比例就会大幅跃升。这些就是”堆点”。

现在，这本身可能不是什么大问题。在\(X\)轴的某些点上有大量观测值会严重影响你拟合的回归线。但平均而言，没什么坏处，也没什么错。真正的问题出现在堆点是非随机的时候。注意\(图20.20（b）\)，那里不仅堆点很常见，而且它们也有所不同。堆点是未填充的圆圈，它们的平均结果值比其他的更高。

只要存在堆点，就可能出现非随机堆点。毕竟，四舍五入回答的人和不四舍五入回答的人（或医院）真的是同一类吗？

对于非随机堆点，我们有运行变量中极具影响力的值（因为它们比其他值常见得多），但这些值也不具代表性。这会使我们的估计偏离真实效应，甚至可能改变其符号！你可以想象，如果临界值正好在其中一个堆点上，情况会特别糟糕 —— 临界值的一侧是堆点，另一侧不是。当堆点本身就存在差异时，这就会导致偏差。

解决这个问题的一种常见方法是”甜甜圈孔断点回归”，也就是简单地去掉临界值附近的观测值，以清除临界值附近的堆点。在断点回归中去掉近临界值观测值似乎有违直觉，因为断点回归正是非常需要那些近临界值观测值的。但这就是我们面临的问题。

不过，\(Barreca、Lindo \space \& \space Waddell（2016）\)研究了堆点问题和甜甜圈孔解法。他们表明，即使堆点远离临界值，它也是个问题。所以，除非你消除了唯一的堆点，否则甜甜圈孔方法无法解决问题。相反，他们建议根据一个值是否为堆点来拆分样本，并分别分析每一侧。这消除了处于堆点中和结果之间存在明显后门路径的问题。

Barreca, Alan I., Jason M. Lindo, and Glen R. Waddell. 2016. “Heaping-Induced Bias in Regression-Discontinuity Designs.” Economic Inquiry 54 (1): 268–93.

20.3.4 处理带宽问题

离临界值多远时，你仍然能在其两侧获得具有可比性的观测值？ 这是一个重要的问题，并且存在真正的权衡。如果选择一个围绕临界值过宽的带宽，你会纳入不具可比性的观测值（并且更多地依赖你选择了正确函数形式的想法），这会让你的估计更不可信且更有偏差。如果选择的带宽过窄，你最终会在几乎没有数据的情况下估计效应，导致估计值存在噪声。这就是”效率 - 稳健性权衡”⁴⁸。你想要一个基于大量观测值的精确有效估计吗？那么，你就会引入一些偏差，这会使你的估计在存在后门路径或函数形式选择不佳的情况下，稳健性更差。

带宽应该是多少呢？你可以采用三种主要方法：

直接选一个带宽。这是一种相当常见的方法，尽管随着时间的推移，它可能变得不那么常见了。对于这种方法，你只需查看数据，然后思考在不觉得越界的情况下，能延伸多远。理想情况下，你能证明为什么这个带宽是合理的。或者，你可能只是使用你拥有的所有数据，并认为你所拥有的范围已经足够好了。
选一堆带宽。这种基于敏感性测试的方法也相当常见 —— 本章前面提到的 \(Fetter（2013）\)以及\(Bana、Bedard \space \& \space Rossin-Slater（2020）\) 的论文都采用了这种方法（此外还有下面的一点选项 \(3\) 的内容）。对于这种方法，你选择你认为对研究设计有意义的最大带宽，然后开始缩小它。每次都要估计你的断点回归模型。

一旦你从每个带宽中得到了估计值，你就会呈现所有结果，有可能用 \(图20.21\) 这样的图表，我用 \(Fetter（2013）\)的数据，通过一个非常基础的非局部线性断点回归，改变带宽大小制作了这个图。我们可以在图上看到带宽过小带来的噪声风险。实际上，你看不到…… 我不得不把最小的带宽（\(1-3\) 个月）从图中删掉，因为它们的估计值太不稳定了，置信区间也太大，以至于你根本看不到其他估计值。你能看到的是，随着带宽扩大，置信区间缩小。这很合理 —— 更大的带宽意味着更多的观测值，进而意味着更小的标准误。

\(图20.21\) 大致是你使用这种方法时想要看到的。在合理范围内，当我使用越来越宽的带宽时，估计值似乎没有太大变化。所以我们选择哪一个可能没那么重要（在原文中，\(12\) 是主要分析中的带宽）。

让数据为你选择带宽。或者更准确地说，确定”一个好的带宽是什么样的”这一目标，然后利用你的数据，根据该标准找出最适合的带宽宽度。

有几种方法可以做到这一点。第一种，在\(Imbens \space \& \space Lemieux（2008）\)中讨论过，是交叉验证。我在第 \(13\) 章介绍过交叉验证。基本思路是，你试图获得最佳的样本外预测能力。所以你把数据随机分成若干块。对于每一块，你用所有其他块来估计你的模型，然后对现在处于样本外的那一块（你留下的那块）进行预测。对每一块都重复这个过程，然后看看总体上你的样本外预测效果如何。接着，对每个潜在的带宽重复整个过程，看看哪一个效果最好。

Imbens, Guido W., and Thomas Lemieux. 2008. “Regression Discontinuity Designs: A Guide to Practice.” Journal of Econometrics 142 (2): 615–35.

交叉验证方法曾经非常流行，现在你也仍然经常能看到它。\(Fetter（2013）\)的论文就是这么做的。然而，我们真正想要的是在临界值处有良好的预测。标准的交叉验证方法可能会过于关注拟合远离临界值的点。

第二种方法是”最优带宽选择”方法，源于\(Imbens \space \& \space Kalyanaraman（2012）\)。这种方法的目标是在临界值处（就在其两侧）获得最佳预测。当然，我们不知道临界值处的真实值是多少，所以我们会尽最大努力去估计它。

Imbens, Guido W., and Karthik Kalyanaraman. 2012. “Optimal Bandwidth Choice for the Regression Discontinuity Estimator.” The Review of Economic Studies 79 (3): 933–59.

\(Cataneo \space \& \space Cazquez-Bare（2016）\)描述了几种不同的选择最优带宽方法的流程，就这类论文而言，它不算太难读。简而言之，通过追求在临界值处获得最佳预测这一目标，最终你会希望，对于更大的样本，带宽更小（因为你能承受损失观测值）；使用的多项式项越多，带宽越大（因为你可以通过扩大范围更好地拟合你引入的任何非线性）；除此之外，还基于错误函数形式所引入的偏差程度，以及许多其他细节。

Cattaneo, Matias D., and Gonzalo Vazquez-Bare. 2016. “The Choice of Neighborhood in Regression Discontinuity Designs.” Observational Studies 3 (2): 134–46.

\(Cataneo \space \& \space Cazquez-Bare（2016）\)的论文还指出，使用最优带宽选择方法会搞乱标准误和置信区间。解决这个问题的办法是应用偏差校正，估计错误函数形式所带来的偏差有多大，并据此调整效应的估计值和方差的估计值。

他们进一步深入，讨论了带宽选择程序，这些程序不是试图在临界值处获得最佳预测，而是试图最小化覆盖误差—— 即估计的抽样分布与实际抽样分布之间的不匹配。这个过程往往会得到更窄的带宽，而在撰写本文时，运行 \(rdrobust\) 包中的 \(rdrobust\) 函数就会得到这样的结果。

这些带宽选择程序相当重要。再以 \(Fetter（2013）\)为例，原文中交叉验证方法导致选择了 \(12\) 个季度的带宽，用我们的方法得到的估计值是 \(0.17\)（用原文中的方法是 \(0.177\)）。我们也看到，\(rdrobust\) 方法选择了 \(3.39\) 个季度的带宽，得到的估计值是 \(0.518\)—— 非常不同。在这种情况下，确实是带宽选择造成了差异。强制\(rdrobust\)使用 \(12\) 个季度的带宽，会使估计值达到 \(0.2\)，非常接近原文的结果。

在圣地亚哥 / 蒂华纳的例子中，实际上涉及多种处理，因为在边境有很多事情都发生了变化。但边境显然起到了某种作用！↩︎
断点回归几乎总是用于涉及二元处理变量的情况。本章不会费心去调整措辞以暗示任何其他用法，尽管存在涉及非二元处理的断点回归方法，而且你会在下面 “专业人士怎么做” 的 “回归弯折” 部分看到一点相关迹象。↩︎
在此示例及许多其他示例中，若高于断点则接受处理，低于则否。但在许多场景中，若低于断点亦可接受处理。其逻辑原理相同。↩︎
或者至少您的处理概率会发生巨大变化——我们将在讨论模糊断点回归时涉及这种变体。↩︎
更准确地说，我们需要结果变量与运行变量之间的关系具有连续性。例如，如果结果变量与运行变量之间存在正相关关系，我们预期断点右侧的结果会高于左侧。这完全没问题，即使两组并不完全可比。只要在没有处理的情况下，这种正相关关系不会出现断点。大致上，您可以这样理解：“在调整了结果变量与运行变量之间的整体关系后，两侧基本上是随机分配的。”↩︎
正如本书多次提到的，此图表并非严格意义上的因果图。箭头允许任何形式的关系，因此”运行变量仅在断点处影响处理”可以通过”运行变量→处理”箭头来表示，同时将断点纳入”限制图”中——即将”运行变量→断点”本身作为一个节点。但我认为这种方式比我的做法更复杂，故采用当前形式。如果您持续使用因果图，更规范的版本将会出现。↩︎
您可能已注意到图20.1与第17章事件研究中的因果图相似。实际上它们几乎相同，只是事件研究以时间作为运行变量。本章您还会发现与事件研究（尤其是中断时间序列分析）的其他相似之处。主要区别在于：使用非时间运行变量时，许多必要假设会变得更合理，其原因希望后文能得以阐明。↩︎
尽管实际估计通常比此处所述更复杂——我们将在”操作实施”章节详细讨论。↩︎
见鬼，我们本可以凭直觉猜测这些数值。不过学术期刊往往不认可这种做法。那群象牙塔里的势利眼全都一个样。↩︎
带宽选择的基本矛盾在于：带宽范围越大，就越难以确保”断点两侧个体可比”这一识别假设的可靠性；但带宽范围越小，可用观测值就越少，从而导致估计精度下降。↩︎
请注意这些图中\(y\)轴表示的是接受处理的观测值比例，而非平均结果。↩︎
发现这种消费下降对经济学家来说是个谜题，因为标准模型预测人们通常不愿经历消费的突然大幅下降，而退休这种情况本应能够通过提前调整消费来避免这种下降（如果他们愿意的话）。这种事乍看像是经济学家做出的非常愚蠢且不现实的预测，但越想越觉得有道理。↩︎
我们将在”操作方法”章节深入探讨更多细节。↩︎
这是一个关于反事实的假设——我们无法真正验证它。我们只需要深入理解研究背景，并仔细思考在断点处可能发生的其他变化。↩︎
在修正”趋势”之后。↩︎
存在一些理解断点回归的方法，并不真正要求”在断点附近随机”。但本章将重点保持这种最简单的理解方式。↩︎
进行断点回归前最好先进行功效分析（如第\(15\)章所述），以确认是否有足够观测值进行精确估计。在\(R\)、\(Stata\)和\(Python\)中，rdpower包专门用于断点回归的功效分析。在\(R\)和\(Python\)中可正常安装该包，但\(Stata\)中存在另一个功能不同的rdpower包，因此需访问https://rdpackages.github.io/rdpower/获取安装指引。↩︎
当我们开始估计模糊断点回归时，它们将不再相同。↩︎
在这里，具体如何实现异方差稳健标准误确实变得有些复杂。您可能会看到一些论文使用在运行变量上聚类的标准误——我们过去认为这是个好主意，但现在不再这么认为了。我将在”专业人士如何操作”部分讨论这个问题。↩︎
话虽如此，有时添加控制变量可以通过减少未解释变异来提高估计的精确度。不过，如果您正在进行局部回归，添加控制变量的过程可能会有些复杂，正如我即将描述的那样。但确实存在相关方法，例如Calonico等人（2019）提出的方法，这些方法已在预打包的断点回归命令中实现。

Calonico, Sebastian, Matias D. Cattaneo, Max H. Farrell, and Rocio Titiunik. 2019. “Regression Discontinuity Designs Using Covariates.” Review of Economics and Statistics 101 (3): 442–51.↩︎
您必须将运行变量中心化处理才能使此结论成立。如果忘记对运行变量进行中心化处理，虽然仍能从模型中得出断点回归估计值，但难度会大幅增加。因此请务必将运行变量中心化！↩︎
局部回归也被称为”非参数回归”的一种形式，因为它不对整个数据范围强制设定特定形状。但我需要指出，它确实对数据子集强制设定了形状。因此在我看来，它更像是”较少参数回归”。↩︎
灵活性和对断点的关注是优势，但就像几乎任何丢弃数据或放松假设的情况一样，通过对数据施加较少结构，您会获得较低精度的估计。如果您需要施加的假设非常错误，这是一个很好的权衡！但如果您想使用局部回归，最好在断点附近有大量观测值，否则估计将不够精确而无法使用。↩︎
社会紧急国家援助计划（Plan de Atención Nacional a la Emergencia Social，简称PANES）。这些支付金额也很可观——每月约70美元，相当于典型受益者此前收入的一半左右。此外，任何有子女的家庭还会获得一张食品卡。还有一些其他福利。↩︎
此处的部分结果与原始论文略有不同，这是为了简化数据准备过程或突出一些您可能做到而他们未实施的潜在操作。↩︎
您是否注意到本章中有多少图表？断点回归设计是一种极具视觉魅力的研究方法。↩︎
特别说明：与此书所有图表一样，我使用了ggpubr包中的theme_pubr()主题，并在theme函数中添加了一些额外的文本选项。花时间学习如何使用您选择的绘图包制作_美观_的图表是非常值得的。↩︎
使用原始论文中的特定方法，他们的估计结果更接近11到13个百分点。↩︎
参见 https://rdpackages.github.io/ 及其引用的众多计量经济学论文。↩︎
系数为正，但其假设处理发生在断点右侧，因此这转化为被处理的负面效应——而在此案例中处理实际发生在断点左侧。↩︎
在这个图中，\(Z\)可能代表这里的一堆不同变量 —— 它甚至可能是不同的变量，带有指向运行变量和处理（干预）的箭头。其逻辑仍然成立。↩︎
由于我们进行的那些交互作用，我们所做的事情与基本形式并不完全相同，但逻辑能很好地延续下来。↩︎
有些人是不被允许的。比如，当时他们并没有招募很多女性。但即使在符合条件的男性中，也不是每个人都会选择入伍或者最终被征召入伍。↩︎
但如果只把自己局限在最干净的分析中，会让我们根本无法研究很多事物。↩︎
另一种思考方式是，抓住上下两张图，把它们拉得更高，直到图 20.13（a）中的跃升从 0% 变为 100%。最终在图 20.13（b）中得到的跃升就是模糊断点回归（RDD）的估计值 —— 也就是如果每个人都从未接受处理（干预）变为接受处理（干预）（从 0% 到 100%），你预期会出现的跃升幅度。↩︎
一般来说，对于你在分析中必须选择的任何参数，比如带宽，尝试几个不同的参数并报告结果如何变化，这并非是个糟糕的主意。也许结果不会有太大变化，这会让人安心。又或者结果会有变化，而了解这一点并且诚实地报告出来，也是好的。↩︎
只要它不是处理（干预）本应该影响的事物就行。不过话说回来，我们从第 8 章就知道，反正我们通常也不会想要控制处理（干预）会影响的事物。↩︎
丹尼可能是自作自受。↩︎
这些中的每一个都需要仔细思考究竟该如何对它们进行估计，但从概念上讲，它们应该是可行的。我曾去寻找是否有金融学论文使用过断点回归的标准差版本，结果很惊讶，没找到一篇。这对你们这些金融学专业的学生来说是个机会。（或者可能只是我没找得足够仔细。）↩︎
换句话说，就是第 17 章中的事件研究。这里提到的设计是不明智的，因为它并没有真正满足第 17 章中所有的事件研究假设，但它仍然能把意思传达清楚。↩︎
然后，你可以将在结果中看到的斜率变化除以在处理（干预）中看到的斜率变化，也就是说，进行工具变量分析，以得到处理（干预）的效应。↩︎
他们的模型实际上并不允许存在跃升，而是将其排除了。但他们在估计中纳入了其他几个细节，这使得那个假设是合理的。↩︎
你可以先允许存在跃升进行一次估计，看看它是否存在，如果没什么大问题，再在不考虑跃升的情况下进行估计以提高精度。↩︎
不过现在还有一个额外的 “模糊性” 层面需要考虑 —— 如果处理（干预）本身的弯折没有被完美实施呢？在上面引用的 Card 等人（2015）的研究中有关于这一点的说明。我就说过会在这里引用它的。Card, David, David S. Lee, Zhuan Pei, and Andrea Weber. 2015. “Inference on Causal Effects in a Generalized Regression Kink Design.” Econometrica 83 (6): 2453–83.↩︎
注意到这个叫 “Cattaneo” 的家伙是怎么一直出现的吗？↩︎
斯科特・坎宁安的说唱主题因果推断教材《因果推断：混音带》确实介绍了rdrobust和rddensity包的使用，但跳过了rdmulti及其rdmc函数，错过了一个绝佳的 Run-DMC（美国著名说唱组合）典故的机会。↩︎
基尔和蒂图尼克（2015）讨论了地理断点回归的双运行变量性质（这本身就可以作为本章的一整节内容），以及地理断点回归的其他一些特征。这些特征包括运行变量存在操纵的可能性非常高（人们往往会根据选择在边界的某一侧或另一侧定居），还有跨越边界几乎总是会带来多种处理（干预），而不只是一种（从圣地亚哥到墨西哥的跨越所产生的影响，不能只归因于从美元到比索的转换；那里还有很多其他情况在发生）。

Keele, Luke J., and Rocio Titiunik. 2015. “Geographic Boundaries as Regression Discontinuities.” Political Analysis 23 (1): 127–55.↩︎
令人惊讶的是，这只是我们在其他章节中讨论过的老掉牙的偏差 - 方差权衡在断点回归中的另一种说法。↩︎

第20章 断点回归