Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
73,813 characters · 37 sections · 42 citation commands
Inference in Regression Discontinuity Designs with Clustered Data This version: . We thank Debopam Bhattacharya, Morten Nielsen, Zhuan Pai and numerous seminar and conference participants for helpful comments and suggestions. The second author gratefully acknowledges financial support from the European Research Council ERC through grant SH-1852332. Author contact information: Claudia Noack, Department of Economics, University of Bonn. E-Mail: [email removed]. Website: https://claudianoack.github.io. Tomasz Olma, Department of Statistics, Ludwig Maximilian University of Munich. E-Mail: [email removed]. Website: https://tomaszolma.github.io. Christoph Rothe, Department of Economics, University of Mannheim. E-Mail: [email removed]. Website: http://www.christophrothe.net.
\pagestyle{plain} \onehalfspacing
\onehalfspacing
Regression discontinuity (RD) designs are widely applied in economics and other social sciences. In these settings, treatment is assigned whenever a unit's realization of the running variable crosses a known cutoff; for instance, a candidate is elected only if their vote share exceeds 50%. Under continuity conditions on the conditional expectations of the potential outcomes, the average treatment effect at the cutoff is identified by the jump in the conditional expectation of the observed outcome given the running variable at the cutoff. This jump is typically estimated as the difference of the local linear estimates on either side of the cutoff.
Most theoretical results on estimation and inference in RD designs are derived in settings where the researcher observes a sample of independent and identically distributed (i.i.d.) observations drawn from a large population hahn2001identification, imbens2012optimal, calonico2014robust, armstrong2020simple. In practice, however, applied researchers often regard the i.i.d.\ assumption as unrealistic and routinely report clustered standard errors. For example, we revisited the survey of recent RD studies published in the journals of the American Economic Association conducted by noack2021flexible and found that clustered standard errors were used in around 80% of the surveyed articles. Despite the prevalence of clustered sampling in applied RD studies, formal results for local linear RD estimators -- and more generally local polynomial regression -- under such dependence structures remain limited. In particular, the existing literature offers little guidance on the conditions under which these estimators are asymptotically normal across the range of clustering patterns encountered in practice, or on how to conduct valid inference in such settings.
This paper makes two main contributions. First, we provide an asymptotic theory for the local linear RD estimator for clustered data with a large number of independent groups. From a statistical perspective, RD estimators are weighted averages of the outcome variable where the weights depend on the running variable, the kernel, and the bandwidth. When units are clustered, the interaction between local weighting and within-cluster dependence alters the asymptotic behavior of the local linear RD estimator. We derive high-level conditions under which the local linear RD estimator is asymptotically normally distributed. These conditions are formulated in terms of the weights assigned to units from different clusters and translate into restrictions on cluster sizes within the estimation window. To relate our high-level conditions to empirically relevant designs, we introduce four stylized asymptotic frameworks motivated by empirical RD applications. These frameworks capture how the asymptotic behavior of local linear RD estimators depends on (i) the effective number of units per cluster within the estimation window, (ii) the dependence structure of the running variable within clusters, and (iii) assumptions on the within-cluster covariance structure of the outcome. Within these frameworks, we find that distinct convergence rates and optimal bandwidth choices are qualitatively different from i.i.d. settings. Our results complement the analysis of hansen2019asymptotic by considering nonparametric models. In contrast to their analysis, where the estimators are based on the full sample, the localization of the RD estimator leads to non-standard convergence rates.
Second, we consider estimation of the conditional variance of the local linear RD estimator. We show that the conventional clustered regression residual-based standard error is consistent under the same cluster size conditions that ensure asymptotic normality of the RD estimator. However, in settings with independent data, the regression residual-based standard errors are known to exhibit less finite-sample bias than the so-called nearest-neighbors standard errors. A naive adaptation of the nearest-neighbors standard error for independent data to clustered settings is invalid, and to our knowledge, no valid nearest-neighbors-type standard error for clustered RD designs currently exists. As the second main contribution, we propose a novel clustered nearest-neighbors (CNN) standard error for RD estimators. Our proposed method chooses nearest neighbors taking into account the clustering structure and exploiting independence between clusters. We establish consistency of our proposed CNN standard error under our high-level assumptions on the cluster sizes.\footnote{Although our CNN standard error is developed for RD designs in this paper, the general idea behind it extends naturally to other conditional inference problems under misspecification, such as those studied by abadie2014inference.}
We complement our theoretical analysis with empirical applications that assess the finite-sample performance of the proposed standard error relative to existing alternatives.
Cluster-robust inference is routinely employed in empirical RD designs. Despite this fact, formal results remain limited, even for the standard nonparametric regression using local polynomial estimators. For RD designs with clustering, bartalotti2017regression show asymptotic normality of local polynomial RD estimators under the assumption that all clusters are of the same size and the realizations of the running variable are on the same side of the cutoff within each cluster. Clustered standard errors have been implemented in popular RD packages {\tt rdrobust} and {\tt RDHonest} without much theoretical foundation. Our paper contributes to this literature by providing a unified theory for all common RD variants with arbitrary clustering.
lin2000nonparametric, wang2003marginal, and bhattacharya2005asymptotic study general local polynomial estimators under bounded cluster sizes. In these regimes, as the bandwidth converges to zero, the probability of having more than one unit within the estimation window converges to zero for any given cluster, and in consequence, the clustering does not affect the asymptotic distribution. shimizu2024nonparametric studies nonparametric density and local polynomial estimation under clustered sampling with heterogeneous and potentially unbounded cluster sizes. However, his framework imposes restrictions that rule out several empirically relevant RD settings. In particular, cluster sizes within the estimation window are required to remain uniformly bounded in expectation, the covariates are not allowed to be perfectly dependent within a cluster, and the within-cluster covariance structure of the residuals takes on a very specific form. Furthermore, his proposed variance estimator relies on parametric assumptions. By contrast, our framework accommodates weaker conditions on cluster sizes, allows for more general dependence structures in the running variable within clusters, and imposes milder assumptions on the covariance functions, while still delivering asymptotic normality. Moreover, our proposed nearest-neighbor standard error is fully nonparametric and does not rely on parametric modeling assumptions.
In settings of regular parameters, there is a vast literature about cluster-robust inference. liang1986longitudinal, white2014asymptotic, arellano1987computing provide the foundational large-\(G\) theory for cluster-robust covariances for bounded cluster sizes; see cameron2015practitioner, mackinnon2023cluster for reviews. djogbenou2019asymptotic, hansen2019asymptotic, bugni2025inference, hansen2025jackknife extend the results to unequal, possibly large clusters. abadie2023should discuss clustering adjustments from a design-based perspective. chiang2025genuinely point out that in many empirical applications the size of the largest cluster is not negligible relative to the total sample size, and they provide an alternative bootstrap inference procedure.
This paper is the first to propose and formally study a cluster-robust nearest-neighbors-type standard error. In this regard, we extend the work of abadie2014inference, who showed consistency of the nearest-neighbors standard errors under i.i.d.\ sampling (in a more general class of misspecified models).
Section (ref) introduces the model and gives a preview of our main results. Section (ref) states the high-level assumptions and establishes asymptotic normality of the local linear RD estimator. Section (ref) studies the asymptotic frameworks. Section (ref) introduces our proposed standard error and we show that it is consistent. Section (ref) contains the empirical applications. All proofs are collected in the Appendix.
We consider a sharp RD design, in which a unit receives the treatment if and only if their running variable exceeds a known cutoff value, which we normalize to zero. The observed data is divided into $G$ clusters, and in each cluster, we observe $n_g$ units. Let $X_{gi}$ and $Y_{gi}$ denote the running variable and the outcome variable of observation $i$ in cluster $g$, respectively. The total sample size is given by $n = \sum_{g \in [G]} n_g$, where $[G] = (1,\ldots,G)$. In our asymptotic analysis, we will treat $G$ and $(n_g)_{g \in [G]}$ as deterministic sequences indexed by the sample size $n$.
Observations in different clusters are independent, but they can be dependent within a cluster. The outcome is generated according to the model
where $\mu(x)= \mathbb E[Y_{gi}|X_{gi}=x]$, $\mathbb{E}[\varepsilon_{gi}|\mathcal X_g] = 0$, $\mathcal X_g = (X_{gi})_{i \in I_g}$, and $I_g = (1,\ldots,n_g)$. The error terms $\varepsilon_{gi}$ can be arbitrarily dependent within a cluster. We denote their variance and covariances, conditional on $\mathcal X_g$ by
and we denote the covariance matrix of $\mathcal Y_g = (Y_{gi})_{i \in I_g}$ conditional on $\mathcal X_g$ by $\Sigma_g$.
Under continuity assumptions on the conditional expectation of the potential outcomes, the jump in the conditional expectation $\mu(x)$ at the cutoff identifies the average treatment effect of units at the cutoff hahn2001identification. Our parameter of interest is therefore given by $$ \tau = \mu(0^+) - \mu(0^-), $$ where for a generic function $f$, $f(0^+)$ and $f(0^-)$ denote the right and left limits of the function $f$ at zero.
In practice, it is common to estimate the RD parameter via local linear RD regressions. This estimator is defined as
where $V_i =(T_i, X_i, T_i X_i,1)^\top$, $k_h(v)=k(v/h)/h$ with $k(\cdot)$ a kernel function and $h>0$ a bandwidth, and $e_1 = (1,0,0,0)^\top$ is the first unit vector. The weights $w_{gi}(h)$ depend only on the running variable and the bandwidth; the exact expressions for the weights are given in Appendix (ref). In our setting, the variance of the RD estimator conditional on the running variable $\mathcal X_n = (\mathcal X_g)_{g \in [G]}$ equals
Our first main result establishes the large-sample behavior of the RD estimator under clustered sampling. We derive it under high-level conditions on the weights $w_{gi}(h)$ that translate into restrictions on cluster sizes within the estimation window. Under these general high-level conditions, we show that the local linear RD estimator is asymptotically normal: \[ se(h)^{-1}\big(\widehat{\tau}(h) - \mathbb{E}[\widehat{\tau}(h)\mid \mathcal{X}_n]\big) \xrightarrow{d} \mathcal{N}(0,1). \] The convergence rate of the conditional variance $se^2(h)$ depends on the covariance of outcomes within each cluster, $\Sigma_g$, the joint distribution of the realizations of the running variable within each cluster, $\mathcal X_g$, and the cluster sizes, and it can be slower than the convergence rate of $(nh)^{-1}$ obtained in i.i.d.\ settings. For example, we show in Section (ref) that if the joint distribution of $\mathcal X_g$ admits a bounded density and some mild regularity conditions hold, then \[ se^2(h) = O_P\left(\frac{1+\lambda_n}{nh}\right), \] where $$ \lambda_n = \frac{h}{n}\sum_{g \in [G]}n_g(n_g-1). $$ If, in turn, the distribution of $\mathcal X_g$ is degenerate, meaning that the realizations of the running variable are equal within each cluster, then \[ se^2(h) = O_P\left(\frac{1+\lambda_n/h}{nh}\right). \]
The above results provide bounds on the convergence rate of the conditional variance $se^2(h)$. Whether these rates are binding depends on the exact form of the conditional covariance matrix $\Sigma_g$. In Section (ref), we give an example where these bounds are achieved, but we note that the convergence rate of $se^2(h)$ can be faster. For example, in the case of bounded joint density, if the residuals are uncorrelated within each cluster, then $se^2(h) = O_P((nh)^{-1})$, as in the i.i.d.\ case, even if $\lambda_n$ diverges to infinity.
Reliable and efficient inference requires a variance estimator that is both consistent and well-behaved in finite samples. Natural estimates of $se^2(h)$ are of the form
with $\widehat\sigma_{g,ij}$ being some estimate of $\sigma_{g,ij}$. Setting $\widehat\sigma_{g,ij}$ to the product of residuals from the local linear RD regression associated with observations $i$ and $j$ in cluster $g$ yields a clustered regression residual-based standard error, analogous to the cluster-robust standard errors proposed for the linear regression liang1986longitudinal. We show that this standard error is valid under the same high-level conditions on the weights as those ensuring asymptotic normality.
Even though the regression residual-based standard error is consistent under mild assumptions, it has been long recognized in i.i.d.\ settings that this approach may be overly conservative in finite samples, and the nearest-neighbors standard error has been proposed to alleviate this issue abadie2014inference. At its core, the nearest-neighbors approach estimates the conditional variance of the outcome variable for any given observation using the outcome variability among its nearest neighbors.
A direct way to adapt the nearest-neighbors variance estimation approach to clustered data is to set $\widehat\sigma_{g,ij}$ in equation (ref) to
where $\mathcal N^\text{naive}_{gi}$ is the set of $J$ units that are closest to unit $i$ in cluster $g$ in terms of the running variable.\footnote{A standard error of that type was proposed by calonico2019regression and is implemented in the software package {\tt rdrobust}. We note that this approach also assumes that residuals are uncorrelated across observations on opposite sides of the cutoff.} However, this procedure does not take into account the clustering structure, and there is no reason to expect that such a variance estimator is consistent in general. First, if the choice of neighbors is associated with cluster membership, then an additional bias term can be present. Second, the bias-correction factor is devised for variance estimation, but it turns out it is not suitable for covariance estimation. We give two illustrative examples showing these problems in Section (ref).
To remedy the deficiencies of the naive approach, we propose a different clustered nearest-neighbors standard error where $\widehat\sigma_{g,ij}$ is set to
where $\mathcal N^{\,1}_{gi}$ and $\mathcal N^{\,2}_{gi}$ are carefully chosen sets of neighbors of observation $i$ in cluster $g$. Crucially, we require that the observations in $\cup_{i \in I_g}\mathcal N^{\,1}_{gi}$ and $\cup_{i \in I_g}\mathcal N^{\,2}_{gi}$ belong to two disjoint sets of clusters not including $g$. Under this requirement and additional assumptions we show this standard error is consistent, and one can see in simulations that our procedure has favorable finite-sample properties, exhibiting a smaller bias relative to the regression residual-based approach is settings where the curvature of $\mu$ is substantial.
In this section, we show asymptotic normality of the local linear RD estimator under high-level assumptions on the weights.
The first assumption controls the cluster sizes in terms of the weights assigned to units in each cluster.
Assumption (ref) puts restrictions on the size of the terms $\sum_{i,j \in I_g} |w_{gi}(h) w_{gj}(h)|$. In our asymptotic analysis, it is used to control the contributions of different clusters to the conditional variance $se^2(h)$. This assumption allows the number of non-zero weights within each cluster to grow with the sample size, as long as none of the clusters dominates all others in terms of the sums of the absolute value of cross-products of the weights. This assumption is very general and we provide a range of different asymptotic frameworks in which it is satisfied in Section (ref). However, we note that it is a sufficient condition to obtain asymptotic normality, but it is not necessary. In particular, it does not cover settings where the number of clusters does not increase with the sample size. For example, if there is one cluster of weakly-dependent data (e.g., strongly mixing process, as extensively studied in time-series settings), one might still obtain asymptotic normality, but part (i) of the assumption will generally not hold.
The second assumption controls the moments of the error terms and is standard for results invoking central limit theorems with non-i.i.d.\ data.
Our first main result establishes asymptotic normality of the RD estimator. To control its bias, we introduce the H{\"o}lder-type class of real functions that are potentially discontinuous at zero, are twice differentiable on either side of the threshold, and whose second derivatives are uniformly bounded by some constant $M>0$: $$\mathcal{F}_H(M)=\{f_1(x)\mathbf{1}\{x\geq 0\}- f_0(x)\mathbf{1}\{x< 0\}: \|f_w''\|{\infty} \leq M, w=0,1\}.$$
Part (i) of the theorem shows that $\widehat\tau(h)$, appropriately recentered and rescaled, is asymptotically normally distributed. Part (ii) shows that the conditional bias is bounded by a quantity that is the product of the bound on the second derivative of the conditional expectation function and an expression that depends only on the weights and the running variable. The bias bound is the same as in the i.i.d.\ setting armstrong2020simple, noack2024bias since the local linear RD estimator is a linear operator and its conditional expectation is not affected by the dependence across outcomes.
Given the standard deviation of the local linear RD estimator and a bound on its bias, we will consider the conditional worst-case mean squared error over data-generating process with $\mu \in \mathcal{F}_H(M)$: \[ \overline{MSE}(h) = \bar b(h)^2 + se^2(h). \] We will explicitly derive the limit of $\overline{MSE}(h)$ and characterize the corresponding asymptotically optimal bandwidth under low-level conditions in the next section.
Our high-level Assumption (ref) introduced in Section (ref) accommodates a wide range of empirical clustering structures. To illustrate this flexibility, we formalize four asymptotic frameworks and provide low-level conditions under which the high-level assumption holds. Each of the asymptotic frameworks is calibrated to a class of empirical RD applications. They differ in (i) how quickly the number of near-cutoff units can grow inside a cluster, (ii) the dependence structure of the running variable within clusters, and (iii) assumptions on the conditional covariance matrix of the outcome.
Asymptotic Frameworks I and II are motivated by empirical applications where there are many clusters that can be potentially large, but the number of units from any cluster within the estimation windows is relatively small. They cover two distinct dependence patterns in the running variable. While Asymptotic Framework I assumes that the joint distribution of the running variable within each cluster admits a bounded joint density, Asymptotic Framework II does not put any restrictions on this distribution at the cost of imposing stronger restrictions on the rates of the cluster sizes. Asymptotic Frameworks III and IV relax the restrictions on cluster sizes imposed in Asymptotic Frameworks I and II, respectively, while instead requiring additional assumptions on the covariance structure of the outcome residuals. We discuss representative empirical examples of each of the asymptotic frameworks in Section (ref).
We will first present general assumptions that are maintained in all the four asymptotic frameworks.
Part (i) of Assumption (ref) matches standard conditions for local polynomial estimation with continuously distributed regressors fan1996local. We impose continuity (and continuity at the cutoff) only to obtain simple closed-form variance limits; our high-level results can also accommodate discrete, mixed, or cutoff-discontinuous running variables under suitable modifications to the variance calculations. In the same spirit, we assume the density is continuous at the cutoff to simplify the formulas; formally, an RD analysis can proceed without this requirement, and our main results remain unaffected if the running variable is discontinuous at the threshold.
The kernel and bandwidth requirements in parts (ii) and (iii) are standard in the nonparametric regression literature. We note that the assumption that the bandwidth shrinks to zero is not necessary for our high-level conditions to hold. We could accommodate a fixed bandwidth; we maintain the usual condition that the bandwidth shrinks to zero solely to obtain closed-form expressions for the leading bias and variance terms. In the asymptotic results, we will use the following kernel constants. For $j \in \mathbf{N}$, let $\bar\mu_j = \int_0^1 k(v)v^j dv$. Further, define $\bar\mu = (\bar{\mu}_2^2 - \bar{\mu}_1 \bar{\mu}_{3})/(\bar{\mu}_2\bar{\mu}_0-\bar{\mu}_1^2)$ and $\bar\kappa = \int_0^1 \bar k(v)^2 dv$, where $\bar k(v) = k(v)(\bar\mu_2 - \bar\mu_1 v ) / (\bar\mu_2\bar\mu_0 - \bar\mu_1^2)$.
Part (i) is slightly weaker than Assumption (ref). Part (ii) is imposed to rule out degenerate conditional dependence structures in the outcome variable, such as situations where two residuals within the same cluster are perfectly negatively correlated. The variance of the RD estimator could collapse to zero in such cases, precluding asymptotic normality.
The first two asymptotic frameworks apply to settings where the number of units from any cluster within the estimation window is relatively small.
The number of units within the estimation window from any given cluster depends on the total number of units in the cluster as well as the joint distribution of the running variable within the cluster. In Asymptotic Framework I, we assume that the joint distribution of the realizations of the running variable admits a bounded joint density, while Asymptotic Framework II leaves the joint distribution unrestricted.
Assumption (ref)(i) excludes cases of perfect within-cluster correlation in the running variable where all units share the same realization of the running variable; such extreme dependence is allowed under Assumption (ref) below at the cost of more restrictive rate conditions. Assumption (ref) is similar to Assumption 2 of hansen2019asymptotic, who study the convergence of the average of clustered observations and other regular, full-sample estimators. Since we consider estimation using the data close to the cutoff, our conditions are formulated in terms of the “local sample size” $nh$ and the “local cluster sizes” $n_gh$, rather than the full sample size $n$ and the cluster sizes $n_g$ used by hansen2019asymptotic. Another conceptual difference is that part (iii) includes an additional \(\log G\) term. This extra factor is due to the randomness of cluster sizes within the estimation window in our framework, whereas hansen2019asymptotic consider the setting of regular parameters and use all observations within each cluster. We note that this asymptotic framework shares some similarities with the framework of shimizu2024nonparametric specialized to one continuous covariate. He also assumes that the joint distribution of the realizations of the running variable admits a joint density (albeit only for subsets of four observations) and $\lambda_n=O(1)$. However, he imposes a restrictive assumption that $\max_{g \in [G]} n_g h = O(1)$, whereas our framework allows this quantity to diverge.
Assumption (ref) imposes no restrictions on the joint density of the running variable, but the rate conditions on the cluster sizes are more stringent than in Assumption (ref). In particular, it allows for the realizations of the running variable to be equal within each cluster. It turns out that this extreme case determines the restrictions on the cluster sizes that are necessary to verify Assumption (ref).
To illustrate the restrictions imposed in Assumptions (ref) and (ref), we consider two examples that differ in the degree of allowed heterogeneity in cluster sizes. We revisit these examples in the next subsection to facilitate comparisons between asymptotic frameworks.
We note that in Example 1, the expected number of units within the estimation window remains bounded in any given cluster. The next example shows that both Asymptotic Frameworks I and II allow the maximal expected number of observations within the estimation window to diverge for some clusters.
The following proposition presents our key theoretical results for Asymptotic Frameworks I and II.
Proposition (ref) verifies our high-level condition on the weights, shows that the conditional variance of the RD estimator is of order $(nh)^{-1}$, and it characterizes the leading term of the conditional worst-case bias $\bar b(h)$. It follows that worst-case asymptotic mean squared error is of order $h^4 + (nh)^{-1}$, which is minimized for bandwidths of order $n^{-1/5}$. With such a bandwidth, the estimator converges at the rate $n^{-2/5}$, given that our assumptions hold for this bandwidth choice. We note that the order of the variance is the same as in the i.i.d.\ case, but its exact form is in general affected by clustering; we derive its limit under additional assumptions in Section (ref).
While Asymptotic Frameworks I and II cover many relevant clustering patterns, the conditions on the growth rates of the cluster sizes might be restrictive in some applications. In this section, we show these rate conditions can be significantly relaxed under direct assumptions on the standard error. We then argue in Section (ref) that these assumptions are reasonable in many settings.
Asymptotic Framework III considers cases where the joint distribution of the realizations of the running variable is continuous, as in Asymptotic Framework I, but it relaxes the rate conditions imposed on the cluster sizes. This asymptotic framework is motivated by settings where the clusters have a large number of units even in a shrinking neighborhood of the cutoff.
Part (i) of Assumption (ref) coincides with part (i) of Assumption (ref). However, the rate conditions on cluster sizes in part (ii) are considerably weaker than those imposed in parts (ii) and (iii) of Assumption (ref). In particular, $\lambda_n$ is allowed to diverge to infinity within this framework.
Asymptotic Framework IV applies to settings where each cluster may contain many units whose realizations of the running variable might be highly correlated. As a result, even within a shrinking neighborhood of the cutoff, some clusters can contribute a large number of units concentrated in a narrow region of the support of the running variable. Such settings occur naturally, for example, if the outcomes are measured at the individual level, whereas the running variable is assigned at the cluster level.
The rate conditions on cluster sizes are significantly weaker than those in Assumption (ref) as we combine these conditions with an additional assumption on the standard error. We note that the rate conditions in Assumption (ref) are stronger than those in Assumption (ref); this is needed because we do not impose any assumptions on the joint distribution of the running variables. We next illustrate the implications of these rate conditions in our examples studied above.
The following proposition presents our key theoretical results for Asymptotic Frameworks III and IV.
Proposition (ref) verifies our high-level Assumption (ref) and characterizes the leading term of the conditional worst-case bias $\bar b(h)$. We note that the rate of the optimal bandwidth and the resulting convergence rate of the estimator are different than in the i.i.d.\ case or the case of small effective cluster sizes discussed in the previous section. Specifically, in Framework III, we have that \[ \widehat \tau(h) - \tau = O_P\left(h^2 + \frac{1}{\sqrt{nh}} + \sqrt{ \frac{\sum_{g \in [G]}n_g^2}{n^2} } \right). \] In contrast to Asymptotic Framework I, the third term in the remainder on the right-hands side can dominate the other terms. Since the bandwidth $h$ appears only in the first two terms, the convergence rate is optimized if $h \sim n^{-1/5}$, in which case we obtain $\widehat \tau(h) - \tau = O_P(n^{-2/5} + \big(\sum_{g \in [G]}n_g^2 / n^2 \big)^{1/2} )$. We note that if $n^{-2/5} = o\big( \big( \sum_{g \in [G]}n_g^2 / n^2 \big)^{1/2} \big)$ and $h\sim n^{-1/5}$, then the variance dominates the bias under the optimal bandwidth choice, such that \[ \frac{\widehat\tau(h) - \tau}{se(h)} \xrightarrow{d} \mathcal N(0,1), \] assuming that Assumption (ref) holds for this bandwidth choice.
In Framework IV, in turn, we have that \[ \widehat \tau(h) - \tau = O_P\left(h^2 + \frac{1}{\sqrt{nh}} \left(1 + \sqrt{ \frac{\sum_{g \in [G]}n_g^2}{n} } \right) \right). \] The fastest possible convergence rate is achieved for $ h \sim \big(\big(1 + \sum_{g \in [G]}n_g^2/n\big)/n \big)^{1/5}$, yielding $\widehat \tau(h) - \tau = O_P\big( n^{-2/5} + \big(\sum_{g \in [G]}n_g^2/n^2\big)^{2/5} \big) $, given that this bandwidth choice satisfies Assumption (ref).
To illustrate the behavior of the local linear RD estimator with clustered data, we will further consider a simplified setup, where exact characterization of the limit of the variance is possible. The following assumption formalizes the setup.
Assumption (ref) imposes a common covariance between any two units and across clusters, and it imposes continuity of the conditional variance and covariance functions.\footnote{This type of assumption can be rationalized by models studied in the functional data literature zhang2007statistical. shimizu2024nonparametric also imposes such an assumption.}
We first study the standard error under Asymptotic Frameworks I and III.
Lemma (ref) imposes only the weak rate conditions on the cluster sizes of Asymptotic Framework III, which, in particular, cover Asymptotic Framework I. First, the lemma provides an upper bound on the convergence rate of the conditional variance $se^2(h)$, and then it derives its exact limit in a special case. Part (ii) demonstrates that the conditional variance is asymptotically equal to the sum of two terms. The first one is driven by the variances of individual units and is the same as in the i.i.d.\ case. The second component is due to the covariance between outcomes within a cluster; we note this part depends neither on the bandwidth nor the choice of the kernel function. If $\lambda_n$ converges to a positive constant, then the two terms are of the same order. If $\lambda_n$ converges to zero, then the effect of clustering on variance becomes asymptotically negligible; a similar result was obtained by shimizu2024nonparametric under stronger rate restrictions on the cluster sizes. If $\lambda_n$ diverges to infinity and $\sum_{\star, \diamond \in \{+,-\}} \mathbf{1}^\pm_{\{\star = \diamond\}} \sigma(0^\star, 0^\diamond) \neq 0$, the covariance part dominates.
Under the assumptions of part (ii) of Lemma (ref), the worst-case mean squared error of $\widehat\tau(h)$ satisfies \[ \overline{MSE}(h) = \left(M^2\bar\mu^2h^4 + \frac{V_1}{nh} + \frac{\sum_{g \in [G]} n_g^2}{n^2} \sum_{\star, \diamond \in \{+, -\}} \mathbf{1}^\pm_{\{\star=\diamond\}} \sigma(0^\star,0^\diamond) \frac{ f(0, 0)}{f_X(0)^2} \right)(1+o_P(1)) , \] where $V_1 = \frac{\bar \kappa}{f_X(0)} \sum_{\star \in \{+,-\}} \sigma^2(0^\star)$. The bandwidth minimizing the leading term is given by $$ h^*_1 = \left(\frac{V_1}{4M^2\bar\mu^2}\right)^{1/5} n^{-1/5}, $$ given that this bandwidth choice satisfied the assumptions of part (ii) of Lemma (ref). The optimal bandwidth $h^*_1$ is the same as the AMSE-optimal bandwidth in the i.i.d.\ case. We emphasize that this result relies on Assumption (ref); under more general covariance structures, it need not hold.
We now study the standard error under Asymptotic Frameworks II and IV. To derive a closed-form expression for the limit of the conditional variance in this setting, we need to impose an additional assumption on the joint distribution of the running variable. For illustration, we consider the setting where all the realizations of the running variable are the same within cluster.
Lemma (ref) imposes only the weak rate conditions on the cluster sizes of Asymptotic Framework IV, which, in particular, cover Asymptotic Framework II. First, the lemma provides an upper bound on the convergence rate of the conditional variance $se^2(h)$, and then it derives its exact limit in a special case.
Under the assumptions of Lemma (ref), the worst-case mean squared error of $\widehat\tau(h)$ satisfies \[ \overline{MSE}(h) = \left(M^2\bar\mu^2h^4 + \frac{V_{2}}{nh}\right)(1+o_P(1)), \] where $ V_2 = \frac{\bar\kappa}{f_X(0)} \sum_{\star \in \{+,-\}} \big( \sigma^2(0^\star) + \sigma(0^\star,0^\star) \sum_{g \in [G]}n_g(n_g-1)/n \big)$. In this setting, the bandwidth minimizing the leading term of the AMSE-optimal bandwidth is given by $$ h^*_2 = \left(\frac{V_2}{4M^2\bar\mu^2}\right)^{1/5} n^{-1/5}, $$ assuming the assumptions of Lemma (ref) hold for this bandwidth choice.
In this section, we study variance estimation based on the nearest-neighbors and regression residual-based approaches.
We first show why the naive adaptation of the nearest-neighbors approach devised for i.i.d.\ settings is in general not valid with clustered data. We then introduce our proposed clustered nearest-neighbors (CNN) standard error and show its consistency.
To highlight the main problems with the naive clustered nearest-neighbors standard error described in Section (ref), we consider a simple setup with $\mu(X) = 0$, $\operatorname{Var}(Y_{gi}|\mathcal X_g) = \sigma^2$, and $J=1$ nearest neighbor. We discuss two examples that differ in the assumptions on the joint distribution of the realizations of the running variable.
First, suppose that each cluster consists of only two observations and $X_{g1}=X_{g2}$ for all $g \in [G]$, such that $\mathcal N^{\text{naive}}_{g1} = \{(g,2)\}$ and $\mathcal N^{\text{naive}}_{g2} = \{(g,1)\}$. Then
Clearly, this standard error cannot be consistent. While this is a very specific example, the same type of problem arises more generally whenever the realizations of the running variable are highly concentrated within clusters.
Second, suppose that the joint distribution of the running variable within each cluster admits a bounded joint density. Consider two observations $i,j \in I_g$, $i \neq j$, and let $(g_1, i')$ and $(g_2, j')$ be their respective nearest neighbors. Assume further that the clusters $g$, $g_1$, and $g_2$ are pairwise different, which occurs with high probability in this setup if all clusters are relatively small. Then, while it is easy to see that the conditional variance estimate is correctly centered, $\mathbb E\big[ \big(Y^{\Delta, \text{naive}}_{gi}\big)^2 |\mathcal X_n\big] = \sigma^2$, the conditional covariance estimates are biased:
It follows that $\widehat{se}_{\text{naive}}^2(h)$ is not correctly centered in general.\footnote{We note that the problem arising in the naive conditional covariance estimation can be easily resolved by leaving out the correction factor when $i \ne j$. However, a complete proof of approximate unbiasedness of the standard error would still require controlling the probability of the respective clusters being distinct across all observations, limiting the applicability of this method to relatively small clusters with a bounded joint density of the realizations of the running variable.}
It is evident from the preceding discussion that selecting neighbors from distinct clusters induces certain independence restrictions that are instrumental for establishing conditional unbiasedness of the standard error. We leverage this insight to construct our proposed standard error, explicitly enforcing the desired independence structure.
For every cluster $g\in [G]$, we define two sets of its “companion clusters”, $\mathcal R_g^1 \subset [G]$ and $\mathcal R_g^2 \subset [G]$. Next, for every $i \in I_g$ and $d \in \{1,2\}$, we define $\mathcal N_{gi}^{d}$ as the set of at least $J$ nearest neighbors of unit $i$ in cluster $g$ in terms of the running variable that are on the same side of the cutoff as $X_{gi}$ and belong to a cluster in the set $\mathcal R_g^d$. Our proposed clustered nearest-neighbors (CNN) standard error is then defined as:
To show consistency of this standard error, we require the sets of nearest neighbors to satisfy two properties. First we require that the neighbors in the two sets $\mathcal N_{gi}^{\,1}$ and $\mathcal N_{gi}^{\,2}$ are independent of the units in cluster $g$ and of each other. This is achieved by selecting companion clusters that do not include cluster $g$, $g \notin \mathcal R^1_g \cup \mathcal R^2_g$, and are disjoint, $\mathcal R^1_g \cap \mathcal R^2_g = \emptyset$. These properties ensure that our standard error is asymptotically conditionally unbiased.
Second, we need to control the dependence between $\sum_{i,j \in I_g} w_{gi}(h) w_{gj}(h) Y^{\Delta_1}_{g i} Y^{\Delta_2}_{g j}$ across different clusters. We achieve that by imposing that each cluster can be a companion cluster for at most $R$ clusters for some fixed number $R$, i.e., for all $g \in [G]$, we require
These general requirements on the choice of companion clusters can be satisfied by many different selection methods. We provide one concrete algorithm in Appendix (ref).
As is standard for nearest-neighbors-type estimators, our method relies on the assumption that the nearest neighbors selected in the construction of $\widehat{se}_{\scriptscriptstyle CNN}^2(h)$ are uniformly close to the respective units.
This assumption is automatically satisfied whenever the nearest neighbors are selected within the estimation window, as is the case for the algorithm given in Appendix (ref), and the bandwidth converges to zero. Since all matched units then lie within distance $h$ of each other, we have $D(h) = O(h)$ by construction. In general, we expect $D(h)$ to converge to zero at a much faster rate. For illustration, consider a setting where the marginal density of $X_{gi}$ is bounded and bounded away from zero in a neighborhood of the cutoff. If the realizations of the running variable are constant within each cluster and $Gh \to \infty$, then we expect that $D(h) = O_P(\log(Gh)/G)$. As a different example, consider a setting where the joint density of $\mathcal X_g$ is continuous for all $g \in [G]$ and $n_{\min}h\equiv\min_{g \in [G]} n_gh \to \infty$, then we expect that $D(h) = O_P(\log(n_{min}h)/n_{min})$.
Our second main result states that $\widehat{se}_{\scriptscriptstyle CNN}^2(h)$ is consistent for $se^2(h)$.
In this subsection, we show consistency of the clustered regression residual-based (CRR) standard error for the local linear RD estimator. For $\star \in \{+,-\}$, define $b_0^\star = \mu(0^\star)$ and $b_1^\star = \mu'(0^\star)$, and let $\hat b_0^\star$ and $\hat b_1^\star$ denote the intercept and slope coefficient on the respective side of the cutoff in the local linear RD regression in equation (ref). The CRR standard error is defined as \[ \widehat{se}_{\scriptscriptstyle CRR}^2 = \sum_{g \in [G]} \sum_{i,j \in I_g} w_{gi}(h)w_{gj}(h) \hat\varepsilon_{gi} \hat\varepsilon_{gj},\quad \hat\varepsilon_{gi} = Y_{gi} - \widehat\mu(X_{gi}), \] where $\widehat\mu(x) =(\hat b_0^- + \hat b_1^- x) \mathbf{1}\{x < 0\} + (\hat b_0^+ + \hat b_1^+ x) \mathbf{1}\{0 \leq x\} $.
Theorem (ref) imposes the same high-level assumption on the weights $w_{gi}(h)$ as we used to show consistency of the CNN standard error. The consistency requirements for $\hat b_0^\star$ and $\hat b_1^\star$ are very mild and they are satisfied in all the considered asymptotic frameworks given our smoothness assumption. The main difference in assumptions relative to the result for the CNN standard error is that Theorem (ref) requires the bandwidth to converge to zero, while the assumptions of Theorem (ref) may hold even for an asymptotically fixed bandwidth.
In this section, we apply our clustered nearest-neighbors (CNN) standard error in four empirical applications, and we compare it to four alternative standard errors. The first two, the classical nearest-neighbors (NN) and Eicker-Huber-White (EHW) standard errors, do not account for clustering. The third approach is the naive clustered nearest-neighbors (Naive CNN) approach described in Section (ref), and the fourth is the clustered regression residual-based (CRR) standard error described in Section (ref).
In order to connect the asymptotic frameworks introduced in Section (ref) to observable features of the data, we first propose a simple diagnostic rule of thumb for assessing whether clusters are sufficiently small and reasonably balanced for Asymptotic Framework I or II to provide accurate approximations to the underlying data-generating process. We then revisit four recent RD applications, each exemplifying one of the asymptotic frameworks.\footnote{For each application, we use the data provided in the respective replication package and consider one of the main RD specifications reported in the paper. For simplicity, we ignored any additional covariates that were included to improve estimation precision. We fix the bandwidth at the value used by the authors of the original study. To simplify the analysis, when the original paper used two different bandwidths on each side of the cutoff, we chose the bigger one.}
In this section, we provide a practical rule of thumb to assess whether our high-level Assumption (ref) on the weights can be plausibly satisfied in a given empirical application. When our conditions hold, the data-generating process may be well approximated by Asymptotic Frameworks I or II under mild assumptions on the covariance of the residuals. If, in turn, these conditions are violated, one may need to justify stronger assumptions on the dependence structure of the residuals to ensure that the standard error converges at a suitable rate, as discussed in our Asymptotic Frameworks III and IV.
To describe our proposed criterion, for $g \in [G]$, define \[ w_{g,\textnormal{ratio}}(h) = \frac{\sum_{i,j \in I_g} \lvert w_{gi}(h) w_{gj}(h) \rvert} {\sum_{g \in [G]} \sum_{i \in I_g} w_{gi}(h)^2}, \] and let \[ w_{\max}(h) \equiv \max_{g \in [G]} w_{g,\textnormal{ratio}}(h), \qquad w_{\textnormal{sum}}(h) \equiv \sum_{g \in [G]} w_{g,\textnormal{ratio}}(h). \] If the conditional covariance matrix of the residuals has eigenvalues bounded away from zero, the conditions $w_{\max}(h) = o_{P}(1)$ and $w_{\textnormal{sum}}(h) = O_{P}(1)$ are sufficient to ensure that our high-level Assumption (ref) is satisfied. To operationalize these asymptotic conditions, in finite samples, one needs to choose some threshold values $\eta_{\text{max}}$ and $\eta_{\text{sum}}$, and check whether $w_{\max}(h) \le \eta_{\max}$ and $w_{\textnormal{sum}}(h) \le \eta_{\textnormal{sum}}$. We consider $\eta_{\max} = 0.1$ and $\eta_{\textnormal{sum}} = 10$ to be reasonable benchmark values in practice.\footnote{To give a heuristic interpretation of these threshold values, we note that, asymptotically, $w_{\max}(h) \approx \max_{g \in [G]} n_{g,h}^2 / n_h$ and $w_{\text{sum}}(h) \approx \sum_{g \in [G]} n_{g,h}^2/n_h$, where \( n_{g,h} \) and $n_h$ denote the number of units from cluster \( g \) and the total number of units within the estimation window, respectively. Now suppose that within the estimation window, there are 100 clusters with 10 observations each, a setting in which asymptotic normality can plausibly be a good approximation. Then $w_{\max}(h) \approx 0.1$ and $w_{\text{sum}}(h) \approx 10$.}
In this section, we revisit four empirical applications that motivated the asymptotic frameworks introduced in Section (ref). Figure (ref) illustrates the number of clusters, the distribution of cluster sizes, and the distribution of the running variable within a cluster. Each point represents an individual observation. For discrete outcome variables, values are jittered to improve visual clarity and mitigate overplotting. In selected clusters, all observations are displayed in a common color and use the same marker symbol to emphasize cluster membership. The dashed lines mark the bandwidths used.
delvalle2020rules study the impact of Mexico's indexed disaster fund (Fonden) on post-disaster economic recovery. The outcomes are constructed as changes in log night lights in the year following a disaster, measured using satellite-based night lights at the municipality level. They leverage a fuzzy regression discontinuity design where the eligibility for disaster transfers depended on whether the realized rainfall exceeded a pre-specified cutoff. The data contains information on municipal requests for Fonden funding in the period between 2004 and 2012. The authors cluster the standard errors at the municipality level. The estimation window contains around 1000 municipalities with an average of 1.5 requests per municipality.
granzier2023coordination study French two-round elections to evaluate how candidates' first-round ranking affects their second-round results. Specifically, they measure the impact of barely achieving a higher rank on remaining in the race or winning, and the running variable is the first-round vote margin between adjacent candidates. For concreteness, we focus on the effect of ranking 1st vs 2nd in the first round on the probability of running in the second round. The data contain information on electoral races from several decades of local and parliamentary elections. The standard errors are clustered at the district level. There are about 2300 clusters in the estimation window, with an average of 3 observations per cluster. By construction, the realizations of the running variable are symmetric around the cutoff: for every observation with running variable value $X_{gi}$, there is a corresponding observation with running variable value $-X_{gi}$. Such degenerate distributions of the running variable are allowed under our Asymptotic Framework II.
wasserman2021up studies the causal effect of an electoral defeat on subsequent political participation using a close-election RD design. The analysis focuses on first-time candidates for US state legislative offices. The running variable is the candidate's margin of victory, and the primary outcome of interest is whether the candidate runs again for any state legislative office within four years of the initial run. The data covers state legislative elections over several decades across the United States, and the standard errors are clustered at the state level. This yields clusters with a large number of observations spread across the support of the running variable, with approximately 250 observations per state on average within a local neighborhood of the cutoff.
johnson2020regulation studies the deterrence effects of a policy under which the Occupational Safety and Health Administration (OSHA) issues press releases about violations of workplace safety and health regulations that exceed a penalty threshold. In this RD design, the running variable is the penalty amount assigned at inspection, with a discontinuity at the press-release threshold. The primary outcome is the count of violations recorded in subsequent inspections. In the dataset, the facilities are organized into “peer groups”--groups of facilities in the same sector located within a 5 km radius of a facility where a penalty was levied---such that all facilities within a peer group share the same penalty value, while the outcome is measured at the facility level. The standard errors are clustered at the peer group level. There are 707 peer groups of facilities, each containing roughly 16 facilities on average in a local neighborhood around the cutoff.
The results for all four empirical applications are reported in Table (ref). Before discussing the values of the different standard errors, we first note that, according to the rule-of-thumb diagnostics reported in the last two columns of the table, the studies of delvalle2020rules and granzier2023coordination represent settings with relatively small and fairly balanced cluster sizes. These settings appear to fit into our Asymptotic Frameworks I and II well, suggesting that asymptotic normality is likely a reasonable approximation to the finite-sample distribution of the RD estimator. By contrast, in the applications studied by wasserman2021up and johnson2020regulation, cluster sizes are either large or markedly unbalanced. In such cases, justifying asymptotic normality requires an additional assumption on the convergence rate of the conditional variance $se^2(h)$; Assumption (ref) provides one illustrative sufficient condition for it to hold.
The standard errors computed under the assumption of i.i.d.\ sampling are substantially smaller than their clustered counterparts in applications with large clusters, which suggest that they fail to account for relevant within-cluster dependence. Our proposed CNN standard error is close in magnitude to the conventional CRR approach, while the Naive CNN standard error is markedly smaller in the fourth application, consistent with the discussion in Section (ref).
This paper proposes a general framework for sharp RD designs with clustered data. Under general high-level conditions, we establish the asymptotic normality of the local linear RD estimator, and we illustrate the high-level conditions in empirically motivated asymptotic frameworks. Furthermore, we develop a novel nearest-neighbors-type standard error tailored to clustered samples. Our approach is easily extendable: with minor modifications, it readily accommodates fuzzy RD and kink designs as well as settings with covariate adjustments.