A unified test for regression discontinuity designs

A unified test for regression discontinuity designs

\if0 {2pt} {2pt} \fi

Introduction

Diagnostic tests are vital to regression discontinuity (RD) design. The typical procedure evaluates the testable restrictions which are expected to hold under a non-testable identifying restriction hahnIdentificationEstimationTreatment2001. All such testable restrictions are expected to hold whenever identification holds. Hence, one must test the null hypothesis of all restrictions hold against the alternative hypothesis of at least one restriction fails.

However, these restrictions are often tested separately without appropriate size control. Throughout this study, we consider the continuity-based framework. \footnote{In Cattaneo_Idrobo_Titiunik_2019(Cattaneo_Idrobo_Titiunik_2019,Cattaneo_Idrobo_Titiunik_2023), the continuity-based framework is a contrasting concept against the alternative local randomization framework for a finite sample analysis in which explicit randomization is considered within a small range as studied in Cattaneo_Frandsen_Titiunik_2015, Cattaneo_Titiunik_Vazquez-Bare_2016, cattaneoComparingInferenceApproaches2017.} In this framework, most empirical studies test these restrictions as null hypotheses for the density test mccraryManipulationRunningVariable2008 and the balance or placebo test leeRandomizedExperimentsNonrandom2008. These tests are well-established separately but are rarely evaluated jointly. Hence, the multiple-testing problem plagues size control. Consequently, empirical studies may over-reject the underlying null hypothesis that all restrictions hold under identification.

In this study, we document the severe over-rejection of the diagnostic restriction in a meta-analysis of $59$ empirical RD studies published recently in the top five economics journals. \footnote{See Online Appendix B for the search criterion.} We found that $19$ out of $59$ studies rejected at least one testable restriction, leading to the underlying restriction being questioned in $32\%$ of the studies.

Nevertheless, each test appeared valid. Among $787$ tests run separately, only $4.3\%$ of the restrictions were rejected, and their p-values uniformly distributed over $[0,1]$, suggesting that the test statistics are drawn from the null hypothesis that all testable restrictions hold. Hence, the underlying identification was plausible; nonetheless, we ended up with an incorrect conclusion for $19$ studies rather than three.

We conclude that this over-rejection was caused by multiple testing problems. In $59$ studies, the median number of tests was 12 per study. Furthermore, we found only five studies out of $59$ reporting some type of joint testing. No study has reported a joint test that combines the density and balance or placebo tests because no test has been proposed to handle them jointly.

In this study, we propose a unified diagnostic test to resolve this over-rejection problem. This involves two major challenges. First, no tests consider the joint nature of the density and balance tests. Second, the test statistics are nonparametric estimates. Only one of the five studies considered the nonparametric nature of the test statistics. To address these two challenges, we derived the joint asymptotic normality of the nonparametric density estimation and nonparametric conditional mean estimation. We find that nonparametric density estimates and conditional mean estimates are asymptotically orthogonal. Based on this finding, we provide the first unified test with density and balance test statistics based on our new theoretical results from their nonparametric bias-corrected estimates.

Our unified test aggregates the vector of the density and balance test statistics; however, possible aggregation is not unique. An appropriate aggregation depends on the underlying alternative hypothesis. We considered two different aggregations: (1) a Wald (chi-squared) statistic of the sum of the squared statistics of the vector and (2) a Max test statistic of the maximum squared statistics of the vector. The two test statistics complement each other in detecting different types of alternative hypotheses.

The Wald statistic detects the alternative hypothesis that a large fraction of the original vector deviates from the null hypothesis. The Max statistic detects the alternative hypothesis that only a small fraction of the vector deviates from the null hypothesis. Nevertheless, the conventional Wald statistic did not attain the nominal size, possibly because of poor inverse estimation of the variance matrix. We bypass this problem by standardizing only the diagonal elements instead of multiplying the inverse of the variance matrix and propose it as the standardized Wald (hereafter, sWald) statistic.

In numerical simulations, the sWald and Max statistics exhibit appropriate size control properties for a moderate number of covariates when the sample size is as large as $1000$. In particular, the sWald test is superior to other methods in size control performance for any number of covariates up to $25$ in cases where even the Bonferroni correction may fail to achieve size control. \footnote{We conjecture that the reason for this failure is that Bonferroni correction with $25$ covariates evaluates the thin tail of the nonparametric test statistic distribution, which could be imprecisely estimated.}

We verified the superior power properties of the proposed methods. The Max test is superior to the Bonferroni correction in its power against the alternative hypothesis that only one covariate sees a jump in the conditional mean. The sWald test is superior to the Bonferroni correction in its power against the alternative hypothesis that all covariates see jumps in their conditional means. In conclusion, we recommend the Max test if a user has a few particularly concerning covariates. Otherwise, we recommend the sWald test for better size control and power in other scenarios. \footnote{We also compare our methods with the resampling based method of romanoExactApproximateStepdown2005, romanoEfficientComputationAdjusted2016 using a stata package rwolf2 clarkeRwolf2ImplementationFlexible2021 in Online Appendix D. We note that the applicability of romanoExactApproximateStepdown2005 is not necessarily guaranteed for the RD context, and rwolf2 can accept the balance tests only. Hence, our simulation comparisons with rwolf2 are limited to the balance tests only. In the simulation, the Romano and Wolf correction exhibits a bimodal performance: it achieves a good power for the large effect sizes with high correlations, but the power is much lower than our methods when the effect sizes are small. Hence, we continue to support our methods for the purpose of detecting possibly small violations in diagnostic testings.}

This study contributes to the literature on diagnostic procedures for RD design. \footnote{See the latest textbooks Cattaneo_Idrobo_Titiunik_2019(Cattaneo_Idrobo_Titiunik_2019,Cattaneo_Idrobo_Titiunik_2023) for an extensive survey of RD literature.} mccraryManipulationRunningVariable2008 has developed the density test. Several theoretical studies have improved these procedures. \footnote{Several theoretical studies also provide other forms of testing procedures in regression discontinuity designs. For example, hsuTestingTreatmentEffect2019 and hsuTestingMonotonicityConditional2021 provide tests for the presence of treatment effect heterogeneity across covariate values and for monotonicity in the conditional treatment effects; Bertanha_Imbens_2020 provide tests for exogeneity and external validity in fuzzy designs.} otsuEstimationInferenceDiscontinuity2013 have proposed an empirical likelihood-based density test, cattaneoSimpleLocalPolynomial2019 have considered a test based on a local polynomial estimation of the density function, and bugniTestingContinuityDensity2020a have proposed an approximate sign test for the density test. Nevertheless, none of these studies have considered the joint procedure for these tests. Diagnostic test practices follow a similar pattern: separate multiple testing of these diagnostic tests or some joint procedure without considering the nonparametric nature of the test statistics. In this study, we provide unified testing for the null hypothesis that all these diagnostic restrictions hold, suggesting that identification is plausible. We also provide a statistical software rdtest, which is an explicit wrapper of two standard packages, rdrobust (Calonico_Cattaneo_Farrell_Titiunik_2017) and rddensity Cattaneo_Jansson_Ma_2018. Therefore, we provide an easy-to-implement unified diagnostic procedure for the widely used continuity-based framework.

Moreover, we contribute to the asymptotic theory of nonparametric estimates. We show the joint asymptotic normality of multiple local polynomial estimates for the conditional mean and density functions. calonicoRobustNonparametricConfidence2014a show the asymptotic normality of a single local polynomial estimate for a conditional mean function at a boundary point after bias correction. cattaneoSimpleLocalPolynomial2019 also show the asymptotic normality of a local polynomial estimate of a density function at a boundary point. We developed the joint normality results based on these results. The marginal bias and variance expressions are identical to those in calonicoRobustNonparametricConfidence2014a and cattaneoSimpleLocalPolynomial2019, and we impose no restrictions on the choice of bandwidth apart from the original rate conditions. Consequently, we succeeded in accommodating the same mean squared error (MSE) optimal estimates for the balance and placebo tests calonicoRobustNonparametricConfidence2014a \footnote{Recently, another bandwidth selector of coverage error optimal bandwidths has been proposed by Calonico_Cattaneo_Farrell_2020 along with theoretical background calonicoCoverageErrorOptimal2022. In our numerical analyses and implementation, we used the default option for MSE optimal bandwidths.} and density tests cattaneoSimpleLocalPolynomial2019.

The remainder of the paper proceeds as follows. Section (ref) demonstrates the multiple testing problem in the meta-analysis, and Section (ref) explains the joint asymptotic normality of the nonparametric estimates. Subsequently, Section (ref) presents our joint tests, Section (ref) demonstrates the performance of our tests in the simulation, and Section (ref) considers the two extensions for our procedure: equivalence testing which flips the null and alternative hypotheses, and pretesting analysis in view of roth2022pretest. Finally, Section (ref) concludes the article by recommending practices and proposing future research questions. In Online Appendix A , we provide all the formal results and the proofs.

Current practices of diagnostic tests

Most empirical studies diagnose their designs using two continuity restrictions: the density function of the running variable and the conditional mean functions of the pretreatment covariates. In this study, we focus on the two major practices that are routinely evaluated with well-established packages. There are other popular diagnoses such as the placebo cutoff analysis. See cattaneoRegressionDiscontinuityDesigns2022 and Cattaneo_Idrobo_Titiunik_2019 for further diagnostic concepts.

The logic behind the diagnoses is that “the RD design offers an array of empirical methods that, under reasonable assumptions, can provide useful evidence about the plausibility of its assumptions.” Cattaneo_Idrobo_Titiunik_2019 In the continuity-based framework, identification holds when the conditional mean function of the potential outcome is continuous in the conditioning running variable at the cutoff. In other words, the units near the cutoff must be similar in the unobserved characteristics of the potential outcomes.

However, units often know the cutoff value and hence they may manipulate the running variable to be eligible for treatment. A typical idea to assess such manipulation is that “if units lack the ability to precisely manipulate the score value they receive, there should be no systematic differences between units with similar values of the score” Cattaneo_Idrobo_Titiunik_2019. This idea conveys a testable restriction that the observed characteristics must be balanced near the cutoff.

Similar idea is that “if units do not have the ability to precisely manipulate the value of the score that they receive, the number of treated observations just above the cutoff should be approximately similar to the number of control observations below it” Cattaneo_Idrobo_Titiunik_2019. This latter idea conveys another testable restriction that the observed density function of the running variable must be continuous at the cutoff. The former restriction is often conducted as the balance or placebo test, and the latter restriction is evaluated as the density test.

All those testable restrictions stem from the ideas that support identification. Hence, if we claim the plausibility of identification, all those diagnostic restrictions should hold, and none of them must be violated. \footnote{Without additional restrictions, these restrictions have no direct connection with identification. For example, the density test is neither necessary nor sufficient for identification mccraryManipulationRunningVariable2008. See Ishihara_Sawada_2022 for conditions under which the null hypothesis of the density test implies identification.}

In other words, we must test the null hypothesis that all testable restrictions hold against the alternative hypothesis that at least one restriction fails to diagnose the underlying designs. However, most studies have tested these restrictions separately as density tests mccraryManipulationRunningVariable2008, cattaneoSimpleLocalPolynomial2019 and balance or placebo tests leeRandomizedExperimentsNonrandom2008, leeRegressionDiscontinuityDesigns2010. Running these tests separately would over-reject the underlying null hypothesis that all testable restrictions hold because of the multiple testing problem.

We quantified the size control problem by conducting a meta-analysis of RD studies published in the top five economics journals. Among $59$ papers that satisfied our search criteria up to November 2021, we collected $787$ reported test statistics. \footnote{Many studies aim to demonstrate their robustness with different specifications of the same tests with different bandwidths, order of polynomials, or kernels. We screen the main specification for each test, and our reported results are the lower bounds of their over-rejections. See Online Appendix B for details.} Our meta-analysis found two facts: (1) the reported test statistics appear to conform to the null hypothesis that all restrictions hold; (2) nonetheless, $32\%$ ($19$) studies had at least one rejected test.

Fact 1: Test statistics are likely drawn from the null distributions

First, we document the reported test statistics likely to be drawn from null distributions. Table (ref) presents the percentage of tests rejected from the reported ones. Among $787$ tests reported in $59$ studies, we found that less than $5\%$ rejected null hypotheses. Figure (ref) shows the quantile plot of the p-values. The quantile plot aligns with the $45$-degree line, and the distribution is likely to be uniform over $[0,1]$. \footnote{This is also true among the rejected studies. See Appendix Figure 7 for the histogram of p-values among the rejected studies.} Both lines of evidence support the null hypothesis that all restrictions hold and hence identification is plausible.

Fact 2: The underlying null hypothesis is over-rejected because of multiple testing

From Fact 1, approximately three studies out of $59$ should have shown any test rejection because the null hypotheses appear valid. However, Table (ref) suggests that $32\%$ of the studies ($19$ of $59$) rejected at least one null hypothesis. Among $19$ studies, $8$ studies rejected more than once. As shown in Table (ref) and Figure (ref), many studies ran more than $10$ diagnostic tests without appropriate size control. \footnote{Also, studies with many rejections tend to run excessively many tests. See Appendix Figure 8 for the scatter plot of the number of rejected tests and the number of tests.} Hence, we blame the multiple testing problem for over-rejecting the null hypothesis of all restrictions hold.

How does empirical studies respond to their rejections

Upon rejection of their diagnostic tests, most studies defended their statistical significance in diagnostic tests by making some supporting arguments such as stating economical insignificance of their rejected test statistics. Those supporting arguments are indeed desirable. However, the majority of the studies did not address the root cause of their rejections, the multiple testing problem.

Notably, only five studies reported joint testing for balance tests. All but one study did not consider the nonparametric nature of the test statistics. \footnote{We confirm some sort of joint testing in the following five studies: meyerssonIslamicRuleEmpowerment2014, abdulkadirogluEliteIllusionAchievement2014, johnsonRegulationShamingDeterrence2020, Fort_Ichino_Zanella_2020, jagerLaborBoardroom2021. Fort_Ichino_Zanella_2020 use a nonparametric procedure by canayApproximatePermutationTests2018.} Furthermore, no studies have reported joint testing of the balance and density tests because no method is available for the joint test.

We emphasize that only a few studies made explicit argument that testing many hypotheses can result in their rejections of the diagnostic tests. Furthermore, their arguments were all informal without any multiple testing corrections. For example, tyreforshinnerichDemocracyRedistributionPolitical2014 stated that “Only one of the 44 estimates is significant at the 5 percent level. However, that is to be expected since, if 100 specifications are tested, it is likely that five will be statistically significant by chance” (page 981). With a unified diagnostic test, these studies could have formally concluded that their rejections were by chance in their data.

Unified test based on joint asymptotic normality

As shown in Section (ref), the multiple-testing problem causes the over-rejection of the underlying restriction which suggests plausible identification. We introduce a notation to formally state the problem. Let $(Z_{1,i},\ldots,Z_{d,i},X_i)'$, $i=1,\ldots,n$ denote the observed random sample, where $X_i$ is the running variable with density $f(x)$ with respect to the Lebesgue measure and $Z_i=(Z_{1,i},\ldots,Z_{d,i})'$ denote the pre-treatment covariates. Given a known threshold $\Bar{x}$, which we set to $\Bar{x}=0$ without loss of generality, let the running variable $X_i$ determine whether unit $i$ is assigned a treatment ($X_i\geq0$) against ($X_i<0$). In this notation, empirical applications frequently consider two different types of restrictions: the density function $f(x)$ is continuous at $x = 0$, and the conditional mean functions of $Z_k$, $\mu_{Z_k}(x)=E[Z_{k,i}|X_i=x]$ for $k=1,\ldots,d$, are continuous at $x = 0$. Specifically, the following $d+1$ restrictions must hold simultaneously \[ f_+ = f_-, \mu_{Z_1+} = \mu_{Z_1-}, \ldots, \mu_{Z_d+} = \mu_{Z_d-} \] where $f_+ = \lim_{x\to0^+} f(x)$, $f_- = \lim_{x\to0^-} f(x)$ and $\mu_{Z_k+}=\lim_{x\to0^+}\mu_{Z_k}(x)$, $\mu_{Z_k-}=\lim_{x\to0^-}\mu_{Z_k}(x)$. Hence, one must conduct the following test:

to determine whether identification is plausible where $\tau_f = f_+ - f_-$ and $\tau_{Z_k} = \mu_{Z_k+} - \mu_{Z_k-}$. However, if one runs $d+1$ tests separately:

the null hypothesis $H_0$ is over-rejected as documented in Section (ref).

We resolved the size control problem by presenting a unified test in two steps. First, we present the theoretical basis for our unified test. Second, we propose two test statistics for the unified test against multiple testing.

New joint normality result

First, we derive the joint asymptotic normality of the nonparametric local polynomial estimators for $f$ and $(\mu_{Z_1}, \ldots, \mu_{Z_d})$ evaluated at the boundary of the support.

Following calonicoRobustNonparametricConfidence2014a and cattaneoSimpleLocalPolynomial2019, we use local polynomial regressions to approximate the unknown functions $f(x)$ and $\mu_{Z_k}(x)$ for $k=1,\ldots,d$ flexibly near the cutoff $\Bar{x}=0$; see fan1996local for a review of local polynomial regressions. Based on the testing problem, we must obtain two distinct estimands for each function using the two subsamples $\{X_i: X_i\geq 0\}$ (approximation from the right) and $\{X_i: X_i< 0\}$ (approximation from the left) at a boundary point $x = 0$. These estimators are recommended owing to their excellent boundary properties.

Using a polynomial expansion, we construct smooth local one-sided approximations of the sample average of $Z_{k,i}$ and obtain the estimators $\hat \mu_{Z_k+,l}(h_{k,n})$ and $\hat \mu_{Z_k-,l}(h_{k,n})$ as intercepts in the local polynomial regression. Specifically, we estimate $\tau_{Z_k}$ using the following $l$th-order local polynomial estimators: for some $l\geq1$ and for each $k=1,\ldots,d$,

where

$e_0=(1,0,\ldots,0)'$ is the first $(l+1)$-dimensional unit vector, $r_l(x)=(1,x,\ldots,x^l)'$ is a $l$th order polynomial expansion, $K(\cdot)$ is a kernel function, and $h_{k,n}$ is a positive-bandwidth sequence.

Similarly, using a polynomial expansion, we construct smooth local one-sided approximations of the empirical distribution of $X_i$ and obtain the estimators $\hat f_{+,p}(h_{f,n})$ and $\hat f_{-,p}(h_{f,n})$ as the slope coefficients in the local polynomial regression. Specifically, we estimate $\tau_{f}$ using the following $p$th-order local polynomial estimator. For some $p\geq2$,

where

$e_1=(0,1,0,\ldots,0)'$ is the second $(p+1)$-dimensional unit vector; $\tilde{F}(x) \equiv n^{-1}\sum_{i=1}^n$\\$ 1\{X_i \leq x\}$ is the empirical distribution function of $X_i$, and $h_{f,n}$ is the positive-bandwidth sequence.

We derive the asymptotic distributions of the local polynomial estimators based on the assumptions in calonicoRobustNonparametricConfidence2014a and cattaneoSimpleLocalPolynomial2019, as follows.

Assumption (ref) is a set of basic regularity conditions for the data generation process. Assumption (ref) (a) ensures that the density of the running variable is well defined and sufficiently smooth and that the observed values are arbitrarily close to the cutoff for large samples. Assumption (ref) (b) ensures that the conditional expectations and variances of the pre-treatment covariates are well defined and have sufficient smoothness with the existence of the fourth moments.

Assumption (ref) is a set of standard conditions for the kernel function and is satisfied for kernels such as triangular and uniform kernels, commonly used in empirical work. Our results can be extended to accommodate kernels with unbounded supports or asymmetric kernels with more complex notations.

We provide a joint asymptotic distribution for the local polynomial estimators. In the following, let $l,p\in\mathbb{N}$; $h_n=(h_{1,n},\ldots,h_{d,n})'$; $h_{max,n}=\max\{h_{1,n},\ldots,h_{d,n}\}$; $h_{min,n}=\min\{h_{1,n},\ldots,h_{d,n}\}$; we assume that all the limits are $n\to\infty$ (and $h_{k,n}, h_{f,n}\to0$), unless stated otherwise. Whenever there was no confusion, we dropped the sample size $n$ sub-index notation. For derivatives $\mu_{Z_k}^{(s)}(x)$ and $f^{(s)}(x)$, let $\mu_{Z_k+}^{(s)}=\lim_{x\to0^+}\mu_{Z_k}^{(s)}(x)$, $\mu_{Z_k-}^{(s)}=\lim_{x\to0^-}\mu_{Z_k}^{(s)}(x)$, $f_+^{(s)} = \lim_{x\to0^+} f^{(s)}(x)$, and $f_-^{(s)} = \lim_{x\to0^-} f^{(s)}(x)$.

Theorem (ref) shows that the joint asymptotic distribution of $\hat\tau_{f,p}(h_f)$ and $\hat\tau_{Z_k,l}(h_k)$ is a multivariate normal distribution under the same set of conditions as in calonicoRobustNonparametricConfidence2014a for $\hat\tau_{Z_k,l}(h_k)$ and cattaneoSimpleLocalPolynomial2019 for $\hat\tau_{f,p}(h_f)$. Hence, the proposed joint distributional approximation is a natural generalization of the distributional approximation of each local polynomial estimator in these studies.

We emphasize that the covariance matrix of $(\hat\tau_{Z_1,l}(h_1),\ldots,\hat\tau_{Z_d,l}(h_d))'$ and the variance of $\hat\tau_{f,p}(h_f)$ are asymptotically orthogonal. Orthogonality simplifies the estimation of the overall asymptotic covariance matrix because the asymptotic covariance of $\hat\tau_{f,p}(h_f)$ and $\hat\tau_{Z_k,l}(h_k)$ is known to be zero. We derive this orthogonality from the well-known fact that $\operatorname{Cov}(Z_{k,i} - \mu_{Z_k}(X_i) ,g(X_i))=0$ for any function $g$.

Orthogonality affects the nature of the multiple testing problem in various ways. When the test statistics are independent, the multiple testing problem is likely worse. Hence, the size control problem in the naive test procedure becomes worse. In turn, a typical multiple test correction such as Bonferroni correction may become less conservative. For our procedure, we have one less parameter to estimate, and the orthogonality may help stabilizing our procedure.

Proposed unified tests

Given the joint asymptotic normality result, we introduce a unified test for the plausibility of identification. We constructed the test statistics based on the bias-corrected local polynomial estimators $\hat\tau_{Z_k,l+1}(h_k)$ and $\hat\tau_{f,p+1}(h_f)$, where the bandwidth was chosen as MSE-optimal for $\hat\tau_{Z_k,l}(h_k)$ and $\hat\tau_{f,p}(h_f)$. Our unified test has two forms: one aggregates the vector of the density and balance test statistics with the $L^2$-norm, and the other aggregates them via the $L^\infty$-norm.

The Wald (chi-squared) test statistic with the $L^2$-norm has the following form.

where

This conventional test statistic suffers from severe size distortions in finite samples. In particular, the inverse of the asymptotic covariance matrix estimator becomes unstable for a sufficiently large dimension $d$.

Instead, we propose a sWald test statistic that modifies the Wald test statistic. Specifically, we standardize only the diagonal elements instead of multiplying the entire matrix by the inverse of the asymptotic covariance matrix. The sWald test statistic has the following form:

where

and $\operatorname{diag}(a_1,\ldots,a_n)$ is the $(n\times n)$ diagonal matrix with diagonal elements $a_1,\ldots,a_n$.

We obtain the asymptotic distribution of $\tilde \chi^2_{l+1,p+1}(h)$ from the joint asymptotic distribution for the local polynomial estimators, with their diagonal elements standardized by their standard errors $(\hat V_{f,p+1}(h_f))^{1/2}$ and $(\hat V_{Z_k,l+1}(h_k))^{1/2}$. We construct the standard errors of $\tau_{Z_k}$ as $\hat V_{Z_k,l+1}(h_k)\stackrel{P}{\to} V_{Z_k,l+1}$, where $V_{Z_k,l+1}$ is the asymptotic variance of $\sqrt{n h_{k}}(\hat\tau_{Z_k,l+1}(h_k)-\tau_{Z_k})$. The asymptotic covariance matrix of the standardized local polynomial estimators is the asymptotic correlation matrix of the local polynomial estimators for $\tau_f$ and $\tau_{Z_k}$. The correlation matrix $V^*_{l+1,p+1}$ is \[ V^*_{l+1,p+1} =

, \] where $V^*_{Z,l}=\operatorname{diag}(V_{Z_1,l}^{-1/2},\ldots,V_{Z_d,l}^{-1/2})V_{Z,l}\operatorname{diag}(V_{Z_1,l}^{-1/2},\ldots,V_{Z_d,l}^{-1/2})$; our estimator of $V^*_{l+1,p+1}$ is \[ \hat V^*_{l+1,p+1}(h) =

, \] where

In the following proposition, we show that the finite-sample distribution of $\tilde \chi^2_{l+1,p+1}(h)$ is approximated by the distribution of the sum of the squares of the $(d+1)$-dimensional multivariate normal random vector.

Proposition (ref) establishes the asymptotic validity and consistency of the $\alpha$-level testing procedure, which rejects $H_0$ if $\hat c_{l+1,p+1}(\alpha)$ exceeds $\tilde \chi^2_{l+1,p+1}(h)$. The critical value $\hat c_{l+1,p+1}(\alpha)$ is obtained numerically by generating $(d+1)$-dimensional standard normal random vectors using Monte Carlo simulation.

The Max test statistic with $L^{\infty}$ norm takes the following form.

Similar to the sWald statistic, the Max statistic is constructed by standardizing diagonal elements. The following statement shows that the finite sample distribution of $\hat M_{l+1,p+1}(h)$ is approximated by the distribution of the maximum value of the $(d+1)$-dimensional multivariate normal random vector. \footnote{This construction via simulation is similar to the plug-in sup-t implementation in oleaSimultaneousConfidenceBands2019.}

Proposition (ref) establishes the asymptotic validity and consistency of the $\alpha$-level testing procedure, which rejects $H_0$ if $\hat M_{l+1,p+1}(h)$ exceeds $\hat m_{l+1,p+1}(\alpha)$. The critical value $\hat m_{l+1,p+1}(\alpha)$ is obtained numerically by generating $(d+1)$-dimensional standard normal random vectors using Monte Carlo simulation.

Simulation evidences for the unified tests

Given the established theoretical properties of the proposed tests, their performance was evaluated using the Monte Carlo experiment. We conducted $3000$ replications to generate a random sample $\{(X_i, Z_{1,i}, \ldots, Z_{d,i})':i=1,\ldots,n\}$ with size $n= \{500, 1000\}$ for each of them. \footnote{The detailed simulation data generating process is in Online Appendix C.}

We specify the distribution of running variable $X_i$ as the weighted average of two truncated normal distributions, where $\bar{p} \geq 0.5$ determines the weights. The density of $X_i$ is shown in Figure (ref). The density $f(x)$ is continuous at $x=0$ if and only if $\bar{p} = 0.5$ but jumps at $x=0$ when $\bar{p} > 0.5$.

For covariates $Z_i$, we consider two specifications: first, one of the $d$ covariates sees a jump; second, all $d$ covariates see a jump, with each jump size divided by $d$. For the first specification, we have $\mu_{Z_k}(x)=E[Z_{k,i}|X_i=x], k = 1\ldots, d$ as follows. \[ \mu_{Z_k}(x) = \lambda(x) for k=1,\ldots,d-1\quad and \quad \mu_{Z_d}(x) =

\] where only the $d$th covariate sees a jump, and the functional form of $\lambda(x)$ is obtained from the simulated data close to the data in leeRandomizedExperimentsNonrandom2008 following calonicoRobustNonparametricConfidence2014a, where $a\geq0$ is the level of discontinuity of $\mu_{Z_d}(x)$ at $x=0$. For the second specification, we have $\mu_{Z_k}(x)=E[Z_{k,i}|X_i=x], k=1,\ldots,d$ as follows. \[ \mu_{Z_k}(x) =

\] where $a/d$ denotes the discontinuity level for each $\mu_{Z_k}(x)$ at $x=0$, For both cases, let $1>\rho\geq0$ denote the conditional correlation between $Z_{j,i}$ and $Z_{k,i}$, given $X_i=x$. In this data-generating process, the null hypothesis for the joint manipulation test, $H_0: (\tau_f,\tau_{Z_1},\ldots,\tau_{Z_d})=0$ is true if and only if $\bar{p} = 0.5$ and $a = 0$; $\tau_f > 0$ when $\bar{p} > 0.5$; $\tau_{Z_d}=a>0$ when $a>0$.

We conducted the tests with size $\alpha=0.05$, employing local polynomial estimators $(\hat{\tau}_{Z_{1},2}(h_1),….,\hat{\tau}_{Z_{d},2}(h_d), \hat{\tau}_{f,3}(h_f))$ described in Section (ref) with $K(u)=(1-|u|)\vee0$, where the bandwidths are chosen to be MSE-optimal for $\hat\tau_{Z_k,1}(h_k)$ and $\hat\tau_{f,2}(h_f)$. The selection of the kernel function and polynomial order follows calonicoRobustNonparametricConfidence2014a and cattaneoSimpleLocalPolynomial2019. \footnote{As mentioned in the previous sections, the polynomial orders are chosen so that the estimates are bias-corrected for their valid inference. This choice of polynomial order is the default for rddensity and available as $\rho = 1$ option for rdrobust.} Below, we report the empirical sizes and powers of the conventional and proposed testing methods.

We compared five testing methods: (i) naive testing; (ii) Bonferroni correction of naive testing; (iii) Wald (chi-square) test; (iv) Max test proposed in Section (ref); (v) sWald test proposed in Section (ref). Naive testing separately conducts $d+1$ asymptotic $t$ tests, as described in calonicoRobustNonparametricConfidence2014a and cattaneoSimpleLocalPolynomial2019 with a size $\alpha$ for each of the $d+1$ restrictions of the null hypothesis, $\tau_f=0, \tau_{Z_1}=0,\ldots,\tau_{Z_d}=0$. Naive testing rejects the null hypothesis if at least one test rejects it. The Bonferroni correction conducts asymptotic $t$ tests separately, but the size is changed to $\alpha/(d+1)$ to restrict type 1 error of the null hypothesis that all $d+1$ restrictions hold. \footnote{We do not benefit from rejecting more than one test for evaluating the underlying identifying restriction. Thus, Holm correction for Bonferroni test does not improve the following results because the probability of rejecting at least one null hypothesis remains the same in two corrections.}

In Tables (ref)–(ref), we compare the empirical size of each test under the null hypothesis (where $\bar{p} = 0.5$ and $a = 0$) with different correlation coefficients $\rho=0,0.5,0.9$. Notably, naive testing suffers from size distortion, and over-rejection worsens when $d$ is large. Similarly, the Wald test without correction suffers from size distortion as $d$ increases. As discussed in Section (ref), the effective sample size may be too small compared to a large $d$. When $d$ is large, the covariance matrix increases, and inverse matrix estimation becomes unstable when the effective sample size is limited.

The performances of the other methods vary according to sample size, dimension of covariates, and correlation coefficients. For the relatively small sample size of $n = 500$, the superiority of these methods is ambiguous. For a relatively small number of covariates $dim \in \{1,3,5\}$, the sWald test controls for size by $6\%$, whereas the Max test and Bonferroni correction control for size by approximately $7\%$. For a larger number of covariates, the Bonferroni correction and Max test may fail to control size by $10\%$, whereas the sWald test controls size by $7\%$ throughout. For a larger sample size of $n = 1000$, the sWald test dominates the Bonferroni correction for size control. The sWald test controls size by $5.5\%$ for any correlation and number of covariates up to $25$. \footnote{We consider that the failure in the Bonferroni correction comes from challenges in evaluating the tail of the nonparametric test statistics.} For a moderate number of covariates up to $dim = 5$ or $10$, the Max test can control the size by $6\%$ and is not substantially worse than the Bonferroni correction. Therefore, when the sample size is not small, the sWald test can control the size for any dimension up to dim = $25$; the Max test and Bonferroni correction can control the size up to dim = $5$ or $10$.

Given the results of size control, we evaluate their empirical power properties. We compare the empirical power of the (ii) Bonferroni correction, (iv) Max test, and (v) sWald test under $\bar{p}=0.575$ and $\tau_{Z_d}(=a)=0,0.5,1,1.5,2$ for different dimensions $d=1,3,5$. Figure (ref) presents the results for the scalar covariate. Figures (ref), and (ref) present the results when one of the $d$ covariates sees a jump. Further, Figures (ref) and (ref) show the results when all $d$ covariates see a jump. However, each jump size is divided by $d$ for $d=3,5$. In Figures (ref)–(ref), each figure presents the results for different correlations $\rho = 0.5, 0.9$.\footnote{In Online Appendix D , we present the results for weaker correlations and results containing negative correlations.}

In all cases, the Max test exhibits higher power than the Bonferroni correction. When one of the covariates jumps (Figures (ref) and (ref)), the Max test outperforms the Bonferroni correction and sWald test, and the gap among the three methods expands as the correlation and dimension increase. While the sWald test has relatively less power to detect small jumps in the covariate, the power of the sWald test catches up quickly as the jump size increases.

When all covariates have seen jumps (see Figure (ref) and (ref)), the sWald test outperforms the Max test and Bonferroni correction in most cases. When the correlation is weak, the sWald test is substantially superior to the other two methods; however, the Max test approaches the sWald test as the correlation and dimensions increase. In Online Appendix D , we report qualitatively similar simulation results for other specifications, including lower or negative correlation coefficients.

To summarize, we recommend the Max test when one suspects that a few covariates exhibit jumps but not the others among a moderate number of covariates. We recommend the sWald test for other cases because it has superior size control and power properties for detecting jumps among many covariates.

Extension

Our unified test solves the over-rejection problem prevalent in empirical studies. Below, we consider two extensions of our analysis to supplement our unified test.

Equivalence Testing

The diagnostic test evaluates the null hypothesis that suggests the plausibility of the underlying identification. Hence, the diagnostic procedure tests whether the desired null hypothesis can be rejected, and we control the probability of falsely rejecting the desired null hypothesis of no imbalance or discontinuity. In other words, the roles of type I and type II errors are flipped.

Several studies propose equivalence or non-inferiority testing as alternatives to flipped testing for the desirable null hypothesis. See, for example, wellek2010equivalence and hartman2018equivalence for further details. In RD designs, hartman2021equivalence propose equivalence tests for diagnostic testing of each covariate and density. As before, the existing equivalence tests are uni-variate, and a unified alternative is necessary. Hence, we developed a unified equivalent test for multivariate balance tests and a density test in terms of the max test based on wellek2010equivalence with a rectangular equivalence domain.

Recall that we want to verify whether the following $d+1$ restrictions hold simultaneously \[ f_+ = f_-,\quad \mu_{Z_1+} = \mu_{Z_1-}, \quad \ldots, \quad \mu_{Z_d+} = \mu_{Z_d-}. \] An appropriate set of hypotheses for the equivalence test is as follows:

where $\varepsilon > 0$ determines the range within which the imbalance or discontinuity is sufficiently small.

To conduct an equivalence test, the researcher must select an appropriate equivalence range $(-\varepsilon,\varepsilon)$. We follow hartman2018equivalence and hartman2021equivalence to focus on the equivalence confidence interval, which is the shortest equivalence range $(-\varepsilon,\varepsilon)$ such that the alternative hypothesis with the equivalence range would not be detected at the prespecified level $\alpha$ for the observed data. Hence, we achieve an equivalence confidence interval without manually selecting $\epsilon$.

Denote \[

\] and \[

\] The following proposition demonstrates the asymptotic validity of the proposed equivalence test for a prespecified level $\alpha$.

Pretesting Analysis

Following roth2022pretest, we analyze the bias of the treatment effect estimator induced by conducting a pretest. Let $\hat{\tau}_{Y}, \hat{\tau}_{Z}, \hat{\tau}_{f}$ be the local polynomial estimators of $\lim_{x \to 0+} E[Y_i|X_i=x] - \lim_{x \to 0-} E[Y_i|X_i=x], (\mu_{Z_1+}-\mu_{Z_1-}, \ldots, \mu_{Z_d+}-\mu_{Z_d-})', f_{+}-f_{-}$, respectively.\footnote{To apply Theorem (ref), we use the bias-corrected estimators as $\hat{\tau}_{Y}, \hat{\tau}_{Z}, \hat{\tau}_{f}$.} For simplicity, suppose that $d=1$ and $(\hat{\tau}_{Y}, \hat{\tau}_{Z}, \hat{\tau}_{f})'$ is multivariate normally distributed; that is, \[ \left(

\right) \ \sim \ N \left( \left(

\right) + \left(

\right), \left(

\right) \right), \] where $\tau_{f}$ is the treatment effect and $\delta = (\delta_Y, \delta_Z, \delta_f)'$ is the unconditional bias from the imbalance at the cutoff. Because Theorem (ref) implies that the asymptotic covariance between $(\hat{\tau}_Y, \hat{\tau}_Z)'$ and $\hat{\tau}_f$ is zero, we assume that $\hat{\tau}_f$ is independent of $\hat{\tau}_Y$ and $\hat{\tau}_Z$.

Let $\mathcal{T}(\Sigma) \subset \mathbb{R}^{2}$ be the acceptance region, which potentially depends on the covariance matrix $\Sigma = (\sigma_{ZZ},\sigma_{ff})$. From Proposition 1 in roth2022pretest, the expectation of $\hat{\tau}_Y$ conditional on passing the pretest can be written as

Equation ((ref)) implies that the expectation of $\hat{\tau}_{Y}$ conditional on passing the pretest is the sum of the treatment effect $\tau_{Y}$, the unconditional bias $\delta_{Y}$, and the pretest bias $\frac{\sigma_{YZ}}{\sigma_{ZZ}} \left\{ E\left[ \hat{\tau}_{Z} | (\hat{\tau}_Z, \hat{\tau}_{f})' \in \mathcal{T}(\Sigma) \right] - \delta_{Z} \right\}$. When considering the max test, the acceptance region can be written as

where $c_Z$ and $c_f$ are positive constants that depend on $\Sigma$. Because $\hat{\tau}_f$ is independent of $\hat{\tau}_Z$, ((ref)) becomes

Hence, when we use the Max test as the pretest, the density test statistic $\hat{\tau}_{f}$ does not affect the pretest bias, which disappears if $\delta_Z=0$.

To discuss the sign of the pretest bias, we consider the following model:

where $\mu_Y$ is a continuous function, and $\epsilon_i$ is an error term satisfying $E[\epsilon_i|X_i,Z_i]=0$. Then, we have \[ \lim_{x\to0^+} E[Y_i|X_i=x] - \lim_{x\to0^+} E[Y_i|X_i=x] \ = \ \tau_Y + \beta \delta_Z, \] where we use $\lim_{x\to0^+} E[Z_i|X_i=x] - \lim_{x\to0^+} E[Z_i|X_i=x] = \delta_Z$. This implies that the unconditional bias becomes $\delta_Y = \beta \delta_Z$. Hence, if $\beta > 0$ and $\sigma_{YZ} >0$, then the sign of the pretest bias is opposite to that of the unconditional bias $\delta_Y$. Therefore, in this case, the pretest can reduce bias. In Online Appendix E, we conduct a numerical analysis using simulated and calibrated datasets. The numerical results also support our theoretical analysis that pretesting may reduce bias.

An analogous analysis for the variance implies that the variance of the treatment effect estimate is likely to be smaller after pretesting. Similar to Proposition 3 of roth2022pretest, we obtain

If the max test is conducted, the acceptance region is convex. Hence, Proposition 4 in roth2022pretest implies that \[ V\left[ \hat{\tau}_{Z} | (\hat{\tau}_Z, \hat{\tau}_{f})' \in \mathcal{T}(\Sigma) \right] - \sigma_{ZZ} \leq 0. \] Hence, the variance is likely to be smaller after pretesting. \footnote{The above analysis indicates that pretesting can reduce the bias and variance of the treatment effect estimates. Its implication for the power of the test for the treatment effect is indefinite. For example, pretesting reduces the variance of estimator when the bias is zero $\delta_Y = \delta_Z = 0$. The smaller variance may increase the power when $\tau_Y$ is away from $0$ and decrease the power when $\tau_Y$ is close to $0$.} \FloatBarrier

Conclusion

Diagnostic tests are vital for applied studies of RD designs; however, their inappropriate use frustrates their purpose. The diagnostic procedure for RD design aims to test the underlying null hypothesis that all pre-determined covariates balance and the density function is continuous. Current practices test these restrictions separately without appropriate size control---they over-reject the null hypothesis that the identifying restriction for the RD design is plausible.

We demonstrate an extraordinary over-rejection of the null hypothesis, where each testable restriction is valid. Among the $59$ studies published in the top five economics journals, $32\%$ ($19$ studies) rejected at least one testable restriction. Nevertheless, less than $5\%$ of the $787$ tests were rejected. Hence, we conclude that the plausibility of identification is highly over-rejected where the underlying null hypotheses appear to hold. We blame the multiple testing problem for the over-rejection because most studies run more than $10$ covariates balance tests separately.

We provide a solution to this over-rejection problem. We demonstrate the joint asymptotic normality of RD design estimators and provide alternative unified testing procedures. We show two different unified tests: the standardized Wald (sWald) test and the Max test. In the numerical simulations, the sWald and Max statistics achieved appropriate size-control properties for a moderate number of covariates when the sample size was $1000$. The sWald test is particularly superior to other methods in size control, as it can control for any dimension of covariates up to $25$ where the Bonferroni correction may fail to control the size. The power properties vary according to the underlying alternative hypotheses. The Max test is superior to the Bonferroni correction when only one covariate exhibits a jump in the conditional mean. The sWald test is superior to the Bonferroni correction when all the covariates jump in their conditional means.

In conclusion, we recommend the Max test for a design with a moderate number of covariates when only a few covariates which critically matter for the plausibility of identification such as placebo outcomes. We recommend the sWald test for other designs. This test is appropriate for most practical designs where there is no ex-ante concerning covariate and testing many covariate balancing is the norm.

Our unified test provides a formal judgment for whether the rejected diagnostic tests were spurious and by luck or not. In RD designs, rejections in the diagnostic tests can be genuine signals for their designs being implausible while they may be spurious. Hence, our procedure offers a transparent decision rule that is essential to the users if the diagnostic tests are part of the assessment in RD designs.

Furthermore, our procedure is complementary with the other current practices. For example, practitioners often conducts an additional inspection for a covariate with an apparent graphical evidence for a jump. We continue to recommend graphical analyses and those careful inspections because such procedure is complementary with our unified tests. With our unified test, practitioners may have stronger argument for a plausible design in the following combination of arguments: our unified test suggests that the rejections were spurious; furthermore, other anecdotes conclude that the spurious rejections were immaterial to the plausibility of the design. Hence, our procedure shall be an essential support for effectively justifying designs with a few rejected diagnostic tests when the underlying design is genuinely valid.

Based on our findings, we consider several topics for future research. First, finding another use for testable restrictions may be preferable. To validate the plausibility of identification, we must focus on the underlying null hypothesis rather than detecting which restrictions fail. Nevertheless, if we detect a particular subgroup of covariates that threats the plausibility of identification, such a finding may provide valuable information. Diagnostic procedures using such information may improve design validation. Second, although we focused on the statistical inference of testable restrictions, a graphical analysis of testable restrictions is also important. Recent studies, such as Calonico_Cattaneo_Titiunik_2015_JASA, have considered the optimal procedure for visualization, and Korting_Lieberman_Matsudaira_Pei_Shen_2023 have demonstrated that tuning parameters in visualization is critical. We recommend our unified testing along with graphical analysis, with multiple testing problems in graphical analyses as a future research topic. Finally, evaluating the joint distribution of testing statistics that are nonparametric estimates is challenging. We will consider these issues as future research questions.

\FloatBarrier

Acknowledgements

The study was supported by JSPS KAKENHI Grant Numbers JP20J20046, JP23K18799 (Fusejima), JP22K13373 (Ishihara), and JP21K13269 (Sawada). This study is derived from a chapter of Koki Fusejima's dissertation. We thank Eiji Kurozumi, Yukitoshi Matsushita, Katsumi Shimotsu, and Kohei Yata, the anonymous referees and associate editor as well as seminar participants at Hitotsubashi University, The University of Tokyo, Kansai Keiryo Keizaigaku Kenkyukai, Japanese Joint Statistical Meeting, International Conference on Econometrics and Statistics, and the Asian Meeting of the Econometric Society in East and South-East Asia. Moreover, we thank Masaki Oguni and Yuri Sugiyama for their excellent research assistance.