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.
111,967 characters · 15 sections · 42 citation commands
Algorithm as Experiment: Machine Learning, Market Design, and Policy Eligibility Rules
Today's society increasingly resorts to algorithms for decision making and resource allocation. For example, judges in the US make legal decisions aided by predictions from supervised machine learning algorithms. Supervised learning is also used by governments to detect potential criminals and terrorists, and by banks and insurance companies to screen potential customers. Tech companies like Facebook, Microsoft, and Netflix allocate digital content by reinforcement learning and bandit algorithms. Retailers and e-commerce platforms engage in algorithmic pricing. Similar algorithms are encroaching on high-stakes settings, such as in education, healthcare, and the military.
Other types of algorithms also loom large. School districts, college admissions systems, and labor markets use matching algorithms for position and seat allocations. Objects worth astronomical sums of money change hands every day in algorithmically run auctions. Many public policy domains like Medicaid often use algorithmic rules to decide who is eligible.
All of the above examples share a common trait: a decision-making algorithm makes decisions based only on its observable input variables. Thus conditional on the observable variables, algorithmic treatment decisions are assigned independently of any potential outcome. This property turns algorithm-based treatment decisions into instrumental variables (IVs) that can be used for measuring the causal effect of the final treatment assignment. The algorithm-based IV may produce stratified randomization, regression-discontinuity-style local variation, or some combination of the two.
This paper shows how to use data obtained from algorithmic decision making to identify and estimate causal effects. In our framework, the analyst observes a random iid sample $\{(Y_i,X_i,D_i,Z_i)\}_{i=1}^n$, where $Y_i$ is the outcome of interest, $X_i\in\mathbb{R}^p$ is a vector of pre-treatment covariates used as the algorithm's input variables, $D_i$ is the binary treatment assignment, possibly made by humans, and $Z_i$ is the binary treatment recommendation made by a known algorithm. The algorithm takes $X_i$ as input and computes the probability of the treatment recommendation $A(X_i)=\Pr(Z_i=1|X_i)$. $Z_i$ is then randomly determined based on the known probability $A(X_i)$ independently of everything else conditional on $X_i$. The algorithm's recommendation $Z_i$ may influence the final treatment assignment $D_i$, determined as $D_i=Z_iD_i(1)+(1-Z_i)D_i(0)$, where $D_i(z)$ is the potential treatment assignment that would be realized if $Z_i=z$. Finally, the observed outcome $Y_i$ is determined as $Y_i=D_iY_i(1)+(1-D_i)Y_i(0)$, where $Y_i(1)$ and $Y_i(0)$ are potential outcomes that would be realized if the individual were treated and not treated, respectively. This setup is an IV model where the IV satisfies the conditional independence condition but may not satisfy the overlap (full-support) condition. This setup nests the classic propensity-score and regression-discontinuity-design (RDD) setups.
Within this framework, we first characterize the sources of causal-effect identification for a class of data-generating algorithms. This class includes all of the aforementioned examples, nesting both stochastic and deterministic algorithms. The sources of causal-effect identification turn out to be summarized by a suitable modification of the Propensity Score. We call it the Approximate Propensity Score (APS). For each covariate value $x$, the Approximate Propensity Score is the average probability of a treatment recommendation in a shrinking neighborhood around $x$, defined as $$ p^{A}(x)\equiv \lim_{\delta\rightarrow 0}\frac{\int_{B(x,\delta)}A(x^*)dx^*}{\int_{B(x,\delta)}dx^*}, $$ where $B(x,\delta)$ is a $p$-dimensional ball with radius $\delta$ centered at $x$. The Approximate Propensity Score provides an easy-to-check condition for what causal effects the data from an algorithm allow us to identify. In particular, we show that the conditional local average treatment effect (LATE; imbens1994) at covariate value $x$ is identified if and only if the Approximate Propensity Score is nondegenerate, i.e., $p^{A}(x)\in (0, 1)$.
The identification analysis suggests an estimator. The treatment effects can be estimated by two-stage least squares (2SLS) where we regress the outcome on the treatment with the algorithm's recommendation as an IV. To make the algorithmic recommendation a conditionally independent IV, we propose to control for the Approximate Propensity Score, as formalized below.\footnote{Code implementing this procedure in Python, R, and Stata is available at \url{https://github.com/rfgong/IVaps}.}
As the main theoretical result, we prove the 2SLS estimator $\hat\beta_1^s$ is a consistent and asymptotically normal estimator of a well-defined causal effect (weighted average of conditional local average treatment effects). Our result clarifies how to estimate other parameters by reweighting observations. We also show that inference based on the conventional 2SLS heteroskedasticity-robust standard errors is asymptotically valid as long as the bandwidth $\delta$ goes to zero and the number of simulation draws $S$ goes to infinity at appropriate rates. We prove the asymptotic properties by exploiting results from differential geometry and geometric measure theory.
Our estimator is applicable even if the algorithm is deterministic and produces multidimensional regression-discontinuity variation only. In contrast to standard multidimensional RDD methods that define the distance to the nearest boundary as a single running variable, our estimator is straightforward to use even when the multidimensional boundary is arbitrarily complex and the distance is hard to compute. Moreover, our method applies to more general settings with stochastic algorithms, deterministic algorithms, and combinations of the two.
In settings that mix stochastic and deterministic algorithms, our estimator exploits both the random-assignment variation and RDD variation, producing precision and representability gains compared to RDD estimators and propensity-score estimators only using either variation.
The practical performance of our estimator is demonstrated through simulation and an original application. We first conduct a Monte Carlo simulation mimicking real-world decision making based on machine learning algorithms. We consider a data-generating process combining stochastic and deterministic algorithms. Treatment recommendations are randomly assigned for a small experimental segment of the population and are determined by a high-dimensional, deterministic machine learning algorithm for the rest of the population. Such combinations of small experiments and deterministic treatment allocations arise in the real world when large-scale experiments are prohibited due to ethical, budget, or legislative constraints. Our estimator is shown to be feasible in this high-dimensional setting and has smaller median absolute errors relative to alternative estimators. In particular, our method produces more efficient and representative estimates than a conventional propensity-score approach. The simulation also illustrates a setting where the RDD boundary is complex and no prior multidimensional RDD method is applicable.
Our empirical application is an analysis of COVID-19 hospital relief funding. The Coronavirus Aid, Relief, and Economic Security (CARES) Act designated \$175 billion for COVID-19 response efforts and reimbursement to health care entities for expenses or lost revenues kakani2020allocation. This policy intended to help hospitals hit hard by the pandemic, as “financially insecure hospitals may be less capable of investing in COVID-19 response efforts” khullar2020COVID. We ask whether this problem is alleviated by the relief funding for hospitals.
We identify the causal effects of the relief funding by exploiting the funding eligibility rule. The government runs an algorithmic rule on hospital characteristics to decide which hospitals are eligible for funding. This fact allows us to apply our method to estimate the effect of relief funding. Specifically, our 2SLS estimators use funding eligibility status as an IV for funding amounts, while controlling for the Approximate Propensity Score induced by the eligibility-determining algorithm. The funding eligibility IV boosts the funding amount by about \$15 million on average.
The resulting 2SLS estimates with Approximate Propensity Score controls suggest that COVID-19 relief funding has little to no effect on outcomes, such as the number of COVID-19 patients hospitalized at each hospital. The estimated causal effects of relief funding are much smaller and less significant than the naive ordinary least squares (OLS) (with and without controlling for hospital characteristics) or 2SLS estimates with no controls. The OLS estimates, for example, imply that a \$1 million increase in funding allows hospitals to accommodate 4.53 more COVID-19 patients. The uncontrolled 2SLS estimates produce similar, slightly smaller effects (2.44 more patients per \$1 million of funding). In contrast, the 2SLS estimates with Approximate Propensity Score controls show no or even negative effects (up to 2.21 fewer patients for every \$1 million of funding).
The null effect of funding persists several months after the distribution of funding. We also find no clear heterogeneity in the null funding effect across different subgroups of hospitals. Our finding provides causal evidence for the concern that funding in the CARES Act might not have been well targeted to the clinics and hospitals with the greatest needs.\footnote{See, for example, kakani2020allocation as well as {\it Forbes}'s article, “Hospital Giant HCA To Return \$6 Billion in CARES Act Money,” at \url{https://www.forbes.com/sites/brucejapsen/2020/10/08/hospital-giant-hca-to-return-6-billion-in-cares-act-money}.}
Our framework integrates the classic propensity-score (selection-on-observables) scenario with a multidimensional extension of the fuzzy RDD. We analyze this integrated setup in the IV world with noncompliance. This general setting appears to have no prior established estimator. armstrong2020finitesample provide an estimator for a related setting with perfect compliance.
When we specialize our estimator to the multidimensional RDD case, our estimator has three features. First, it is a consistent and asymptotically normal estimator of a well-interpreted causal effect (average of conditional treatment effects along the RDD boundary) even if treatment effects are heterogeneous. Second, it uses observations near all the boundary points as opposed to using only observations near one specific boundary point, thus avoiding variance explosion even when $X_i$ has many elements. Third, it can be easily implemented even in cases with many covariates and complex algorithms (RDD boundaries). No existing estimator appears to have all of these properties papay2011extending,zajonc2012regression,keele2015geographic,cattaneo2016interpreting,imbens2019optimized.
A popular approach to the two-dimensional RDD is to use the shortest (Euclidean) distance from each individual to the boundary as a univariate running variable and apply a univariate RDD method black1999border.\footnote{Another common approach to the two-dimensional RDD is to first estimate the conditional average treatment effect $E[Y_i(1)- Y_i(0)|X_i = x]$ for a large number of boundary points $x$ (either by the univariate local polynomial regression using the distance to the point $x$ as a univariate covariate or by the bivariate local polynomial regression). It then computes a weighted average of the estimated conditional average treatment effects over the boundary zajonc2012regression,keele2015geographic. However, identifying boundary points from a general decision algorithm is hard unless it has a known analytical form. Even if we can trace out the boundary, it is not straightforward to select a grid of points along the boundary.} In general multidimensional RDDs where the boundary is complex or its analytical form is unknown, the distance-based approach requires approximation or estimation of the shortest distance to the boundary. In such cases, no existing estimator has been proven to have properties such as consistency and asymptotic normality. We provide the asymptotic properties of our simulation-based estimator for a class of multidimensional RDDs, taking into account the simulation errors.
Our estimator is applicable to a class of data-generating algorithms that includes stochastic and deterministic algorithms used in practice. Our results thus nest existing insights on quasi-experimental variation in particular algorithms, such as supervised learning cowgill2018impact, bundorf2019humans, bandit, reinforcement learning, and market-design algorithms abdulkadirouglu2017research, abdulkadiroglu2019breaking, atila2013, kawai2020using, narita2021theory, narita2018toward. Our framework also reveals new sources of identification for algorithms that, at first sight, do not appear to produce a natural experiment.
The Approximate Propensity Score in this paper is related to the local random assignment interpretation of the RDD, discussed by cattaneo2015randomization, frandsen2017party, sekhon2017interpreting, frolich/huber2018, and abdulkadiroglu2019breaking. These papers consider special cases of this paper's framework.
Our empirical application uses the proposed method to study hospitals receiving CARES Act relief funding. Our empirical finding contributes to emerging work on how health care providers respond to financial shocks duggan2000hospital, dranove2017nonprofits, adelino2022hospital. Our empirical setting is a healthcare crisis, so our work complements prior work on more normal situations. Our analysis also exploits rule-based locally random assignment of cash flows to hospitals. This feature provides our estimates with additional confidence in their causal interpretation.
Our framework is a mix of the conditional independence, multidimensional RDD, and instrumental variable scenarios. In the setup in the introduction, we are interested in the effect of some binary treatment $D_i\in \{0,1\}$ on some outcome of interest $Y_i\in \mathbb{R}$. As is standard in the literature, we impose the exclusion restriction that the treatment recommendation $Z_i\in \{0,1\}$ does not affect the observed outcome other than through the treatment assignment $D_i$. This allows us to define the potential outcomes indexed against the treatment assignment $D_i$ alone. $Y_i(1)$ and $Y_i(0)$ denote potential outcomes when the individual is treated and not treated, respectively.
We consider algorithms that make treatment recommendations based solely on individual $i$'s predetermined, observable covariates $X_i=(X_{i1},...,X_{ip})'\in \mathbb{R}^p$. Let the function $A:\mathbb{R}^p\rightarrow [0,1]$ represent the decision algorithm, where $A(X_i)=\Pr(Z_i=1|X_i)$ is the probability that the treatment is recommended for individual $i$ with covariates $X_i$. The central assumption is that the analyst knows function $A$ and is able to simulate it. That is, the analyst is able to compute the recommendation probability $A(x)$ given any input value $x\in\mathbb{R}^p$. The treatment recommendation $Z_i$ for individual $i$ is then randomly determined with probability $A(X_i)$ independently of everything else. Consequently, the following conditional independence holds.
The codomain of $A$ contains $0$ and $1$, allowing for deterministic treatment assignments conditional on $X_i$. Our framework therefore nests the RDD as a special case. Another special case is the classic conditional independence scenario with the common support condition ($A(X_i)\in (0,1)$ almost surely). In addition to these simple settings, this framework nests many other situations, such as multidimensional RDDs and complex machine learning and market-design algorithms, as illustrated in Sections (ref)-(ref) and Appendix (ref).
In typical machine-learning scenarios, an algorithm first applies machine learning on $X_i$ to make some prediction and then uses the prediction to output the recommendation probability $A(X_i)$, as in the following example.\\
Example. Automated disease detection algorithms use machine learning, in particular deep learning, to detect various diseases and to identify patients at risk Gulshan2016detection. A detection algorithm predicts whether an individual $i$ has a certain disease ($Z_i=1$) or not ($Z_i=0$) based on a digital image $X_i\in\mathbb{R}^p$ of a part of the individual's body, where each $X_{ij}\in \mathbb{R}$ denotes the intensity value of a pixel in the image. The algorithm uses training data to construct a binary classifier $A:\mathbb{R}^p\rightarrow \{0,1\}$. The classifier takes an image of individual $i$ as input and makes a binary prediction of whether the individual has the disease: $Z_i\equiv A(X_i).$
The algorithm's diagnosis $Z_i$ may influence the doctor's treatment decision for the individual, denoted by $D_i\in\{0,1\}$. We are interested in how the treatment decision $D_i$ affects the individual's health outcome $Y_i$.\\
Let $Y_{zi}$ be defined as $Y_{zi}\equiv D_i(z)Y_i(1)+(1-D_i(z))Y_i(0)$ for $z\in\{0,1\}$. $Y_{zi}$ is the potential outcome when the treatment recommendation is $Z_i=z$. It follows from Property (ref) that $Z_i\indep (Y_{1i},Y_{0i})|X_i$.
We put assumptions on the covariates $X_i$ and the algorithm $A$. To simplify the exposition, the main text assumes that the distribution of $X_i$ is absolutely continuous with respect to the Lebesgue measure. Appendix (ref) extends the analysis to the case where some covariates in $X_i$ are discrete. Let ${\cal X}$ be the support of $X_i$, ${\cal X}_0=\{x\in{\cal X}:A(x)=0\}$, ${\cal X}_1=\{x\in{\cal X}:A(x)=1\}$, ${\cal L}^p$ be the Lebesgue measure on $\mathbb{R}^p$, and ${\rm int}(S)$ be the interior of a set $S\subset\mathbb{R}^p$.
Assumption (ref) (ref) allows the function $A$ to be discontinuous on a set of points with the Lebesgue measure zero. For example, $A$ is allowed to be a discontinuous step function as long as it is continuous almost everywhere. Assumption (ref) (ref) holds if the Lebesgue measures of the boundaries of ${\cal X}_0$ and ${\cal X}_1$ are zero. Assumption (ref) (ref) is only for ruling out perverse cases such as the case where $A(x)=1$ if $x\in\mathbb{R}$ is an irrational number and $A(x)\neq 1$ otherwise.
What causal effects can be learned from data $(Y_i, X_i, D_i, Z_i)$ generated by the algorithm $A$? A key step toward answering this question is what we call the Approximate Propensity Score (APS). To define it, we first define the fixed-bandwidth Approximate Propensity Score as follows:
where $B(x,\delta)=\{x^*\in\mathbb{R}^p:\|x-x^*\|<\delta\}$ is the (open) $\delta$-ball around $x\in {\cal X}$.\footnote{ We use a ball for simplicity. When we instead use a rectangle, ellipsoid, or any standard kernel function to define $p^{A}(x;\delta)$, the limit $\lim_{\delta\rightarrow 0}p^{A}(x;\delta)$ may be different at some points (e.g., at discontinuity points of $A$), but the same identification results hold under suitable conditions.} Here, $\|\cdot\|$ denotes the Euclidean norm on $\mathbb{R}^p$. To make a common bandwidth $\delta$ for all dimensions reasonable, we standardize $X_{ij}$ to have variance one for each $j=1,...,p$.
We assume that $A$ is a ${\cal L}^p$-measurable function so that the integrals exist. We then define APS as follows:
APS at $x$ is the average probability of a treatment recommendation in a shrinking ball around $x$. We call this the {\it Approximate} Propensity Score, since this score modifies the standard propensity score $A(X_i)$ to incorporate local variation in the score. APS exists for most covariate points and algorithms (see Appendix (ref)).
Figure (ref) illustrates APS. In the example, $X_i$ is two dimensional, and the support of $X_i$ is divided into three sets depending on the value of $A$. For the interior points of each set, APS is equal to $A$. On the border of any two sets, APS is the average of the $A$ values in the two sets. Thus, $p^{A}(x)=\frac{1}{2}(0+0.5)=0.25$ for any $x$ in the open line segment $AB$, $p^{A}(x)=\frac{1}{2}(0.5+1)=0.75$ for any $x$ in the open line segment $BC$, and $p^{A}(x)=\frac{1}{2}(0+1)=0.5$ for any $x$ in the open line segment $BD$.
Our identification analysis uses the following continuity condition.
Assumption (ref) is a multivariate extension of the local mean continuity condition frequently assumed in the RDD. In the RDD with a single running variable, the point $x$ for which $p^{A}(x)\in (0,1)$ and $A(x)\in \{0,1\}$ is the cutoff at which the treatment probability discontinuously changes. $A(x)\in\{0,1\}$ means that the treatment recommendation $Z_i$ is deterministic conditional on $X_i=x$. If APS at the point $x$ is nondegenerate ($p^{A}(x)\in (0,1)$), however, there exists a point close to $x$ that has a different value of $A$ from $x$'s, which creates variation in the treatment recommendation near $x$. For any such point $x$, Assumption (ref) requires that the points close to $x$ have similar conditional means of the outcome $Y_{zi}$ and treatment assignment $D_i(z)$.\footnote{ In the context of the RDD with a single running variable, one sufficient condition for continuity of $E[Y_{zi}|X_i]$ is a local independence condition in the spirit of hahn2001RDD: $(Y_i(1), Y_i(0), D_i(1), D_i(0))$ is independent of $X_i$ near $x$. A weaker sufficient condition, which allows such dependence, is that $E[Y_i(d)|D_i(1)=d_1, D_i(0)=d_0, X_i]$ and $\Pr(D_i(1)=d_1, D_i(0)=d_0|X_i)$ are continuous at $x$ for every $d\in\{0,1\}$ and $(d_1,d_0)\in\{0,1\}^2$ Dong2018. This assumes that the conditional means of the potential outcomes for each of the four types determined based on the potential treatment assignment $D_i(z)$ and the conditional probabilities of those types are continuous at the cutoff. These two sets of conditions are sufficient for continuity of $E[Y_{zi}|X_i]$ regardless of the dimension of $X_i$, accommodating multidimensional RDDs.} Note that Assumption (ref) does not require continuity of the conditional means at $x$ for which $A(x)\in(0,1)$, since the identification of the conditional means at such points follows from Property (ref) without continuity.
Under the above assumptions, APS provides an easy-to-check condition for whether an algorithm allows us to identify causal effects. We say that a causal effect is identified if it is uniquely determined by the joint distribution of $(Y_i,X_i,D_i,Z_i)$.
Proposition (ref) characterizes a necessary and sufficient condition for identification. Part (a) says that the average effects of the treatment recommendation $Z_i$ on the outcome $Y_i$ and on the treatment assignment $D_i$ conditional on $X_i=x$ are both identified if APS at $x$ is neither 0 nor 1. Nondegeneracy of APS at $x$ implies that there are both types of individuals who receive $Z_i=1$ and $Z_i=0$ among those whose $X_i$ is close to $x$. Assumption (ref) ensures that these individuals are similar in terms of average potential outcomes and treatment assignments. We can therefore identify the average effects conditional on $X_i=x$. In Figure (ref), $p^{A}(x)\in (0,1)$ holds for any $x$ in the shaded region.
Part (b) provides a necessary condition for identification. It says that if the average effect of the treatment recommendation conditional on $X_i$ being in some open set $S$ is identified, then we must have $p^{A}(x)\in (0,1)$ for almost every $x\in S$. If, to the contrary, there is a subset of $S$ of nonzero measure for which $p^{A}(x)=1$ (or $p^{A}(x)=0$), then $Z_i$ has no variation in the subset, which makes it impossible to identify the average effect for the subset.
Proposition (ref) concerns causal effects of treatment {\it recommendation}, not of treatment {\it assignment}. The proposition implies that the conditional average treatment effects and the conditional local average treatment effects (LATEs) are identified under additional assumptions.
Nondegeneracy of APS $p^{A}(x)$ therefore summarizes what causal effects the data from $A$ identify. Note that the key condition ($p^{A}(x)\in(0,1)$) holds for some points $x$ for every standard algorithm except trivial algorithms that always recommend a treatment with probability 0 or 1. Therefore, the data from every nondegenerate algorithm identify some causal effect.
The sources of quasi-random assignment characterized in Proposition (ref) suggest a way of estimating causal effects of the treatment. In view of Proposition (ref), it is possible to nonparametrically estimate conditional average causal effects $E[Y_{1i}-Y_{0i}|X_i=x]$ and $E[D_i(1)-D_i(0)|X_i=x]$ for points $x$ such that $p^{A}(x)\in(0,1)$. This approach is hard to use in practice, however, when $X_i$ has many elements. We instead seek an estimator that aggregates conditional effects at different points into a single average causal effect. Proposition (ref) suggests that conditioning on APS makes algorithm-based treatment recommendation quasi-randomly assigned. This motivates the use of an algorithm's recommendation as an instrument conditional on APS, which we operationalize as follows.
Suppose that we observe a random iid sample $\{(Y_i,X_i,D_i,Z_i)\}_{i=1}^n$ of size $n$ from the population whose data-generating process is as described in the introduction and Section (ref). Consider the following 2SLS regression using the observations with $p^{A}(X_i;\delta_n)\in (0,1)$:
where bandwidth $\delta_n$ shrinks toward zero as the sample size $n$ increases. We drop the constant term if $A(X_i)$ takes only one nondegenerate value in the sample. Let $I_{i,n}=1\{p^{A}(X_i;\delta_n)\in(0,1)\}$, $\mathbf{D}_{i,n}=(1, D_i, p^{A}(X_i;\delta_n))'$, and $\mathbf{Z}_{i,n}=(1, Z_i, p^{A}(X_i;\delta_n))'$. The 2SLS estimator $\hat\beta$ is then given by
Let $\hat\beta_1$ denote the 2SLS estimator of $\beta_1$ in the above regression.
The above regression uses true fixed-bandwidth APS $p^{A}(X_i;\delta_n)$, but it may be difficult to analytically compute if $A$ is complex. In such a case, we propose to approximate $p^{A}(X_i;\delta_n)$ using brute force simulation. We draw a value of $x$ from the uniform distribution on $B(X_i,\delta_n)$ a number of times, compute $A(x)$ for each draw, and take the average of $A(x)$ over the draws. Formally, let $X_{i,1}^*,...,X_{i,S_n}^*$ be $S_n$ independent draws from the uniform distribution on $B(X_i,\delta_n)$, and calculate $$ p^s(X_i;\delta_n)=\frac{1}{S_n}\sum_{s=1}^{S_n}A(X_{i,s}^*). $$ \sloppy We compute $p^s(X_i;\delta_n)$ for each $i=1,...,n$ independently across $i$ so that $p^s(X_1;\delta_n),...,p^s(X_n;\delta_n)$ are independent of each other. For fixed $n$ and $X_i$, the approximation error relative to true $p^{A}(X_i;\delta_n)$ has a $1/\sqrt{S_n}$ rate of convergence.\footnote{More precisely, $|p^s(X_i;\delta_n)-p^{A}(X_i;\delta_n)|=O_{p^s}(1/\sqrt{S_n})$, where $O_{p^s}$ indicates the stochastic boundedness in terms of the probability distribution of the $S_n$ simulation draws.} This rate does not depend on the dimension of $X_i$, so the simulation error can be made negligible even when $X_i$ has many elements.\footnote{When $A(X_i)\in\{0,1\}$ and $p^A(X_i;\delta_n)\in (0,1)$, the approximation error of $1\{p^s(X_i;\delta_n)\in (0,1)\}$ relative to $1\{p^A(X_i;\delta_n)\in (0,1)\}$, $|1\{p^s(X_i;\delta_n)\in (0,1)\}-1\{p^A(X_i;\delta_n)\in (0,1)\}|$, is $O_{p^s}(\sqrt{\max\{p^A(X_i;\delta_n),1-p^A(X_i;\delta_n)\}^{S_n}})$. This is because $\Var(1\{p^s(X_i;\delta_n)\in (0,1)\})=\Pr(p^s(X_i;\delta_n)\in (0,1))(1-\Pr(p^s(X_i;\delta_n)\in (0,1)))$, where $\Pr(p^s(X_i;\delta_n)\in (0,1))=1-p^A(X_i;\delta_n)^{S_n}-(1-p^A(X_i;\delta_n))^{S_n}.$ This convergence rate may become slower with the dimension of $X_i$ since $\max\{p^A(X_i;\delta_n),1-p^A(X_i;\delta_n)\}$ tends to be close to one when $X_i$ has many elements. The rate is, however, faster than $1/\sqrt{S_n}$ regardless of the dimension. Furthermore, the simulation exercise in Section (ref) shows that our method works with a moderate number of simulations ($S_n=400$) even in a high-dimensional setting ($p=100$).} We consider the simulation version of the 2SLS regression ((ref)) and ((ref)), where we use the simulated fixed-bandwidth APS $p^s(X_i;\delta_n)$ in place of $p^{A}(X_i;\delta_n)$. Let $\hat\beta_1^s$ denote the 2SLS estimator of $\beta_1$ in the simulation-based regression.
An alternative 2SLS regression is the one controlling for the variable that equals $A(X_i)$ for observations with $A(X_i)\in (0,1)$ and $p^s(X_i;\delta_n)$ for those with $A(X_i)\in \{0,1\}$. This allows us to avoid any simulation for those with $A(X_i)\in (0,1)$. Moreover, this prevents the finite-sample bias due to using fixed-bandwidth APS instead of the standard propensity score for those with $A(X_i)\in (0,1)$. This modification does not affect the asymptotic results below.
We establish the consistency and asymptotic normality of the 2SLS estimators $\hat\beta_1$ and $\hat\beta_1^s$. We use the following assumptions.
Assumption (ref) is a set of conditions for establishing consistency. Assumption (ref) (ref) assumes that the weighted average effect of the algorithm's recommendation on the treatment assignment is nonzero.\footnote{The average is taken over the covariate values for which $p^A(x)$ is nondegenerate (i.e., $p^A(x)(1-p^A(x))\in (0,1)$). There is a positive mass of such covariate values when $A$ is stochastic ($\Pr(A(X_i)\in (0,1))>0$). When $A$ is deterministic ($\Pr(A(X_i)\in (0,1))=0$), APS is nondegenerate only for boundary points at which the treatment recommendation changes from one to the other. Typically, the Lebesgue measure of the boundary is zero, so we compute the integral with respect to the $(p-1)$-dimensional Hausdorff measure instead. } Under this assumption, the estimated first-stage coefficient on $Z_i$ converges to a nonzero quantity.
Assumptions (ref) (ref)--(ref) are a set of conditions we require for proving consistency and asymptotic normality of $\hat\beta_1$ when $A$ is deterministic and produces only multidimensional regression-discontinuity variation. Assumption (ref) (ref) says that $A$ produces variation in the treatment recommendation.
Assumption (ref) (ref) imposes the differentiability of the boundary of $\Omega^*=\{x\in \mathbb{R}^p: A(x)=1\}$. The conditions are satisfied if, for example, $\Omega^*=\{x\in \mathbb{R}^p:f(x)\ge 0\}$ for some twice continuously differentiable function $f:\mathbb{R}^p\rightarrow \mathbb{R}$ such that $\nabla f(x)=(\frac{\partial f(x)}{\partial x_1},...,\frac{\partial f(x)}{\partial x_p})'\neq \bm{0}$ for all $x\in \mathbb{R}^p$ with $f(x)=0$. $\Omega^*$ takes this form when supervised learning based on smooth models such as lasso is used to construct a binary classifier $A$ such that $A(x)=1\{f(x)\ge 0\}$, $x\in\mathbb{R}^p$.
In general, the differentiability of $\Omega^*$ may not hold. For example, if tree-based algorithms such as Classification And Regression Tree (CART) and random forests are used to construct a classifier $A(x)=1\{f(x)\ge 0\}$, then the function $f$ is not differentiable at some points. Yet, the assumptions approximately hold in that $\Omega^*$ is arbitrarily well approximated by a set that satisfies the differentiability condition.
Part (ref) of Assumption (ref) (ref) says that the boundary of $\Omega^*$ is $(p-1)$ dimensional and that the boundary has nonzero density.\footnote{The boundary of $\Omega^*$ may fail to be $(p-1)$ dimensional in trivial cases where the Lebesgue measure of $\Omega^*$ is zero and hence $A(X_i)=0$ with probability one. For example, when the covariate space is three dimensional ($p=3$) and $\Omega^*$ is a straight line, not a set with nonzero volume nor even a plane, the boundary of $\Omega^*$ is the same as $\Omega^*$, and its two-dimensional Hausdorff measure is zero.} Part (ref) puts a weak restriction on the values $A$ takes on outside the support of $X_i$. It requires that $A(x)=0$ for almost every $x\notin\Omega^*$ outside ${\cal X}$ but in the neighborhood of ${\cal X}$. $A(x)$ may take on any value if $x$ is not close to ${\cal X}$. Without this assumption, it is possible that $p^A(X_i;\delta_n)\in (0,1)$ and observation $i$ is included in the regression even if $X_i$ is not near the boundary of $\Omega^*$.\footnote{For example, suppose that ${\cal X}=[-1,1]$, $A(x)=1$ if $x>0$, $A(x)=0$ if $x\in [-1,0]$, and $A(x)=\frac{1}{2}$ if $x<-1$ (note that $\Pr(A(X_i)\in (0,1))=0$ in this case). In this case, for any sufficiently small $\delta>0$, $p^A(X_i;\delta)\in (0,1)$ if $X_i\in [-1,-1+\delta)$, but these observations should not be included in the regression as they are not close to the boundary of $\Omega^*=(0,\infty)$. Such a case is ruled out by Part (ref) of Assumption (ref) (ref).} These conditions hold in practice. Assumption (ref) (ref) imposes continuity, continuous differentiability, and boundedness on the conditional moments of potential outcomes and the probability density near the boundary of $\Omega^*$. Note that Part (ref) of Assumption (ref) (ref) implies Assumption (ref).
Under the above conditions and technical regularity Assumptions (ref) and (ref) in Appendix (ref), the 2SLS estimators $\hat\beta_1$ and $\hat\beta_1^s$ are consistent and asymptotically normal estimators of a weighted average treatment effect.\footnote{Assumption (ref) is a set of additional regularity conditions (such as the smoothness of $A$) for proving asymptotic normality of $\hat\beta_1$ when $A$ is stochastic ($\Pr(A(X_i)\in (0,1))>0$). Assumption (ref) imposes the condition on the growth rate of the number of simulation draws $S_n$, which we require for proving asymptotic normality of the simulation-based estimator $\hat\beta_1^s$.}
Theorem (ref) says that the 2SLS estimators converge to the limit of a weighted average of causal effects for the subpopulation whose fixed-bandwidth APS is nondegenerate ($p^{A}(X_i;\delta)\in (0,1)$) and who would switch their treatment status in response to the treatment recommendation ($D_i(1)\neq D_i(0)$). The limit $\lim_{\delta \rightarrow 0} E[\omega_i(\delta)(Y_i(1)-Y_i(0))]$ always exists under the assumptions of Theorem (ref). It is possible to estimate other weighted averages and the unweighted average by reweighting different observations appropriately. For example, we can estimate the unweighted average treatment effect by weighting observations by the inverse of fixed-bandwidth APS. Under monotonicity ($\Pr(D_i(1)\geq D_i(0)|X_i)=1$), we could also apply abadie2003semiparametric’s Kappa weighting method using fixed-bandwidth APS instead of the standard propensity score to estimate other weighted averages of treatment effects for compliers (see also sloczynski2020should).
Theorem (ref) also shows that inference based on the conventional 2SLS heteroskedasticity-robust standard errors is asymptotically valid if $\delta_n$ goes to zero at an appropriate rate. The convergence rate of $\hat\beta_1$ is $O_p(1/\sqrt{n})$ if $\Pr(A(X_i)\in (0,1))>0$ and is $O_p(1/\sqrt{n\delta_n})$ if $\Pr(A(X_i)\in (0,1))=0$.
Our consistency result requires that $\delta_n$ go to zero slower than $n^{-1}$. The rate condition ensures that, when $\Pr(A(X_i)\in (0,1))=0$, we have sufficiently many observations in the $\delta_n$-neighborhood of the boundary of $\Omega^*$. Importantly, the rate condition does not depend on the dimension of $X_i$, unlike other bandwidth-based estimation methods such as kernel methods. This is because we use all the observations in the $\delta_n$-neighborhood of the boundary, and the number of those observations is of order $n\delta_n$ regardless of the dimension of $X_i$ if the dimension of the boundary is $p-1$, i.e., the dimension of $X_i$ minus one. When $\Pr(A(X_i)\in (0,1))>0$, this rate condition is not necessary since the effective sample size always goes to infinity at the rate $n$ regardless of the value of $\delta_n$.
The asymptotic normality requires that $\delta_n$ go to zero sufficiently quickly so that $n\delta_n^2\rightarrow 0$. When $\Pr(A(X_i)\in (0,1))>0$, we need to use a small enough $\delta_n$ so that $p^{A}(X_i;\delta_n)$ converges to $p^{A}(X_i)$ fast enough. When $\Pr(A(X_i)\in (0,1))=0$, the asymptotic normality is based on undersmoothing, which eliminates the asymptotic bias by using the observations sufficiently close to the boundary of $\Omega^*$. In both cases, the bias of our estimator is $O(\delta_n)$. The standard deviation is $O(1/\sqrt{n})$ when $\Pr(A(X_i)\in (0,1))>0$ and is $O(1/\sqrt{n\delta_n})$ when $\Pr(A(X_i)\in (0,1))=0$. The condition that $n\delta_n^2\rightarrow 0$ ensures that the bias converges to zero faster than the standard deviation in either case. When $\Pr(A(X_i)\in (0,1))=0$, a weaker condition that $n\delta_n^3\rightarrow 0$ is sufficient for eliminating the asymptotic bias.\footnote{ When $\Pr(A(X_i)\in (0,1))=0$ and $n\delta_n^3$ goes to some nonzero constant, our estimator can have a nonzero asymptotic bias (see Appendix (ref) for the asymptotic distribution).} When $\Pr(A(X_i)\in (0,1))>0$, it is possible to relax the rate condition that $n\delta_n^2\rightarrow 0$ at the cost of strengthening the smoothness of $A$ in Assumption (ref).
Whether or not $\Pr(A(X_i)\in (0,1))=0$, when we use simulated fixed-bandwidth APS, the consistency result requires that the number of simulation draws $S_n$ go to infinity as $n$ increases. The asymptotic normality result requires a sufficiently fast growth rate of $S_n$ (satisfying Assumption (ref)) to make the bias caused by using $p^s(X_i;\delta_n)$ negligible.
Finally, note that the weight $\omega_i(\delta)$ given in Theorem (ref) is negative if $D_i(1)<D_i(0)$, so $E[\omega_i(\delta)(Y_i(1)-Y_i(0))]$ may not be a causally interpretable convex combination of treatment effects $Y_i(1)-Y_i(0)$. This can happen because the treatment effect of those whose treatment assignment switches from 1 to 0 in response to the treatment recommendation (i.e., defiers) negatively contributes to $E[\omega_i(\delta)(Y_i(1)-Y_i(0))]$. Additional assumptions prevent this problem. If the treatment effect is constant, for example, the 2SLS estimators are consistent for the treatment effect.
Another approach is to impose monotonicity imbens1994. Let $LATE(x)=E[Y_i(1)-Y_i(0)|D_i(1)\neq D_i(0),X_i=x]$ be the local average treatment effect (LATE) conditional on $X_i=x$.
Theorem (ref) holds whether $A$ is stochastic ($\Pr(A(X_i)\in (0,1))>0$) or deterministic ($\Pr(A(X_i)\in (0,1))=0$). If we consider these two underlying cases separately, the probability limit of the 2SLS estimators has a more specific expression, as shown in the proof of Theorem (ref) in Appendix (ref). If $\Pr(A(X_i)\in (0,1))>0$,
The 2SLS estimators converge to a weighted average of treatment effects for the subpopulation with nondegenerate $A(X_i)$.
To relate this result to existing work, consider the following 2SLS regression with the (standard) propensity score $A(X_i)$ control:
Under conditional independence, the 2SLS estimator from this regression converges in probability to the treatment-variance weighted average of treatment effects in ((ref)) hull. Not surprisingly, for this selection-on-observables case, our result shows that the 2SLS estimator is consistent for the same treatment effect whether we control for the propensity score, fixed-bandwidth APS, or simulated fixed-bandwidth APS.
Importantly, using fixed-bandwidth APS as a control allows us to consistently estimate a causal effect even if $A$ is deterministic and produces multidimensional regression-discontinuity variation. If $\Pr(A(X_i)\in (0,1))=0$,
The 2SLS estimators converge to a weighted average of treatment effects for the subpopulation who are on the boundary of the treated region.\footnote{We prove this result using techniques from differential geometry and geometric measure theory. Our approach using geometric theory shows that $\hat\beta_1$ converges to an integral of the conditional treatment effect over boundary points with respect to the Hausdorff measure. In contrast, prior studies on multidimensional RDDs express treatment effect estimands in terms of expectations conditional on $X_i$ being in the boundary like $E[Y_{1i}-Y_{0i}|X_i\in \partial\Omega^*]$ zajonc2012regression. However, those conditional expectations are, formally, not well-defined, since ${\cal L}^p(\partial\Omega^*)=0$ and hence $\Pr(X_i\in \partial\Omega^*)=0$. We therefore prefer our expression in terms of an integral with respect to the Hausdorff measure.}
The estimand in ((ref)) nests parameters considered in the Regression Discontinuity (RD) literature. Under the monotonicity condition in Corollary (ref), the estimand in ((ref)) equals
This estimand can be interpreted as the average treatment effect for the subpopulation of compliers who are on the boundary zajonc2012regression. In the univariate RDD with a single cutoff, this estimand further reduces to the average treatment effect for the compliers at the cutoff, the standard parameter in the fuzzy RDD hahn2001RDD. If we instead assume $\Pr(D_i(1)>D_i(0)|X_i=x)=1$ (perfect compliance) for all $x\in \partial\Omega^*$, the estimand in ((ref)) equals
This estimand represents the average treatment effect for the subpopulation who are on the boundary zajonc2012regression,keele2015geographic.
To compare our method with existing ones, first consider the univariate RDD with a single cutoff $c$. In this special case, $p^{A}(X_i;\delta_n)=\frac{X_i-c}{2\delta_n}+\frac{1}{2}$ if $X_i\in[c-\delta_n,c+\delta_n]$ and $p^{A}(X_i;\delta_n)\in\{0,1\}$ otherwise. Therefore, the estimator $\hat\beta_1$ from the 2SLS regression ((ref)) and ((ref)) is numerically equivalent to a version of the regression discontinuity (RD) local linear estimator hahn2001RDD that uses a box kernel and places the same slope coefficient of $X_i$ on both sides of the cutoff.\footnote{It is possible to allow for slope changes at the cutoff by viewing $p^{A}(X_i;\delta_n)$ as a running variable with cutoff $\frac{1}{2}$ and applying standard RD local linear estimators (i.e., adding interaction terms $Z_i(p^{A}(X_i;\delta_n)-\frac{1}{2})$ and $D_i(p^{A}(X_i;\delta_n)-\frac{1}{2})$ to ((ref)) and ((ref)), respectively). Under twice continuous differentiability of $E[Y_{1i}|X_i]$ and $E[Y_{0i}|X_i]$ near the cutoff (which is stronger than our smoothness assumption), this estimator achieves the optimal rate of $n^{-2/5}$. For the multidimensional RDD with a nonlinear boundary, $Z_i$ is not a deterministic function of $p^A(X_i;\delta_n)$. In this case, it is not straightforward to use $p^A(X_i;\delta_n)$ as a single running variable, since no appropriate cutoff value exists. We leave to future research how to allow for more flexible 2SLS specifications in the general multidimensional setting.} When we use the bandwidth $\delta_n$ that converges to zero at the $n^{-1/3}$ rate instead of undersmoothing, our estimator achieves a convergence rate of $n^{-1/3}$. This rate is optimal for the estimation of the conditional LATE at the cutoff under our smoothness condition (continuous differentiability of $E[Y_{1i}|X_i]$ and $E[Y_{0i}|X_i]$ near the cutoff) in Assumption (ref) (ref).\footnote{Note that the continuous differentiability of $E[Y_{1i}|X_i]$ and $E[Y_{0i}|X_i]$ near the cutoff implies that they have bounded derivatives near the cutoff, so the latter condition can be omitted from Assumption (ref) (ref) for the univariate case.}
Our approach is particularly useful in more general scenarios. One such scenario is the multidimensional sharp or fuzzy RDD when the boundary is complex or its analytical form is unknown. Our estimator is computationally feasible for any decision boundary as long as we can simulate the underlying algorithm. More importantly, our estimator is shown to be consistent and asymptotically normal. As far as we know, there appear to be no existing estimators that are computationally feasible and have theoretical validity for a general class of multidimensional RDDs. Moreover, our method is applicable to a more general setting that mixes stochastic and deterministic algorithms. We illustrate such a case in the next section.
This section assesses the feasibility and performance of our method. We do so through a Monte Carlo experiment motivated by decision making by machine learning with high-dimensional data. Consider a government or tech company that applies a machine-learning-based deterministic decision algorithm to a large segment of the population. At the same time, they conduct a randomized controlled trial (RCT) using the rest of the population. They are interested in estimating treatment effects using data from both segments. Our approach offers a way of exploiting not only the RCT segment but also the deterministic algorithm segment. Even if we focus on the deterministic algorithm segment, it is not straightforward to apply existing RDD methods to our simulation setup since the decision boundary is high dimensional and complex. We demonstrate the applicability of our approach in such a setup.
We simulate $1,000$ hypothetical samples from the following data-generating process. Each sample $\{(Y_i,X_i,D_i,Z_i)\}_{i=1}^n$ is of size $n=10,000$. There are $100$ covariates ($p=100$), and $X_i\sim {\cal N}(\bm{0},\mathbf{\Sigma})$. $Y_i(0)$ is generated as $Y_i(0)=0.75X_i'\alpha_0+0.25 \epsilon_{0i}$, where $\alpha_0\in\mathbb{R}^{100}$, and $\epsilon_{0i}\sim {\cal N}(0,1)$. We consider two models for $Y_i(1)$, one in which the treatment effect $Y_i(1)-Y_i(0)$ does not depend on $X_i$ and one in which the treatment effect depends on $X_i$.
The choice of parameters $\mathbf{\Sigma}$, $\alpha_0$ and $\alpha_1$ is explained in Appendix (ref). $D_i(0)$ and $D_i(1)$ are generated as $D_i(0)=0$ and $D_i(1)=1\{Y_i(1)-Y_i(0)>u_i\}$, where $u_i\sim {\cal N}(0,1)$.
To generate $Z_i$, let $q_{0.495}$ and $q_{0.505}$ be the $49.5$th and $50.5$th (empirical) quantiles of the first covariate $X_{i1}$. Let $\tau_{pred}(X_i)$ be a real-valued function of $X_i$, which is constructed by random forests using an independent sample (see Appendix (ref) for the details). $Z_i$ is then generated as
The first case corresponds to the RCT segment while the latter two cases to the deterministic algorithm segment. The algorithm function $A$ is given by
Finally, $D_i$ and $Y_i$ are generated as $D_i=Z_iD_i(1)+(1-Z_i)D_i(0)$ and $Y_i=D_iY_i(1)+(1-D_i)Y_i(0)$, respectively.
{\bf Estimands and Estimators.} We consider four parameters as target estimands: ${\rm ATE}\equiv E[Y_i(1)-Y_i(0)]$; ${\rm ATE(RCT)}\equiv E[Y_i(1)-Y_i(0)|X_{i1}\in [q_{0.495}, q_{0.505}]]$; ${\rm LATE}\equiv E[Y_i(1)-Y_i(0)|D_i(1)\neq D_i(0)]$; and ${\rm LATE(RCT)}\equiv E[Y_i(1)-Y_i(0)|D_i(1)\neq D_i(0),X_{i1}\in [q_{0.495}, q_{0.505}]]$. In the case where the treatment effect does not depend on $X_i$ (Model A), ATE and LATE are the same as ATE(RCT) and LATE(RCT), respectively. In the case where the treatment effect depends on $X_i$ (Model B), the conditional effects are heterogeneous. However, since the RCT segment consists of those in the middle of the distribution of $X_{i1}$, the average effect for the RCT segment is close to the unconditional average effect. As a result, ATE is similar to ATE(RCT), and LATE is similar to LATE(RCT).
We use the data $\{(Y_i,X_i,D_i,Z_i)\}_{i=1}^n$ to estimate the treatment effect parameters. Our main approach is 2SLS with fixed-bandwidth APS controls in Theorem (ref). To compute fixed-bandwidth APS, we use $S=400$ simulation draws for each observation.
We compare our approach with two naive alternatives. The first alternative is OLS of $Y_i$ on a constant and $D_i$ (i.e., the difference in the sample mean of $Y_i$ between the treated group and untreated group) using all observations. The second alternative is 2SLS with $A(X_i)$ controls. This method uses the observations with $A(X_i)\in (0,1)$ to run the 2SLS regression of $Y_i$ on $D_i$ and $A(X_i)$ using $Z_i$ as an instrument for $D_i$ (see ((ref)) and ((ref)) in Section (ref)) and reports the coefficient on $D_i$.
For both models, the 2SLS estimator converges in probability to ${\rm LATE(RCT)}$ (equivalently, the right-hand side of equation ((ref))) whether we control for fixed-bandwidth APS or $A(X_i)$. However, 2SLS with $A(X_i)$ controls uses only the RCT segment while 2SLS with fixed-bandwidth APS controls additionally uses the individuals near the decision boundary of the deterministic algorithm (i.e., the boundary of the region for which $\tau_{pred}(x)\ge 0$). Therefore, 2SLS with fixed-bandwidth APS controls is expected to produce a more precise estimate than 2SLS with $A(X_i)$ controls if the conditional effects for those near the boundary are not far from the target estimand.
We do not apply any multidimensional RD estimators as alternatives, since there appear to be no existing RD estimators applicable to this setup. It is hard to apply distance-based RD methods since it is difficult to compute the distance from each $X_i$ to the high-dimensional random-forests decision boundary. An alternative is to use the individual's predicted effect $\tau_{pred}(X_i)$ as a univariate running variable. However, $\tau_{pred}(X_i)$ may not be a continuous variable since $\tau_{pred}$ is constructed by tree-based methods.
{\bf Performance Measures.} 2SLS with a single instrument has no moments, so we cannot consider the bias, standard deviation, or mean squared error. As an alternative, we calculate the median bias, ${\rm med}(\hat\theta)-\theta$, median absolute deviation from the median, ${\rm med}(|\hat\theta-{\rm med}(\hat\theta)|)$, and median absolute error, ${\rm med}(|\hat\theta-\theta|)$, where $\theta$ and $\hat\theta$ denote the estimand and estimator, respectively.
{\bf Results.} Table (ref) reports the median bias, median absolute deviation from the median (a measure of dispersion), and median absolute error (an overall performance measure). Panels A and B present the results for the cases where the conditional effects are homogeneous and heterogeneous, respectively. OLS with no controls is significantly biased, showing the importance of correcting for omitted variable bias. 2SLS with fixed-bandwidth APS controls achieves this goal, as demonstrated by its smaller biases across models, target parameters, and small values of the bandwidth $\delta$.
2SLS with fixed-bandwidth APS controls shows a consistent pattern; as the bandwidth $\delta$ grows, the bias increases while the absolute deviation from the median declines. For several values of $\delta$, 2SLS with fixed-bandwidth APS controls outperforms 2SLS with $A(X_i)$ controls in terms of the median absolute error. This finding implies that exploiting individuals near the multidimensional decision boundary of the deterministic algorithm can lead to better performance than using only the RCT segment.
We also evaluate our inference procedure based on Theorem (ref). Table (ref) reports the coverage probabilities of the 95% confidence intervals for LATE(RCT) constructed from the estimates and their heteroskedasticity-robust standard errors. The confidence intervals for 2SLS offer nearly correct coverage when $\delta$ is small, which supports the implication of Theorem (ref) that the inference procedure is valid when we use a sufficiently small $\delta$. Overall, Table (ref) shows that our estimator works well in this high-dimensional setting and performs better than alternative estimators.
The COVID-19 pandemic has afflicted millions of people across the country and has imposed historic challenges for hospitals and the health system. The pandemic led to revenue losses coupled with skyrocketing expenses, pushing many already overburdened hospitals further to their financial brink.
To deal with this crisis, as part of the 3-phase Coronavirus Aid, Relief, and Economic Security (CARES) Act, the US government distributed tens of billions of dollars of relief funding to hospitals since April 2020. This funding intended to help health care providers hit hardest by the COVID-19 outbreak. The bill specified that providers may (but are not required to) use the funds for COVID-19-related expenses such as construction of temporary structures, purchasing medical supplies and equipment (including personal protective equipment and testing supplies), increased workforce utilization and training, establishing emergency operation centers, retrofitting facilities, and managing the surge in capacity.
We ask whether this funding had a causal impact on hospital operation and activities in dealing with COVID-19 patients. Answering this question would help the government design better funding policies to respond to future healthcare crises. We focus on an initial portion of this funding (\$10 billion). This portion was allocated to hospitals that qualified as “safety net hospitals” according to a specific eligibility criterion. This eligibility criterion intends to direct funding towards hospitals that “disproportionately provide care to the most vulnerable, and operate on thin margins.” Specifically, an acute care hospital is deemed eligible for funding if the following conditions hold:
Hospitals that do not qualify on any of the three dimensions are funding ineligible. Figure (ref) visualizes how the three dimensions determine funding eligibility. From the original space of the three eligibility determinants, we extract two-dimensional planes to better visualize the structure of quasi-experimental variation. As the bottom two-dimensional planes show, eligibility discontinuously changes as hospitals cross the eligibility boundary in the characteristic space. This setting is a three-dimensional RDD, falling under our framework.
Our treatment is the funding amount, which is calculated as follows. Each eligible hospital is assigned a facility score, which is calculated as the product of DPP and the number of beds in that hospital. This facility score determines the share of funding allocated to the hospital, out of the total \$10 billion. The share received by each hospital is determined by the ratio of the hospital's facility score to the sum of facility scores across all eligible hospitals. The amount of funding that can be received by a hospital is bounded below at \$5 million and capped above at \$50 million. We compute the funding eligibility status as well as the amount of funding received, by using data from the Healthcare Cost Report Information System (HCRIS) for the 2018 financial year.\footnote{We use the methodology detailed in the \href{https://www.hhs.gov/coronavirus/cares-act-provider-relief-fund/general-information/index.html}{CARES Act website} to project funding based on 2018 financial year cost reports. We use the RAND cleaned version of the dataset rand2018. See Appendix (ref) for details on the construction of our dataset.}
A majority of eligible hospitals receive the minimum amount of \$5 million. A small mass of hospitals receive amounts close to the maximum of \$50 million. Figure (ref) in Appendix (ref) shows the distribution of funding amounts received by eligible hospitals.
Our outcomes are a few different versions of the number of COVID patients hospitalized at each hospital. To obtain these outcomes, we use the publicly available COVID-19 Reported Patient Impact and Hospital Capacity by Facility dataset hhs2021. This provides facility-level data on hospital utilization aggregated on a weekly basis, from July 31st 2020 onwards. Summary statistics about hospital outcomes and characteristics are documented in Table (ref). Eligible hospitals have larger numbers of inpatient and ICU beds occupied by COVID-19 patients. Eligible hospitals also have a higher disproportionate patient percentage, higher uncompensated care per bed, lower profit margins, more employees and beds, and shorter lengths of inpatient stay. These patterns are consistent with the funding's goal of helping struggling hospitals.
We first validate our method by evaluating the balancing property of fixed-bandwidth APS conditioning. We calculate fixed-bandwidth-APS-controlled differences in covariate means for hospitals who are and are not eligible for funding. We run the following OLS regression of hospital-level characteristics on the eligibility status using observations with $p^s(X_i;\delta)\in (0,1)$: $$ W_i = \gamma_0+\gamma_1Z_i+\gamma_2p^s(X_i;\delta) + \eta_i, $$ where $W_i$ is one of the predetermined characteristics of the hospital, $Z_i$ is a funding eligibility dummy, $X_i$ is a vector of the three input variables (DPP, UCC, and profit margin) that determine the funding eligibility, and $p^s(X_i; \delta)$ is the simulated fixed-bandwidth APS. We compute fixed-bandwidth APS using $S=10,000$ simulation draws for different bandwidth values.\footnote{We standardize the three input variables when computing fixed-bandwidth APS. Figure (ref) in Appendix (ref) reports fixed-bandwidth APS for several hospitals with varying numbers of simulation draws. We find that $S=10,000$ is sufficient for well stabilizing fixed-bandwidth APS simulation.} The estimated coefficient on $Z_i$ is the fixed-bandwidth-APS-controlled difference in the mean of the covariate between eligible and ineligible hospitals. For comparison, we also run the OLS regression of hospital characteristics on the eligibility status with no controls using the whole sample.
Table (ref) reports the covariate balance estimates. Column 2 shows that, without controlling for fixed-bandwidth APS, eligible hospitals are significantly different from ineligible hospitals. All the relevant hospital eligibility characteristics are strongly associated with eligibility. Once we control for fixed-bandwidth APS with small enough bandwidth $\delta$, eligible and ineligible hospitals have similar financial and utilization characteristics, as reported in columns 3--7 of Table (ref). These estimates are consistent with our theoretical results, establishing the empirical ability of fixed-bandwidth APS controls to eliminate selection bias.
The balancing performance of fixed-bandwidth APS motivates us to estimate the causal effects of funding by using algorithmically-determined funding eligibility as an IV. We study the effect of funding on relevant hospital outcomes, such as the number of inpatient beds occupied by adult COVID patients between July 31st 2020 and August 6th 2020.
We first run the following OLS (reduced-form) regression of each outcome on the binary funding eligibility, while controlling for fixed-bandwidth APS:
where $Y_i$ is a hospital-level outcome and $Z_i$ is the binary indicator for funding eligibility. This OLS (reduced-form) specification is a special case of the 2SLS specification treated in the theoretical analysis.
We then estimate the following 2SLS regression using the funding amount as the treatment and funding eligibility as an instrument. We run the regression on two different hospital-level outcome variables, using hospitals with $p^s(X_i;\delta)\in (0,1)$:
where $D_i$ is the funding amount. We also run the OLS and 2SLS regressions with no controls, as well as OLS regression controlling for the three eligibility determinants. These alternative regressions are computed using the sample of all hospitals, as benchmark estimators.
The first-stage effects of funding eligibility on the funding amount are substantial (columns 4--10 of Table (ref)). In column 4 of Table (ref), funding eligibility increases the amount of funding received by \$15.85 million on average.
OLS estimates of funding effects, reported as the benchmark in column 1 of Table (ref), indicate that funding is associated with a higher number of adult inpatient beds and higher number of staffed ICU beds utilized by COVID patients. For example, the estimates indicate that receiving an additional \$1 million in funding is associated with 4.53 more beds occupied by patients. These uncontrolled OLS estimates show a similar picture as the descriptive statistics in Table (ref). Naive 2SLS estimates with no controls and OLS with covariate controls produce similar significantly positive associations of funding with outcomes.
However, the OLS or uncontrolled 2SLS estimates turn out to be an artifact of selection bias. In contrast with them, our preferred reduced-form and 2SLS estimates with fixed-bandwidth APS controls show a different picture (columns 4--10). The gains in the number of inpatient beds and staffed ICU beds occupied by COVID patients become much smaller and lose significance across all bandwidth specifications. In fact, even the sign of the estimated funding effect is reversed for several combinations of the outcome and bandwidth. Once we control for fixed-bandwidth APS to eliminate the bias, therefore, funding has little to no effect on the hospital utilization level by COVID-19 patients. These results suggest that fixed-bandwidth APS reveals important selection bias in the naively estimated effects of funding.\footnote{The 2SLS estimates in Table (ref) are unlikely to be compromised by differential attrition. Estimates reported in Table (ref) in Appendix (ref) show little difference in outcome availability rates between eligible and ineligible hospitals once we control for fixed-bandwidth APS.}
We also estimate the evolving effects of funding for each week from July 31st, 2020 to April 2nd, 2021 (Appendix (ref)). The estimated dynamic effects are similar to the initial null effects, suggesting that funding has no substantial effect even in the long run. Furthermore, we estimate dynamic effects for different groups of hospitals defined by hospital size and ownership type. We do not find any strong evidence of heterogeneity at any point in time.
The overall insignificance of the estimates suggests that funding by the CARES Act had little effect on hospital utilization during the pandemic. The null effect is widely observed for subgroups of hospitals at different points in time. This finding is consistent with policy and media arguments that CARES Act funding was not well targeted toward needy providers. Unlike the previous media arguments and descriptive analyses, the analysis here provides causal evidence supporting the concern.
As algorithmic decisions become the new norm, the world becomes a mountain of natural experiments. We develop a general method to use these algorithm-produced instruments to identify and estimate causal treatment effects. Our analysis of the CARES Act hospital relief funding uses the proposed method to find that relief funding has little effect on COVID-19-related hospital activities. OLS or uncontrolled 2SLS estimates, by contrast, show considerably larger and more significant effects. The large estimates appear to be an artifact of selection bias; relief funding just went to hospitals with more COVID-19 patients, without helping hospitals accommodate additional patients.
Our analysis provides a few implications for policy and management practices of decision-making algorithms. It is important to record the implementation of algorithms in a replicable way, including what input variables $X_i$ are used to make algorithmic recommendation $Z_i$. Another key lesson is the importance of recording an algorithm's recommendation $Z_i$ even if they are superseded by a human decision $D_i$. These data retention efforts go a long way to exploit the full potential of algorithms as natural experiments.
An important topic for future research is estimation and inference details, such as data-driven bandwidth selection. This work needs to extend imbens2012bandwidth and Calonico2014robustRD's bandwidth selection methods in the univariate RDD to our setting.\footnote{For univariate RDDs, imbens2012bandwidth and Calonico2014robustRD estimate the bandwidth that minimizes the asymptotic mean squared error (AMSE). It is not straightforward to estimate the AMSE-optimal bandwidth in our setting with many running variables and complex IV assignment, since it requires nonparametric estimation of functions on the multidimensional covariate space such as conditional mean functions, their derivatives, the curvature of the RDD boundary, etc.} Inference on treatment effects in our framework relies on large sample reasoning. It seems natural to additionally consider permutation or randomization inference. It will also be challenging but interesting to develop finite-sample optimal estimation and inference strategies such as those recently introduced by Armstrong2018inference, armstrong2020finitesample and imbens2019optimized. Finally, we look forward to empirical applications of our method in a variety of business, policy, and scientific domains.