Introduction

The standard regression discontinuity (RD) design assigns treatment based on whether a single running variable X crosses a known threshold c: units with X ≥ c are treated; those with X < c are not [Imbens and Lemieux, 2008]. The elegance of this design lies in the local randomisation induced near the threshold: units just above and just below c are similar in all ways except treatment status, allowing causal identification under continuity assumptions.  

Many real-world policies, however, assign treatment based on multiple eligibility criteria simultaneously. Students may be eligible for remediation if they fail both a reading test and a maths test. Workers may qualify for a benefit if their income falls and their age exceeds a threshold. Applicants may be admitted to a selective programme based on a composite of multiple scores. When treatment is assigned on the basis of multiple running variables, the single-score RD framework does not apply directly.  

This article surveys the identification and estimation challenges posed by multiple running variable (MRV) designs, introduces the main approaches in the literature, and describes available software.  

1 Two Canonical Settings

1.1 The Bivariate Frontier

Suppose treatment is assigned based on two continuous running variables X₁ and X₂, both observed without manipulation. The simplest case is a bivariate threshold design where the treatment region is:  

creating a rectangular treatment region in the (X₁, X₂) plane. The treatment frontier is the boundary of the treatment region— here, two line segments meeting at the corner point (c₁, c₂).  

An empirically important example is studied by Papay et al. [2011]: students in many U.S. states must pass reading and maths tests to avoid mandatory remediation. A student avoids remediation only by scoring above both thresholds. The relevant frontier is two perpendicular lines in score space.  

1.2 The Single-Index Frontier

A second common structure is a single-index or "combined score" threshold, where treatment is assigned based on a weighted sum of two scores:

This creates a linear frontier in the (X₁, X₂) plane. College admissions, public housing eligibility, and welfare benefit cutoffs sometimes take this form.  

2 The Identification Problem

2.1 Multiple Treatment Effects Along the Frontier

The key challenge in MRV designs is that the treatment effect is not a single number but a function defined along the frontier. Consider the bivariate threshold design. At a point on the frontier where X₁ = c₁ and X₂ > c₂, treatment is assigned by crossing the X₁ threshold while holding X₂ fixed. The treatment effect identified at this frontier point is the effect for units with X₁ ≈ c₁ and a specific value of X₂.  

Formally, let τ(x₁, x₂) denote the treatment effect conditional on running variable values (x₁, x₂). The MRV RD identifies at points along the frontier, not at a single point. The shape of τ(·) along the frontier can reveal heterogeneity— for example, whether the effect of failing a reading test is larger for students who are also close to the maths threshold.  

Reardon and Robinson [2012] formalise this by distinguishing:  

• The frontier average treatment effect (FATE): the average of τ(x₁, x₂) over the frontier, weighted by the density of observations on the frontier.  

• The corner treatment effect: the effect at the corner point (c₁, c₂) where both conditions bind simultaneously.  

2.2 The Corner Identification Problem

At the corner point, the density of observations may be very thin: very few units have X₁ ≈ c₁ and X₂ ≈ c₂ simultaneously. This creates an estimation challenge because identification at the corner requires extrapolation from frontier segments far from the corner. The degree of extrapolation required depends on the joint density of (X₁, X₂) near the corner— if the joint density is low there, estimates will be noisy and extrapolation-dependent.  

3 The Papay et al. [2011] Approach

Papay et al. [2011] propose a practical solution: estimate the treatment effect separately at multiple points along each frontier segment and then aggregate. Specifically, along the segment where X₁ = c₁ (the vertical frontier), they estimate:  ‍


within a bandwidth around X₁ = c₁, conditioning on X₂ᵢ. This is a standard polynomial RD applied locally to the running variable that triggers the threshold, controlling flexibly for the other running variable.  

The treatment effect τ̂₁(x₂) varies as a function of X₂. Averaging over the empirical distribution of X₂ values among frontier-adjacent units yields the frontier-average treatment effect along this segment.  

4 Centering and the Multivariate Running Variable

An alternative approach, applicable when the two running variables are conceptually interchangeable, is to construct a distance to the frontier as a single composite running variable. Let:  

Treatment is assigned when Zᵢ ≥ 0. This transforms the MRV design into a standard single-score RD in Z.  

The challenge is measuring the signed distance. For a rectangular frontier, one natural choice is:  

which equals the minimum margin of eligibility across the two criteria. Units with Zᵢ < 0 fail at least one criterion; units with Zᵢ ≥ 0 pass both. This "centering" approach allows the analyst to use standard RD estimators and bandwidth selectors.  

Keele and Titiunik [2015] apply a version of this idea in geographic RD designs, where the boundary is a county line or legislative district boundary rather than a pair of score thresholds. They show that defining distance to the boundary as a single running variable and applying standard RD methods can recover treatment effects under conditions analogous to the single-score continuity assumption.  

5 The rdmulti Package

Cattaneo et al. [2020] develop the rdmulti package for Stata and R, which implements RD analysis in settings with multiple cutoffs or multiple running variables. For the multiple running variable setting, the package provides:  

1. Estimation of frontier-average treatment effects using local polynomial methods applied along each segment of the treatment frontier.  
2. Optimal bandwidth selection adapted to the bivariate problem.  
3. Bias-corrected robust confidence intervals extending Calonico et al. [2014].  
4. Tools for plotting the estimated treatment effect surface along the frontier.  

The package requires the user to specify the frontier structure and the local polynomial order and bandwidth for each running variable.  

6 Empirical Application: Remediation Thresholds

Papay et al. [2011] apply their framework to data from a large U.S. school district where students must pass both a reading and a maths test to avoid mandatory remediation in 9th grade. Treatment (avoiding remediation) is assigned iff X₁ ≥ c₁ (reading score above threshold) and X₂ ≥ c₂ (maths score above threshold).  

Their key findings:  

• Along the reading frontier (students who just passed reading but varied in maths): avoiding remediation improved high school GPA by 0.08 grade points.  

• Along the maths frontier (students who just passed maths but varied in reading): avoiding remediation improved GPA by 0.12 grade points.  

• The effect at the corner point (students who barely passed both) was approximately 0.10, consistent with a convex combination of the segment effects.  ‍


7 Assumptions and Limitations

7.1 Continuity

As in the single-score RD, identification requires that potential outcomes are continuous in (X₁, X₂) at the frontier. Manipulation of either score would invalidate this assumption; McCrary-style density tests [McCrary, 2008] should be conducted for each running variable separately.  

7.2 No Simultaneous Manipulation

An additional concern unique to MRV designs is that units might manipulate one score in anticipation of how it interacts with the other. For example, a student who knows they score just below the reading threshold might not bother preparing for the maths test, creating a discontinuous relationship between maths score and unobservables at the reading threshold. This type of strategic interaction is difficult to detect and requires institutional knowledge to evaluate.  

7.3 Limited Power at the Corner

As discussed, the corner treatment effect may be poorly identified if the joint density of (X₁, X₂) near the corner is low. Researchers should report the empirical density near the frontier and assess whether corner estimates rely on substantial extrapolation.  

8 Conclusion

Multiple running variable RD designs are increasingly common as researchers study multi-criteria eligibility rules in education, social policy, and finance. The key insight from Papay et al. [2011] and Reardon and Robinson [2012] is that MRV designs identify a frontier function of treatment effects rather than a single point estimate. The rdmulti package provides practical implementation tools. The main limitations are thin data at frontier corners, the potential for cross-variable strategic manipulation, and the difficulty of summarising a multidimensional treatment effect surface in a single interpretable number.  

References

1. Calonico, S., Cattaneo, M. D., and Titiunik, R. (2014). Robust nonparametric confidence intervals for regression-discontinuity designs. Econometrica, 82(6), 2295-2326.  
2. Cattaneo, M. D., Titiunik, R., and Vazquez-Bare, G. (2020). The regression discontinuity design. In V. Chernozhukov, C. Hansen, N. Kallus, M. Spindler, and V. Syrgkanis (Eds.), Handbook of Causal Methods. MIT Press.  
3. Imbens, G. W. and Lemieux, T. (2008). Regression discontinuity designs: A guide to practice. Journal of Econometrics, 142(2), 615-635.  
4. Keele, L. and Titiunik, R. (2015). Geographic boundaries as regression discontinuities. Political Analysis, 23(1), 127-155.  
5. McCrary, J. (2008). Manipulation of the running variable in the regression discontinuity design: A density test. Journal of Econometrics, 142(2), 698-714.  
6. Papay, J. P., Willett, J. B., and Murnane, R. J. (2011). Extending the regression-discontinuity approach to multiple assignment variables. Journal of Econometrics, 161(2), 203-207.  
7. Reardon, S. F. and Robinson, J. P. (2012). Regression discontinuity designs with multiple rating-score variables. Journal of Research on Educational Effectiveness, 5(1), 83-104.