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.
146,776 characters · 15 sections · 127 citation commands
Policy Learning with Observational Data
\ifecma \setstretch{1.4} \fi
The problem of learning treatment assignment policies, or mappings from individual characteristics to treatment assignments, is ubiquitous in applied economics and statistics. \begingroup \footnote{Forthcoming in Econometrica. We circulated an earlier draft of this paper under the title “Efficient Policy Learning”; the current title was chosen following guidance from the review process. We are grateful for helpful conversations with colleagues including Victor Chernozhukov, David Hirshberg, Guido Imbens, Michael Kosorok, Alexander Luedtke, Eric Mbakop, Whitney Newey, Xinkun Nie, Molly Offer-Westort, Alexander Rakhlin, James Robins, Erik Sverdrup, Max Tabord-Meehan and Zhengyuan Zhou, and for feedback from the editor, referees, as well as seminar participants at a variety of universities and workshops. We thank Guido Imbens for sharing the GAIN dataset with us. Financial support was provided by the Sloan Foundation, Office of Naval Research grant N00014-17-1-2131, National Science Foundation grant DMS-1916163, and a Facebook Faculty Award.} \addtocounter{footnote}{-1} \endgroup It arises, for example, in medicine when a doctor must decide which patients to refer for a risky surgery; in marketing when a company needs to choose which customers to send targeted offers to; and in government and policy settings, when assigning students to educational programs or inspectors to buildings and restaurants.
The treatment assignment problem rarely arises in an unconstrained environment. Treatments are often expensive, and so a policy may need to respect budget constraints. Policies may need to be implemented in environments characterized by human or machine constraints; for example, emergency medical professionals or police officers may need to implement decision policies in the field, where a simple decision tree might be used. For internet or mobile services, algorithms may need to determine the set of information displayed to a user very quickly, and a simple lookup table may decrease the time it takes to respond to a user's request. Fairness constraints may require a treatment assignment policy to depend only on particular types of covariates (for example, test scores or income), even when other covariates are observed.
This paper is about using observational data to learn policies that respect the types of constraints outlined above. The existing literature on policy learning has mostly focused on the setting where we want to optimize allocation of a binary treatment using data from a randomized trial, or from a study with a known, random treatment assignment policy. In many problems, however, one may need to leverage richer forms of observational data to learn treatment assignment rules. For example, if we want to learn whom to prescribe a drug to based on data from a clinical trial, we need to have methods that deal with non-compliance and resulting endogenous treatment assignments.\footnote{If we believed that compliance patterns when we deploy our policy would be similar to those in the clinical trial, then an intent-to-treat analysis may be a reasonable way to side-step endogeneity concerns. However, if we suspect that compliance patterns may change (e.g., if patients may be more likely to adhere to a treatment regime prescribed by their doctor than one randomly assigned in a clinical trial), then using an analysis that disambiguates received treatment from assigned treatment is necessary.} Or, if we are interested in offering some customers discounts, then we need methods that let us study interventions to continuous variables (e.g., price) rather than just discrete ones. The goal of this paper is to develop methods for policy learning that don't just work in randomized trials (or related settings), but can instead work with a rich variety of observational designs.
Formally, we study the problem where we have access to observational data and want to use it to learn a policy that maps a subject's characteristics $X_i \in \xx$ to a binary decision, $\pi : \xx \rightarrow \cb{0, \, 1}$. The practitioner has also specified a class $\Pi$ that encodes problem-specific constraints pertaining to budget, functional form, fairness, etc., and requires that our learned policy \smash{$\hpi$} satisfies these constraints, \smash{$\hpi \in \Pi$}. Then, following manski2004statistical,manski2009identification, hirano2009asymptotics, stoye2009minimax,stoye2012minimax and kitagawa2015should, we seek guarantees on the regret \smash{$R(\hpi)$}, i.e., the difference between the expected utility from deploying the learned policy \smash{$\hpi$} over a target population and the best utility that could be achieved from deploying any policy in the class $\Pi$ over the population.
Our paper builds on a rich literature at the intersection of econometrics, statistics and computer science on learning structured treatment assignment rules, including kitagawa2015should, swaminathan2015batch and \citet*{zhao2012estimating}. Most closely related to us, kitagawa2015should study a special case of our problem where treatments are binary and exogenous with known assignment probabilities, and show that an algorithm based on inverse-probability weighting achieves regret that depends optimally on the sample size and the complexity of the policy class $\Pi$.\footnote{kitagawa2015should also consider the case where treatment assignment probabilities are unknown; in this case, however, their method no longer achieves optimal dependence on the sample size.}
Here, we develop a new family of algorithms that achieve regret guarantees with optimal dependence on sample size and on $\Pi$, but under considerably more generality on the sampling design. We consider both the classical case where we want to optimize a binary treatment, and a related setting where we want to optimize infinitesimal nudges to a continuous treatment (e.g., a price). Moreover, our approach can leverage observational data where the treatment assignment mechanism may either be exogenous with unknown assignment probabilities, or endogenous, in which case we require an instrument.
Our approach starts from recent unifying results of \citet*{chernozhukov2016locally} on semiparametrically efficient estimation. As discussed in more detail in Section (ref), chernozhukov2016locally show that in many problems of interest, we can construct efficient estimates of average-treatment-effect-like parameters $\theta$ as
where \smash{$\hGamma_i$} is an appropriate doubly robust score for the target estimand under the intervention of interest. This approach can be used to target the average effect of a binary treatment, the average derivative of a continuous treatment, and other related estimands.
In this paper we find that, whenever one can estimate the average utility of treating everyone\footnote{Throughout this paper, we assume that there is no interference, i.e., assigning one unit to treatment doesn't affect outcomes for others. For a discussion of treatment effect estimation under intereference, see hudgens2008toward, manski2013identification, and references therein.} using an estimator of the type (ref) built via the doubly robust construction of chernozhukov2016locally, we can also usefully learn whom to target with the intervention via a simple procedure: Given a pre-specified policy class $\Pi$ (e.g., linear decision rules or finite-depth decision trees), we propose using the treatment assignment rule $\hpi$ that solves\footnote{If this optimization problem has multiple solutions, we set $\hpi$ to an arbitrary maximizer of the objective. Our formal results apply simultaneously to all solutions of (ref).}
where \smash{$\hGamma_i$} are the same doubly robust scores as used in (ref). Our main result is that, under regularity conditions, the resulting policies \smash{$\hpi$} have regret \smash{$R(\hpi)$} bounded on the order of \smash{$\sqrt{\VC(\Pi)/n}$} with high probability. Here, $\VC(\Pi)$ is the Vapnik-Chervonenkis dimension of the class $\Pi$ and $n$ is the sample size. We also highlight how the constants in this bound depend on fundamental quantities from the semiparametric efficiency literature.
Our proof combines results from semiparametrics with carefully tailored analysis tools that build on classical ideas from empirical process theory. The reason we obtain strong guarantees for the approach (ref) is closely tied to robustness properties of the estimator (ref). In the setting where we only want to estimate a single average effect parameter, it is well known that non-doubly robust estimators can also be semiparametrically efficient \citep*{hirano2003efficient}. Here, however, we need convergence results that are strong enough to withstand optimization over the whole class $\Pi$. The fact that doubly robust estimators are fit for this task is closely related to their ability to achieve semiparametric efficiency under general conditions, even if nuisance components are estimated via black-box machine learning methods for which we can only guarantee fast enough convergence in mean-squared error chernozhukov2016double,van2011targeted.
We spell out our general framework in Section (ref). For intuition, however, it is helpful to first consider this approach in the simpler case where we want to study the effect of a binary treatment $W_i \in \cb{0, \, 1}$ on an outcome $Y_i \in \RR$ interpreted as a utility and are willing to assume selection on observables (unconfoundedness): We have potential outcomes $\cb{Y_i(0), \, Y_i(1)}$ such that $Y_i = Y_i(W_i)$ and $\cb{Y_i(0), \, Y_i(1)} \indep W_i \cond X_i$ imbens2015causal. Then, the utilitarian regret of deploying a policy \smash{$\pi \in \Pi$} is manski2009identification
and we can construct our estimator (ref) using the well known augmented inverse-propensity weighted scores of \citet*{robins1994estimation},\footnote{See Section (ref) for a detailed discussion of how to implement our policy learner (ref) based on these augmented inverse-propensity weighted scores in practice.}
where \smash{$\he(x)$} and \smash{$\hatm(x, \, w)$} denote non-parametric estimates of $e(x)$ and $m(x, \, w)$ respectively. In this setup, our result implies that---under regularity conditions---the estimator (ref) with scores (ref) has regret (ref) bounded on the order of \smash{$\sqrt{\VC(\Pi)/n}$}.
Even in this simplest case, our result is considerably stronger than results currently available in the literature. The main result of kitagawa2015should is that, if treatment propensities $e(X_i)$ are known, then a variant of inverse-propensity weighted policy learning achieves regret on the order of \smash{$\sqrt{\VC(\Pi)/n}$}. However, in observational studies where the treatment propensities are unknown, the bounds of kitagawa2015should depend on the rate at which we can estimate $e(\cdot)$, and will generally decay slower than $1/\sqrt{n}$. The only other available $1/\sqrt{n}$-bounds for policy learning in observational studies with a binary treatment that we are aware of are a result of \citet*{van2006cross} for the case where $\Pi$ consists of a finite set of policies whose cardinality grows with $n$, and a result of kallus2017balanced in the special case $m(\cdot, \, w)$ is assumed to belong to a reproducing kernel Hilbert space. The idea of using doubly robust scores to learn optimal treatment assignment of a binary treatment has been previously discussed in \citet*{langford2011doubly} and \citet*{zhang2012estimating}; however, neither paper provides a regret bound for this approach.
In the more general case where the observed treatment assignments $W_i$ may be continuous and/or we may need to use instrumental variables to identify causal effects, both the methods and regret bounds provided here are new. By connecting the policy learning problem to the semiparametric efficiency literature, we are able to develop a general framework that applies across a variety of settings.
The literature on optimal treatment allocation has been rapidly expanding across several fields. In the econometrics literature, the program of learning regret-optimal treatment rules was started by manski2004statistical,manski2009identification. One line of work considers the case where the policy class is unrestricted, and the optimal treatment assignment rule simply depends on the sign of the conditional average treatment effect for each individual unit. In this setting, hirano2009asymptotics show that when $1/\sqrt{n}$-rate estimation of the conditional average treatment effect function is possible, then treatment assignment rules obtained by thresholding an efficient estimate of the conditional average treatment effect are asymptotically minimax-optimal. Meanwhile, stoye2009minimax derives finite sample minimax decision rules in a class of problems where both the response surfaces and the policies $\pi$ may depend arbitrarily on covariates. Further results are given in armstrong2015inference, bhattacharya2012inferring, chamberlain2011bayesian, dehejia2005program, kasy2016partial, stoye2012minimax and tetenov2012statistical.
Building on this line of work, kitagawa2015should study policy learning in a non-parametric setting where the learned policy $\hpi$ is constrained to belong to a structured class $\Pi$ and show that, in this case, we can obtain regret bounds relative to the best policy in $\Pi$ that scale with the complexity of the class $\Pi$. A key insight from kitagawa2015should is that, when propensity scores are known and $\Pi$ has finite VC dimension, it is possible to get $1/\sqrt{n}$-rate regret bounds for policy learning over a class $\Pi$ even if the conditional average treatment effect function itself cannot be estimated at a $1/\sqrt{n}$-rate; in other words, we can reliably find a nearly best-in-class policy without needing to accurately estimate a model that describes all causal effects. As discussed above, our paper builds on this work by considering rate-optimal regret bounds for best-in-class policy learning in observational studies where propensity scores are unknown and treatment assignment may be endogenous, etc.
One difference between our results and those of kitagawa2015should is that the latter provide finite sample regret bounds, whereas our results are asymptotic in the sample size $n$. The reason for this is that our bounds rely on results from the literature on semiparametric estimation \citep*{bickel,chernozhukov2016locally,chen2008semiparametric,hahn1998role,newey1994asymptotic,robins1}, which themselves are asymptotic. Recently, armstrong2017finite showed that, in a class of average treatment effect estimation problems, finite sample conditionally minimax linear estimators are asymptotically efficient, thus providing a connection between desirable finite sample guarantees and asymptotic optimality. It would be interesting to examine whether similar connections are possible in the policy learning case.
Policy learning from observational data has also been considered in parallel literatures developed in both statistics \citep*{luedtke2016statistical,qian2011performance,zhang2012estimating,zhao2012estimating} and machine learning \citep*{beygelzimer2009offset,langford2011doubly,kallus2017balanced,swaminathan2015batch}. Two driving themes behind these literatures are the development of performant algorithms for solving the empirical maximization problems (and relaxations thereof) that underlie policy learning, and the use of doubly robust objectives for improved practical performance. kallus2017balanced, swaminathan2015batch and zhao2012estimating also prove regret bounds for their methods; however, they do not achieve a $1/\sqrt{n}$ sample dependence, with the exception of kallus2017balanced in the special case of the reproducing kernel Hilbert space setting described above. Finally, luedtke2017faster propose a class of regret bounds that decay faster than $1/\sqrt{n}$ by exploiting non-uniform asymptotics; see Section (ref) for a further discussion.
The problem of optimal treatment allocation can also be seen as a special case of the broader problem of optimal data-driven decision making. From this perspective, our result is related to the work of rudin2018big and bertsimas2014predictive, who study data-driven rules for optimal inventory management and related problems. Much like in our case, they advocate learning with a loss function that is directly tied to a utility-based criterion. Finally, we note a growing literature on estimating conditional average treatment effects, including athey2015machine, \citet*{athey2018generalized}, nie2017learning, and references therein. Although the goal is similar to that of learning optimal treatment assignment rules, the specific results themselves differ; they focus on squared-error loss rather than utilitarian regret.
Our goal is to learn a policy $\pi \in \Pi$ that maps a subject's features $X_i \in \xx$ to a treatment decision: $\pi : \xx \rightarrow \cb{0, \, 1}$. In order to do so, we assume that we have independent and identically distributed samples $(X_i, \, Y_i, \, W_i, \, Z_i)$, where $Y_i \in \RR$ is the outcome we want to intervene on, $W_i$ is the observed treatment assignment, and $Z_i$ is an (optional) instrument used for identifying causal effects. In cases where $W_i$ is exogenous, we simply take $Z_i = W_i$. Throughout our analysis, we interpret $Y_i$ as the utility resulting from our intervention on the $i$-th sample, e.g., $Y_i$ could measure the benefit accrued by a subject minus a potentially personalized cost of treatment (in Section (ref) we demonstrate inclusion of linear costs in the context of an application). We then seek policies that make the expected value of $Y_i$ large.
We define the causal effect of the intervention $\pi(\cdot)$ in terms of the potential outcomes model neyman1923applications,rubin1974estimating, whereby the $\cb{Y_i(w)}$ correspond to utilities we would have observed for the $i$-th sample had the treatment been set to $W_i = w$, and $Y_i = Y_i(W_i)$. When instruments are present, we always assume that the exclusion restriction holds so that this notation is well specified. We consider both examples with a binary treatment $W_i \in \cb{0, \,1}$ and with a continuous treatment $W_i \in \RR$.
In the case where $W_i$ is binary, we follow the existing literature hirano2009asymptotics,kitagawa2015should,manski2004statistical,stoye2009minimax, and study interventions that directly specify the treatment level. In this case, the utility of deploying a policy $\pi(\cdot)$ relative to treating no one is manski2009identification
and the corresponding policy regret relative to the best possible policy in the class $\Pi$ is
As discussed in the introduction, in this binary setting, kitagawa2015should show that if $W_i$ is exogenous with known treatment propensities, then we can use inverse-propensity weighting to derive a policy \smash{$\hpi$} whose regret \smash{$R(\hpi)$} decays as $1/\sqrt{n}$, with
Here, we develop methods that can also be used in observational studies where treatment propensities may be unknown, and where we may need to use instrumental variables to identify $V(\pi)$ from (ref).
Meanwhile, when $W_i$ is continuous, we study infinitesimal interventions on the treatment level motivated by the work of \citet*{powell1989semiparametric}. We define the utility of such an infinitesimal intervention as
and then define regret in terms of $V(\pi)$ as in (ref). One interesting conceptual difference that arises in this case is that, now, our interventions $\pi(X_i) \in \cb{0, \, 1}$ and observed treatment assignments $W_i \in \RR$ may take values in different spaces. This can arise, for example, if we want to target customers with personalized discounts and have access to past prices $W_i$ that take on a continuum of values, but are restricted to considering a class of interventions that only allow us to make a binary decision $\pi(X_i) \in \cb{0, \, 1}$ on whether to offer each customer a small discount or not. The fact that we can still learn low-regret policies via the simple strategy (ref) even when these two spaces are decoupled highlights the richness of the policy learning problem.\footnote{Another interesting question one could ask is how best to optimize the assignment of $W_i$ globally rather than locally (i.e., the case where we can set the treatment level $w$ to an arbitrary level, rather than simply nudge the pre-existing levels of $W_i$). This question would require different formal tools, however, as the results developed in this paper only apply to binary decisions.}
With both binary and continuous treatments, the regret of a policy $\pi$ can be written in terms of a conditional average treatment effect function,
such that $V(\pi) = \EE{\pi(X_i) \tau(X_i)}$ and regret $R(\pi)$ is as in (ref). Our analysis pertains to any setup with a regret function $R(\pi)$ that admits such a representation. Given these preliminaries, recall that our goal is to learn low regret policies, i.e., to use observational data to derive a policy $\hpi \in \Pi$ with a guarantee that $R(\hpi) = \oo_P\p{1/\sqrt{n}}$. In order to do so, we need to make assumptions on the observational data generation distribution that allow for identification and adequate estimation of $V(\pi)$, and also control the size of $\Pi$ in a way that makes emulating the best-in-class policy a realistic objective. The following two subsections outline these required conditions; our main result is then stated in Section (ref).
In order to learn a good policy $\hpi$, we first need to be able to evaluate $V(\pi)$ for any specific policy $\pi$. Our main assumption, following \citet*{chernozhukov2016locally}, is that we can construct a doubly robust score for the average treatment effect $\theta = \EE{\tau(X_i)}$. At the end of this section we discuss how this approach applies to three important examples, and refer the reader to chernozhukov2016locally for a more general discussion of when such doubly robust scores exist.
Given this setup, chernozhukov2016locally propose first estimating $g(\cdot)$ and $m(\cdot)$, and then consider
They show that this estimator is $\sqrt{n}$-consistent and asymptotically unbiased Gaussian for $\theta$, provided that the nuisance estimates $\hg(\cdot)$ and $\hatm(\cdot)$ converge sufficiently fast and that we use cross-fitting chernozhukov2016double,schick1986asymptotically. This estimator is also semiparametrically efficient under general conditions newey1994asymptotic.\footnote{Our results don't depend on efficiency of (ref); rather, we only use $\sqrt{n}$-consistency. In cases where (ref) may not be efficient, our regret bounds still hold verbatim; the only difference being that we can no longer interpret the terms of the form \smash{$\EE{\Gamma_i^2}$} appearing in the bound as related to the semiparametric efficient variance for $\theta$.}
Our approach to policy learning builds on these foundations. We again start by estimating nuisance components and by forming doubly robust scores as in (ref). However, instead of just averaging the \smash{$\hGamma_i$} to estimate $\theta$, we use these scores for policy learning by plugging them into (ref). Our main result will establish that we can get strong regret bounds for learning policies under conditions that are similar to those used by chernozhukov2016locally to show asymptotic normality of (ref) and, more broadly, that build on assumptions often made in the literature on semiparametric efficiency \citep*{bickel,chen2008semiparametric,hahn1998role,newey1994asymptotic,robins1}.
As in the recent work of chernozhukov2016double on double machine learning or that of van2011targeted on targeted learning, we take an agnostic view on how the nuisance estimates $\hg(\cdot)$ and $\hatm(\cdot)$ are obtained, and simply impose high level conditions on their rates of convergence. Given sufficient regularity, we can construct estimators that satisfy the rate condition (ref) via, e.g., sieve-based methods \citep*{chen2007large} or kernel regression \citep*{caponnetto2007optimal}. Moreover, in applications, we may want to consider several different machine learning methods for each component, or potentially combinations thereof, and then use cross-validation to choose which method to use. For completeness, we allow problem specific quantities to change with the sample size $n$, and track this dependence with a subscript $n$, e.g., $m_n(x, \, w) = \EE[n]{Y_i(w) \cond X_i = x}$, etc.
We end this section by verifying that Assumption (ref) in fact covers several settings of interest, and is closely related to several standard approaches to semiparametric inference. In cases with selection on observables we do not need an instrument (or can simply set $Z_i = W_i$), so for simplicity of notation we replace all instances of $Z_i$ with $W_i$.
\paragraph{Binary treatment with selection on observables.} Most existing work on policy learning, including kitagawa2015should, has focused on the setup where $W_i$ is binary and unconfounded, i.e., $\cb{Y_i(0), \, Y_i(1)} \indep W_i \cond X_i$. In this case, weighting by the inverse propensity score lets us recover the average treatment effect, i.e., $g(x, \, w) = (w - e(x)) / (e(x)(1 - e(x)))$ with $e(x) = \PP{W_i = 1 \cond X_i = x}$ identifies the conditional average treatment effect \smash{$\tau_m(x) = m(x, \, 1) - m(x, \, 0)$} via (ref). The estimation strategy (ref) yields
and recovers augmented inverse propensity weighting \citep*{robins1994estimation}.
\paragraph{Continuous treatment with selection on observables.} In the case where $W_i$ is continuous and unconfounded $\cb{Y_i(w)} \indep W_i \cond X_i$, we can derive a representer $g(\cdot)$ via integration by parts \citep*{powell1989semiparametric}. Under regularity conditions, the $\tau$-function \smash{$\tau_m(x, \, w) = [{d}/{d\nu} \ m(x, \, w + \nu)]_{\nu = 0}$} can be identified via (ref) using
where $f(\cdot \cond x)$ denotes the conditional density of $W_i$ given $X_i = x$. The resulting doubly robust estimator was to our knowledge first derived via the general approach of chernozhukov2016locally, which in turn is closely related to an approach proposed by ai2007estimation.
\paragraph{Binary, endogenous treatment with binary treatment and instrument.} Instead of unconfoundedness, now suppose that $Z_i$ is a valid instrument conditionally on features $X_i$ in the sense of Assumption 2.1 of abadie2003semiparametric. Suppose moreover that treatment effects are homogenous, meaning that the conditional average treatment effect matches the conditional local average treatment effect imbens1994late,\footnote{As discussed above, our notation has potential outcomes $Y_i(W_i)$ that only depend on treatment $W_i$, and do not involve the instrument $Z_i$. This is only meaningful when the exclusion restriction holds.}
Then we can use a weighting function $g(\cdot)$ defined in terms of the compliance score abadie2003semiparametric,aronow2013beyond,
to identify this $\tau$-function using (ref). We note that our formal results all require that $g(\cdot)$ be bounded, which implicitly rules out the case of weak instruments (since if $\Delta$ approaches 0, the $g(\cdot)$-weights blow up).
Next, in order to obtain regret bounds that decay as $1/\sqrt{n}$, we need some control over the complexity of the class $\Pi$ (and again let $\Pi$ potentially change with $n$ for generality). The Vapnik-Chervonenkis (VC) approach vapnik2000nature presents us with a natural way to do so. Recall that the VC-dimension of a class $\Pi$ of binary decision rules is the largest value of $d \in \NN$ such that there exists a set of $d$ points $x_1, \, ..., \, x_d \in \xx$ that is “shattered” by $\Pi$ in the following sense: For each $2^d$ of the binary vectors $v \in \cb{0, \, 1}^d$, there exists a policy $\pi_v \in \Pi$ such that $\pi_v(X_i) = v_i$ for all $i = 1, \, ..., \, d$. Throughout our analysis, we control the complexity of $\Pi_n$ by assuming that its VC-dimension does not grow too fast with the sample size $n$. As is familiar from the literature on classification, we will find that the best possible uniform regret bounds scale as $\sqrt{\VC(\Pi_n)/n}$ vapnik2000nature.
In order to illustrate this assumption, we give two examples of policy classes that have a finite VC dimension, and one that does not. In all three examples below, we assume that the features $X_i$ take values in $\xx = \RR^p$ for some $p \geq 1$.
\paragraph{Linear Rules} The VC-dimension of the class of linear decision rules is wainwright2019high $\VC(\Pi) = p+1$ for $\Pi = \cb{\pi_{v,c} : \pi_{v,c}(x) = 1\p{\cb{v \cdot x \geq c}}, \, v \in \RR^p, \, c \in \RR}$. Thus, our approach applies to linear decision rules in dimension $p_n \leq n^\beta$ for some $\beta < 1/2$.
\paragraph{Decision Trees}
Trees represent decision rules recursively \citep*{breiman1984classification}. A depth-0 decision tree $T_0$ is a trivial decision rule, $T_0(x) = a$ for some $a \in \cb{0, \, 1}$ and all $x \in \xx$. For any $L \geq 1$, a depth-$L$ decision tree $T_L$ is specified via a splitting variable $j \in 1, \, ..., \, p$, a threshold $t \in \RR$, and two depth-$(L-1)$ decision trees $T_{(L-1),A}$ and $T_{(L-1),B}$, such that $T_L(x) = T_{(L-1),A}(x)$ if $x_j \leq t$, and $T(x) = T_{(L-1),B}(x)$ else. See Figure (ref) for an example of a decision tree. The class of depth-$L$ decision trees over $\RR^p$ has VC dimension bounded on the order of $\VC(\Pi) = \too\p{2^L \log(p)}$.\footnote{This bound follows Lemma 4 of \citet*{zhou2018offline}, paired with the alternative characterization of the VC dimension given in Section A of the supplemental material. The notation $f(n) = \too(g(n))$ means that there is a function $h(\cdot)$ that scales poly-logarithmically in its argument for which $f(n) \leq h(g(n))g(n)$.} Thus, our results apply to trees whose depth may grow as $L_n = \lfloor \kappa \log_2(n) \rfloor$ for some $\kappa < 1/2$.
\paragraph{Monotone Rules}
We have $x \in [0, \, 1]^2$ and units get treated if $x_2$ exceeds some increasing function of $x_1$, i.e., $\Pi = \cb{\pi_f : \pi_f(x) = 1\p{\cb{x_2 \geq f(x_1)}}, \, \text{$f$ is monotone increasing}}$. This class has infinite VC dimension, because any set of points $\cb{x_i}_{i = 1}^d$ with $x_i = (\alpha_i, \, \alpha_i^2)$ and $0 < \alpha_1 < \ldots < \alpha_d < 1$ can be shattered using $\Pi$. Thus, our results do not apply to monotone rules over $[0, \, 1]^2$.\footnote{The difficulty here is not a mere technicality: Monotone decision rules can match arbitrary decision rules along the curve $(\alpha, \, \alpha^2)$ for $\alpha \in [0, \, 1]$, and so it is impossible to establish any non-trivial learning rates over monotone decision rules without making further assumptions on the distribution of the features $X_i$. In particular, we need assumptions that guarantee that all observations cannot concentrate around the curve $(\alpha, \, \alpha^2)$. In this paper, we do not consider results that require specific distributional assumptions over the features $X_i$. We note however the recent work by mbakop2016model, who establish polynomial rates of convergence for learning monotone rules under an assumption that the $X_i$ have a bounded density under Lebesgue measure on $[0, \, 1]^2$.}
We are now ready to state our main result on the asymptotic regret of policy learning using doubly robust scores. Following chernozhukov2016double,chernozhukov2016locally we assume that we run our method with scores obtained via cross-fitting, which is a type of data splitting that can be used to verify asymptotic normality given only high-level conditions on the predictive accuracy of the methods used to estimate nuisance components. In particular, cross-fitting allows for the use of black-box machine learning tools provided we can verify that they are accurate in mean-squared error as in Assumption (ref).
We proceed as follows: First divide the data into $K$ evenly-sized folds and, for each fold $k = 1, \, ..., \, K$, run an estimator of our choice on the other $K - 1$ data folds to estimate the functions \smash{$m_n(x, \, w)$} and \smash{$g_n(x, \, z)$}; denote the resulting estimates \smash{$\hatm_n^{(-k)}(x, \, w)$} and \smash{$\hg_n^{(-k)}(x, \, z)$}. Throughout, we will only assume that these nuisance estimates are accurate in the sense of Assumption (ref). Then, given these pre-computed values, we choose \smash{$\hpi_n$} by maximizing a doubly robust estimate of $A(\pi) = 2V(\pi) - \EE{\tau(X_i)}$,
where $k(i) \in \cb{1, \, ..., \, K}$ denotes the fold containing the $i$-th observation. The $K$-fold algorithmic structure used in (ref) was proposed in an early paper by schick1986asymptotically as a general purpose tool for efficient estimation in semiparametric models, and has also been used by other authors including robins2017minimax and zheng2011cross.
Finally, we assume that the weighting function $g_n(x, \, z)$ is bounded uniformly as below. In the case of a binary exogenous treatment, this is equivalent to the “overlap” assumption in the causal inference literature imbens2015causal, whereby $\eta \leq \PP{W_i = 1 \cond X_i = x} \leq 1 - \eta$ for all values of $x$. In our setting, the condition below acts as a generalization of the overlap assumption hirshberg2018balancing.
We also define the following quantities, where $S_n$ bounds the second moment of the scores, and $S_n^*$ is the asymptotic variance for estimating the policy improvement $A(\pi)$ of the best policy in $\Pi_n$ via (ref):\footnote{By expanding the square, we see that policies with higher values have lower variance of their scores, and so $S_n^*$ corresponds to the asymptotic variance for evaluating an optimal policy. Moreover, in the case where arguments from newey1994asymptotic imply that the doubly robust estimator (ref) is efficient, then $S_n^*$ is the semiparametric efficient variance for evaluating an optimal policy.}
We note that, unless we have an exceptionally large signal-to-noise ratio, we will have $S_n^* \geq S_n / 4$ and so the rounded log-term in (ref) below is just 0. A proof of Theorem (ref) is given in the following section.
In the simplest case where the maximizer of \smash{$\hA_n(\pi)$} over $\pi \in \Pi_n$ is unique and $\psi_n = 0$ (i.e., we solve the maximization problem exactly), the statement in (ref) simplifies to a bound on \smash{$\EE{R_n\p{\hpi_n}}$}, where \smash{$\hpi_n$} is as defined in (ref). However, in practice, \smash{$\hA_n(\pi)$} may have many maximizers. Moreover, the optimization problem (ref) is not convex and so---given a reasonable computational budget---we may only be able to solve it to within some tolerance $\psi_n > 0$. The more comprehensive form of our result given above highlights the fact that, in this case, our regret bound in fact applies uniformly over all approximate solutions to (ref).
In this section, we present a series of results that culminate in a proof of Theorem (ref), given in Section (ref). All other proofs are deferred to Section C of the supplemental material. Recall that we study policy learning for a class of problems where regret can be written as in (ref) using a function $V_n(\pi) = \EE[n]{\pi(X_i)\tau_n(X_i)}$, and we obtain \smash{$\hpi_n$} by maximizing a cross-fitted doubly robust estimate of \smash{$A_n(\pi) = 2V_n(\pi) - \EE[n]{\tau_n(X_i)}$} defined in (ref) over the class $\Pi_n$. If we could use $\hA_n(\pi) = A_n(\pi)$, then (ref) would directly yield the regret-minimizing policy in the class $\Pi_n$; but of course we never know $A_n(\pi)$ in applications. Thus, the main focus of our formal results is to study stochastic fluctuations of the empirical process \smash{$\hA_n(\pi) - A_n(\pi)$} for $\pi \in \Pi_n$, and examine how they affect the quality of policies learned via (ref).
We start our analysis by characterizing concentration of an ideal version of the objective in (ref) based on the true influence scores $\Gamma_i$, rather than doubly robust estimates thereof:
The advantage of studying concentration of the empirical process \smash{$\tA_n(\pi) - A_n(\pi)$} over the set $\pi \in \Pi_n$ is that it allows us, for the time being, to abstract away from the estimation tools used to obtain \smash{$\hA_n(\pi)$}, and instead to focus on the complexity of empirical maximization over the class $\Pi_n$.
A convenient way to bound the supremum of this empirical process over any class $\Pi$ is by controlling its Rademacher complexity $\rr_n(\Pi)$, defined as\footnote{Note that, conditionally on $\cb{X_i, \, \Gamma_i}_{i = 1}^n$ and the Rademacher variables $\xi_i$, the sum $\sum_{i = 1}^n \xi_i \Gamma_i \p{2\pi(X_i) - 1}$ can only take $2^n$ distinct values. Thus, the definition of $\rr_n(\Pi)$ does not entail any measure theoretic problems.}
where the $\xi_i$ are independent Rademacher (i.e., sign) random variables $\xi_i = \pm 1$ with probability $1/2$ each bartlett2002rademacher. For intuition as to why Rademacher complexity is a natural complexity measure, note that $\rr_n(\Pi)$ characterizes the maximum (weighted) in-sample classification accuracy on randomly generated labels $\xi_i$ over classifiers $\pi \in \Pi$; thus, $\rr_n(\Pi)$ measures how much we can overfit to random coin flips using $\Pi$.
Following this proof strategy, we bound the Rademacher complexity of “slices” of our policy class $\Pi_n$, defined as
The reason we focus on slices of $\Pi_n$ is that, when we use doubly robust scores, low-regret policies can generally be evaluated more accurately than high-regret policies, and using this fact allows for sharper bounds. Specifically, we can check that \smash{$n\text{Var}[\tA_n(\pi)] = S_n - A_n^2(\pi)$}, and so
where $S_n$ and $S_n^*$ are defined in (ref). This type of slicing technique is common in the literature, and has been used in different contexts by, e.g., \citet*{bartlett2005local} and gine2006concentration.
The following result provides such a bound in terms of the second moments of the doubly robust score, specifically $S_n^\lambda$ and $S_n$. This bound is substantially stronger than corresponding bounds used in existing results on policy learning. kitagawa2015should build their result on bounds that depend on $\max\cb{\Gamma_i}/\sqrt{n}$, which can only be used with scores that are uniformly bounded in order to get optimal rates. Meanwhile, bounds that scale as $\sqrt{S_n^\lambda\log(n)/n}$ are developed by \citet*{cortes2010learning}, maurer2009empirical and swaminathan2015batch; however, the additional $\log(n)$ factor makes these bounds inappropriate for asymptotic analysis.
Then, following the well known approach of bartlett2002rademacher, we use our bound on Rademacher complexity to obtain a uniform concentration bound for \smash{$\tA_n(\pi)$}. We use a refinement of the argument of bartlett2002rademacher based on Talagrand's inequality to get a bound that depends on second moments of $\Gamma_i$ rather than $\sup \abs{\Gamma_i}$.
In our final argument, we will apply Corollary (ref) for different $\lambda$-slices, and verify that we can in fact focus on those slices where $\lambda$ is nearly 0. Before that, however, we also need to control the discrepancy between the feasible objective \smash{$\hA_n(\pi)$} and the oracle surrogate \smash{$\tA_n(\pi)$} studied here.
In the previous section, we established risk bounds that would hold if we could optimize the infeasible value function \smash{$\tA_n(\pi)$}; we next need to extend these bounds to cover the situation where we optimize a feasible value function. As discussed above, we focus on the doubly robust estimator (ref), obtained using cross-fitting as in chernozhukov2016double,chernozhukov2016locally. As preliminaries, we note that the results of chernozhukov2016locally immediately imply that, given Assumption (ref), \smash{$\hA_n(1)$} is an asymptotically normal estimate of $A_n(1)$, where we use “1” as shorthand for the “always treat” policy. Furthermore, it is easy to check that given any fixed policy $\pi$,
meaning that the discrepancy between the two value estimates decays faster than the variance of either.
However, in our setting, the analyst gets to optimize over all policies $\pi \in \Pi_n$, and so coupling results established for a single pre-determined policy $\pi$ are not strong enough. The following lemma extends the work of chernozhukov2016locally to the case where we seek to establish a coupling of the form (ref) that holds simultaneously for all $\pi \in \Pi_n$.
The above result is perhaps surprisingly strong: Provided that the dimension $\VC(\Pi_n)$ of $\Pi_n$ does not grow too fast with $n$, the bound (ref) is the same coupling bound as we might expect to obtain for a single policy $\pi$, and the dimension of the class $\Pi_n$ does not affect the leading-order constants in the bound. In other words, in terms of the coupling of \smash{$\tA_n(\pi)$} and \smash{$\hA_n(\pi)$}, we do not lose anything by scanning over a continuum of policies $\pi \in \Pi_n$ rather than just considering a single policy $\pi$.
The doubly robust form used here is not the only way to construct efficient estimators for the value of a single policy $\pi$---for example, \citet*{hirano2003efficient} show that inverse-propensity weighting with non-parametrically estimated propensity scores may also be efficient---but it plays a key role in the proof of Lemma (ref). In particular, under Assumption (ref), the natural bound for the bias term due to misspecification of the nuisance components in fact holds simultaneously for all $\pi \in \Pi$, and this helps us pay a smaller-than-expected price for seeking a uniform result as in (ref). It is far from obvious that other efficient methods for evaluating a single policy $\pi$, such as that of hirano2003efficient, would lead to equally strong uniform couplings over the whole class $\Pi_n$.
Given that Assumption (ref), (ref), (ref) and (ref) hold with parameters $\beta < \min\cb{\zeta_m, \, \zeta_g}$, a combination of results from Corollary (ref) and Lemma (ref) implies that \smash{$\hA_n(\cdot)$} concentrates around $A_n(\cdot)$ over $\Pi_n^\lambda$. To conclude, it now remains to apply these bounds at two different values of $\lambda$. First we choose $\lambda^* > 0$ such as to satisfy $4(\lambda^*)^2 + 4\lambda^*\sup\cb{A(\pi) : \pi \in \Pi_n} \leq S_n^*$, so that the following holds via (ref): $$ S_n^{\lambda^*} + 4(\lambda^*)^2 \leq S_n^* + 4(\lambda^*)^2 + 4\lambda^*\sup\cb{A(\pi) : \pi \in \Pi_n} \leq 2S_n^*. $$ Then, by Corollary (ref) and Lemma (ref), we find that the limsup of the following expression is bounded by 1 as $n$ goes to infinity:
Now, recall that if any two functions $h(\cdot)$ and \smash{$\hat{h}(\cdot)$} are uniformly coupled as \smash{$|h(u) - \hat{h}(u)| \leq b$} for all $u \in U$ and \smash{$\hat{h}(\hat{u}) \geq \sup\{\hat{h}(u) : u \in U\} - \psi$}, then $$h(\hat{u}) \geq \hat{h}(\hat{u}) - b \geq \hat{h}(u) - b - \psi \geq h(u) - 2b - \psi$$ for any $u \in U$. Thus, the above implies that (recall that $A_n(\pi)$ scales with $2R_n(\pi)$)
and we note that $\psi_n$ decays fast enough by assumption that it can be omitted from (ref) without altering the result. In other words, if we knew that our learned policy approximately maximizes \smash{$\hA_n(\pi)$} and has regret less than $\lambda^*$, then we could guarantee that its regret decays at the desired rate.
To prove our result, it remains to show that all approximate maximizers of $\hA_n(\cdot)$ have regret bounded by $\lambda^*$ enough for (ref) to capture the leading-order behavior of regret. To do so, we apply a similar argument as above, but at a different value of $\lambda$. Consider \smash{$\lambda_+ = 3\limsup_{n \rightarrow \infty} \sup\cb{R_n(\pi) : \pi \in \Pi_n}$}, and by (ref) we see that
Furthermore, note that \smash{$\Pi_n^{\lambda_+} = \Pi_n$} for large enough $n$, and so (ref) in fact also holds with \smash{$\Pi_n^{\lambda_+}$} replaced by $\Pi_n$. Meanwhile, from (ref) paired with Markov's inequality we know that
By combining these two bounds, we see that
and moreover, because $\tau_{m_n}(x, \, w)$ is uniformly bounded, we find that the contribution of events where (ref) fails to hold to (ref) is vanishingly small as $n$ gets large.
To complement the upper bounds given in Theorem (ref), we also present lower bounds on the minimax risk for policy learning. Our goal is to show that our bounds are the best possible regret bounds that flexibly account for the distribution of the observed data and depend on the policy class $\Pi$ through the Vapnik-Chervonenkis dimension $\VC(\Pi)$. For simplicity, we here only consider the case where $W_i$ is binary and unconfounded; lower bounds for other cases considered in this paper can be derived via analogous arguments.
To establish our result, we consider lower bounds over sequences of problems defined as follows. Let $\xx_s := [0, \, 1]^s$ denote the $s$-dimensional unit cube for some positive integer $s$, and let $f(x)$ and $e(x)$ be $\lceil s/2 + 1 \rceil$ times continuously differentiable functions over $\xx_s$. Moreover, let $\sigma^2(x)$ and $\tau(x)$ be functions on $\xx_s$ such that $\sigma^2(x)$ is bounded away from 0 and $\infty$, and $\abs{\tau(x)}$ is bounded away from $\infty$. Then, we define an asymptotically ambiguous problem sequence as one where $\cb{X_i, \, Y_i, \, W_i}$ are independently and identically distributed drawn as
Because of the number of derivatives assumed on $f(x)$ and $e(x)$, it is well known that simple series estimators satisfy Assumption (ref).\footnote{See nickl2007bracketing for an argument that holds for arbitrary distributions $\pp$ supported on $[0, \, 1]^s$. We also note that, for a complete argument, one needs to address the fact that we have not assumed the treatment effect function $\tau(x)$ to be differentiable. To address this issue, note that in our data-generating process (ref) we have $\EE{Y_i|X_i= x} = f(x)$ regardless of $n$. Thus, because both $e(x)$ and $f(x)$ are sufficiently differentiable, we can use standard results about series estimation to obtain $o_P(n^{-1/4})$-consistent estimators \smash{$\he(x)$} and \smash{$\hat{f}(x)$} for these quantities. Next, for the purpose of our policy learner, we simply set \smash{$\hatm(x, 0) = \hatm(x, 1) = \hat{f}(x)$}; and because $\EE{\tau^2(X_i)/\sqrt{n}} = \oo(1/n)$, these regression adjustments in fact satisfy Assumption (ref).} Thus, because the magnitude of the treatment effects shrinks in (ref), \smash{$S_n^*$} and $S_n$ both converge to $S_{\pp}$ as defined below, and so Theorem (ref) immediately implies that, under unconfoundedness,
for any policy class $\Pi$ with finite VC dimension. The following result shows that (ref) is sharp up to a universal constant (whose value is less than 200).\footnote{The strategy of proving lower bounds relative to an adversarial feature distribution $\pp$ is standard in the machine learning literature; see, e.g., devroye1995lower. If we fix the distribution $\pp$ a-priori, then regret bounds for empirical risk minimization over $\Pi$ based on structural summaries of $\Pi$ (such as the VC dimension) may be loose bartlett2006empirical; however, it is not clear how to exploit this fact other than by conducting ad-hoc analyses for specific choices of $\Pi$.}
Here, the fact that we focus on problems where the magnitude of the treatment effect scales as $1/\sqrt{n}$ is important, and closely mirrors the type of asymptotics used by hirano2009asymptotics. If treatment effects decay faster than $1/\sqrt{n}$, then learning better-than-random policies is effectively impossible---but this does not matter, because of course all decision rules have regret decaying as $o(1/\sqrt{n})$ and so Theorem (ref) is loose. Conversely, if treatment effects dominate the $1/\sqrt{n}$ scale, then in large samples it is all but obvious who should be treated and who should not, and it is possible to get regret bounds that decay at superefficient rates luedtke2017faster, again making Theorem (ref) loose. But if the treatment effects obey the $\Theta(1/\sqrt{n})$ scaling of hirano2009asymptotics, then the problem of learning good policies is neither trivial nor impossible, and the value of using doubly robust policy evaluation for policy learning becomes apparent.
Finally, we note that the bounds of kitagawa2015should for inverse-propensity weighting are not asymptotically sharp in the above sense. Even when propensity scores are known, kitagawa2015should assume that $\abs{Y_i} \leq M$ and $\eta \leq e(X_i) \leq 1 - \eta$, and then prove regret bounds that scale as \smash{$M/\eta \sqrt{\VC(\Pi)/n}$} instead of \smash{$\sqrt{S_{\pp}\VC(\Pi)/n}$} in (ref). Now, the bound of kitagawa2015should is of course sometimes sharp, e.g., it is optimal if all we know is that $\abs{Y_i} \leq M$ and $\eta \leq e(X_i) \leq 1 - \eta$, but it is not adaptively sharp for asymptotically ambiguous sequences of problems as in (ref). In particular, the ratio of the upper bound of kitagawa2015should and the lower bound (ref) scales as $M/(\eta \sqrt{S_{\pp}})$, and there exist sequences of type (ref) where this ratio may be arbitrarily large.\footnote{Using the techniques developed in this paper, we can sharpen the bounds of kitagawa2015should and asymptotically replace $M/\eta$ by \smash{$\EE{Y_i^2/(e(X_i)(1 - e(X_i)))}^{1/2}$}. However, even this improved bound may exceed (ref) by an arbitrarily large factor.}
We now illustrate the value of doubly robust scoring techniques for policy learning using both an example from program evaluation and simulation studies. In Section (ref) we revisit a randomized evaluation of California's GAIN program, while Section (ref) presents a simulation study with endogenous treatment assignment. We present additional simulation results on nudge interventions to a continuous treatment variable in Section B of the supplemental material.
Recall that our approach to policy involves a 3-step algorithm. We start with a set of $n$ independent and identically distributed training examples $(X_i, \, Y_i, \, W_i, \, Z_i)$ and a class $\Pi$ of acceptable policies. Then, we
The main points of freedom left to the analysts involve the choice of estimator for $m(\cdot)$ and $g(\cdot)$ in Step 1, and the implementation of the optimization problem in Step 3. We emphasize that the choice of estimator for $m(x, \, w)$ and $g(x, \, z)$ in Step 1 and the choice of policy class $\Pi$ along with the optimizer used in Step 3 can be made fully independently.
For Theorem (ref) to apply, the main requirement on the method used to estimate $m(x, \, w)$ and $g(x, \, z)$ in Step 1 is that its error decays fast enough in mean-squared error, as detailed in Assumption (ref). Here, one option is to use non-parametric estimators for which we can precisely spell out when they satisfy Assumption (ref), such as sieve-based methods \citep*{chen2007large} or kernel regression \citep*{caponnetto2007optimal}; another is to use more heuristic methods from the statistical learning literature, such as boosting, random forests, or neural networks, in the hope that they will empirically be more accurate in finite samples than sieve or kernel-based methods.\footnote{In a recent advance, \citet*{farrell2018deep} established conditions under which deep neural networks can be shown to provably satisfy the conditions required by Assumption (ref). Thus, depending on the statistical setting and the chosen architecture, deep neural networks could either be seen as a formally validated alternative to sieve-type methods or as heuristic method.} One possible compromise is to run both classical methods known to satisfy Assumption (ref) asymptotically and heuristic statistical learning tools, and then synthesize the output of all models via cross-validation. As argued in \citet*{van2007super}, this approach essentially matches the finite-sample accuracy of the best method under consideration while preserving the asymptotic guarantees of the classical ones.
Meanwhile, the optimization problem in Step 3 is not a convex optimization problem, and so solving it can be computationally challenging. Several authors, including beygelzimer2009offset, kitagawa2015should, \citet*{zhang2012estimating} and \citet*{zhao2012estimating}, have noted that this optimization problem is numerically equivalent to a weighted classification problem,
where we train a classifier $\pi(\cdot)$ with response $H_i$ using sample weights $\lambda_i$. Given this formalism, we can build on existing tools for weighted classification to learn $\hpi$; see \citet*{zhou2018offline} for a further discussion.\footnote{Some popular approaches for solving problems of the form (ref) include best-subset empirical risk minimization chen2016best and optimal trees bertsimas2017optimal. Due to the computational difficulty of solving the problem (ref) exactly, it may also be of interest to consider the empirical performance of alternative methods that solve an approximation to our weighted classification problem, e.g., support vector machines cortes1995support or recursive partitioning \citep*{breiman1984classification}. However, we caution that our formal results only apply to methods that solve the problem (ref) exactly; see wager2020regression for further discussion.} In all our experiments, we set $\Pi$ to be a class of finite-depth decision trees (see Section (ref) for a definition), and solve the optimization problem in Step 3 using our companion R-package policytree sverdrup2020policytree,CRAN; see zhou2018offline for further details and motivation behind the computational strategy taken in this package.
The Greater Avenues for Independence (GAIN) program, started in 1986, is a welfare-to-work program that provides participants with a mix of educational resources and job search assistance. Between 1988 and 1993, the Manpower Development Research Corporation conducted a randomized study to evaluate the program. As described in \citet*{hotz2006evaluating}, randomly chosen registrants were eligible to receive GAIN benefits immediately, whereas others were embargoed from the program until 1993. All experimental subjects were followed for a 9-year post-randomization period and, as documented by hotz2006evaluating, eligibility for GAIN had a significant impact on mean quarterly income averaged over this 9-year period.
Our current question is whether we can find ways to prioritize treatment to some subgroups of GAIN registrants particularly likely to benefit from it. We consider data from four counties, Alameda, Riverside, Los Angeles and San Diego, resulting in $n = 19,170$ observations, and use $p = 28$ covariates, including demographics, education, and per-quarter earnings for 10 quarters preceding treatment. As in hotz2006evaluating, we use average quarterly income over the 9-year post-randomization period (in \$1000s) as our outcome.
Each county participating in the GAIN evaluation conducted its own randomized controlled trial, and the counties had considerable freedom in how they carried out the randomization. In particular, counties had flexibility in choosing whom to enroll in the randomized trial, and which fraction of participants to randomize into treatment. The data reflects this heterogeneity in study specifications: The per-county average outcome for controls varied from 0.64 to 1.04 thousand dollars per quarter, while the per-county fraction of treated units varied from 0.50 to 0.86.
We use this dataset to design a semi-synthetic observational study by pooling the data from all four counties under consideration. Because the mean control outcome and treatment fraction vary from county to county (and are in fact correlated), we expect that an uncorrected analysis of the pooled data would suffer from confounding. In an attempt to correct for the confounding that arises from pooling we pursue a selection-on-observables strategy, and assume that controlling for the $p = 28$ covariates described above is enough to correct for the different study specifications used in different counties.
Our method starts by computing doubly robust scores for the treatment effect, and learning policies by empirical maximization as in (ref). We use the augmented inverse-propensity weighted scores of \citet*{robins1994estimation}, with nuisance component estimates from generalized random forests \citep*{athey2018generalized,breiman2001random},\footnote{The one major deviation between how we compute scores below and the assumptions of Theorem (ref) is that, here, we use leave-one-out (or out-of-bag) estimates for $\tau(X_i)$, etc., whereas Theorem (ref) assumed $K$-fold estimation. The reason for this choice is that, as discussed in breiman2001random, random forests are particularly well suited for leave-one-out estimation, and allow the analyst to obtain such estimates at essentially no additional computational cost.}
where \smash{$\hat{f}(x)$} and \smash{$\he(x)$} are random forest estimates of $\mathbb{E}[Y_i \cond X_i= x]$ and $\mathbb{E}[W_i \cond X_i= x]$ respectively, \smash{$\htau(\cdot)$} is an causal forest\footnote{Random forests are a type of adaptive nearest neighbor estimator that use an ensemble of trees to define a relevant neighborhood function for each query point; see athey2018generalized for a discussion. Causal forests use the adaptive neighborhood function implied by a forest to fit a partially linear model using the method of robinson1988root; see nie2017learning for formal results motivating the use of local partially linear modeling for heterogeneous treatment effect estimation, and Section 1.3 of athey2019estimating for a discussion of how this partially linear modeling is carried out in causal forests. We emphasize that, for our purposes, random forests are simply used as a convenient non-parametric estimator of relevant nuisance components, specifically $f(x)$ and $e(x)$ here, and could seamlessly be replaced with other methods such as boosting or neural networks. The shape of the learned policy $\hpi$ is determined in the optimization step 3, which only depends on the random forests through the predictions used to form doubly robust scores \smash{$\hGamma_i$}.} estimate of the conditional average treatment effect, and $C$ is a parameter measuring the cost of treatment. Tuning parameters for all forests were selected by leave-one-out cross-validation.\footnote{The regression surfaces \smash{$\hat{f}(x)$} and \smash{$\he(x)$} were tuned to optimize mean-squared error. As advocated in nie2017learning, the conditional average treatment effect function was tuned to optimize the error of a local residual-on-residual regression.} Here, we set $C = 0.14$ to roughly match the average treatment effect with the goal of ensuring that the optimal treatment rule is not trivial (i.e., we can only achieve non-zero utility gains by exploiting treatment heterogeneity).
Before starting to optimize policies we first run a brief sanity check on our selection-on-observables strategy, and confirm the ability of estimators that build on this assumption to accurately recover the average treatment effect we would get using a proper randomization-based estimator that does not pool data across counties. The natural doubly robust estimator of the average treatment effect in our setting is \smash{$\htheta_{DR} = \sum_{i = 1}^n \hGamma_i / n$}, with scores \smash{$\hGamma_i$} as in (ref). We compare it to a naive difference-in-means estimator $\htheta_{DM} = \text{avg}\cb{Y_i : W_i = 1} - \text{avg}\cb{Y_i : W_i = 0}$ that does not attempt to correct for bias due to pooling, and to an “oracle” doubly robust estimator that does not estimate propensity scores from covariates but instead uses the true per-county treated fractions: \smash{$\htheta_{DR}^* = \sum_{i = 1}^n \hGamma_i^* / n$} with
where $G_i \in \cb{\text{Alameda}, \, \text{Riverside}, \, \text{Los Angeles}, \, \text{San Diego}}$ denotes the county-membership of the $i$-th sample. Because \smash{$\htheta_{DR}^*$} uses the true per-county treatment fractions $\he_i^*$ and estimates nuisance components using cross-fitting, the point estimates will be $\sqrt{n}$-consistent and the associated confidence intervals asymptotically valid essentially without assumptions \citep*{rothe2018flexible,wager2016high}. The resulting point estimates for the average treatment effect ($\pm 1$ standard error) are: $\htheta_{DR} = 0.141 \pm 0.026$ for the feasible doubly-robust estimator, $\htheta_{DR}^* = 0.146 \pm 0.028$ for the oracle doubly-robust estimator, and $\htheta_{DM} = 0.208 \pm 0.028$ for the naive difference in means. Thus, it appears that pooling county information results in confounding, but that controlling for available covariates helps.
We now move to learning a policy $\hpi$. In doing so, however, we note that caution is warranted because we have measured features pertaining to race, ethnicity, age and gender. On the one hand, there may be legal restrictions on the use of these features for treatment allocation but, on the other hand, they appear to act as counfounders. For example, as shown in Table (ref), white GAIN registrants were randomized to treatment at higher rates than non-white registrants, and also white controls had higher outcomes than non-white controls. Our approach allows us to seamlessly use such sensitive variables for deconfounding without using them for policy allocation: We use these variables when estimating the nuisance components in (ref), but then omit them from the maximization step (ref) that produces the policy.
For our policy class $\Pi$, we consider decision trees of depth either 1 or 2. The learned decision rules are shown in Figure (ref). Interestingly, the depth-1 and 2 trees make the same decisions for the roughly 3/4 of GAIN registrants who were paid 3 quarters prior to randomization, but the depth-2 tree chooses to switch to a different rule for those who weren't paid 3 quarters prior.
In order to choose tree depth and, more broadly, to evaluate the accuracy of the policy learning procedure, we recommend cross-validation. We randomly divide the data into $K$ folds $\set_k$, $k = 1, \, ..., \, K$ and, for each fold, learn a policy $\hpi^{(-k)}(\cdot)$ using all but the data in $\set_k$. Here, we use $K = 10$. Finally, we estimate improvement over a random baseline as
Table (ref) shows the estimated improvement of our depth-1 and -2 trees, as well as two baselines: A variant of the inverse-propensity weighted method of kitagawa2015should using the propensities used to construct (ref), as well as a plug-in policy that does not obey our functional form restriction, and simply treats all samples with \smash{$\htau^{(-i)}(X_i) > C$}. Our depth-2 trees achieve markedly better performance than the depth-1 trees. Interestingly, the depth-2 tree is also competitive with the unconstrained plug-in estimator. Based on this analysis, we prefer the depth-2 tree in Figure (ref).
One potential concern with this analysis is that our evaluation hinges on validity of the selection-on-observables assumption, as well as accuracy of the doubly robust scores \smash{$\hGamma_i$} from (ref). To assuage this concern, we also computed a version of the improvement measure (ref), but with scores \smash{$\hGamma_i^*$} computed using the true per-county treatment fractions as in (ref):
As seen in the second rightmost column of Table (ref), our feasible evaluation discussed above gave the correct ordering for the methods, but was somewhat optimistic in terms of the quality of the learned policies. The formal properties of treatment rules whose complexity is tuned via cross-validation are considered by mbakop2016model.\footnote{Recall that cross-validation is a means of evaluating the quality of the policy learning procedure, not the decision that was produced by a specific realization of the procedure. If we want an accuracy assessment that is valid conditionally on the learned rule \smash{$\hpi(\cdot)$}, one can either use a single test-train split, or use the more sophisticated data carving approach of \citet*{fithian2014optimal}.}
In order to develop a richer quantitative understanding of the behavior of our method, we now turn to a simulation study. Here, we consider a setting with a binary, endogenous treatment $W_i$ and a binary instrument $Z_i$ and assume homogeneity as in (ref). In this case, our method chooses the policy \smash{$\hpi = \argmax\{\frac{1}{n} \sum_{i = 1}^n \p{2\pi(X_i) - 1} \hGamma_i : \pi \in \Pi\}$}, where \smash{$\hGamma_i$} is a cross-fit doubly robust score with estimates of the compliance weights as in (ref):
\sloppy{where \smash{$\Delta(x) = \PP{W_i = 1 \cond Z_i = 1, \, X_i = x} - \PP{W_i = 1 \cond Z_i = 0, \, X_i = x}$} is the conditional average effect of the instrument on the treatment, \smash{$z(x) = \PP{Z_i = 1 \cond X_i = x}$}, \smash{$f(x) = \EE{Y_i \cond X_i = x}$}, \smash{$e(x) = \PP{W_i = 1 \cond X_i = x}$}, and $\tau(x)$ is the conditional average treatment effect as specified in (ref). We estimate all nuisance components via random forest methods with the package grf, and use an instrumental forest for $\tau(\cdot)$, a causal forest for $\Delta(\cdot)$, and a regression forest for $f(\cdot)$, $e(\cdot)$ and $z(\cdot)$.}
In this simulation experiment, we generate data independently as follows, for various choices of $n$ and $\tau(\cdot)$:
Note that $W$ is in fact endogenous, because $Q$ (and thus also $W$) is more likely to be 1 when the noise term $\varepsilon$ is large. Given this setup, we consider $\tau(\cdot)$ functions
In both cases, we learn $\pi(\cdot)$ over the class $\Pi$ of depth-2 trees and note that best non-parametric policy \smash{$\pi^*(x) = 1\p{\cb{\tau(x) > 0}}$} belongs to $\Pi$ in case (ref) but not in case (ref).
In Figure (ref), we display the improvement $A(\pi) = \EE{(2\pi(X_i) - 1)\tau(X_i)}$ of our learned policies relative to a random assignment baseline, for different values of $n$. Over all, we see that the regret of the learned policies improves with $n$, and approaches best-in-class regret as $n$ gets large. We also note an interesting difference in the behavior of the learned rules in settings (ref) and (ref). In the first case, $\tau(\cdot)$ is continuous, and regret improves smoothly with sample size. Conversely, in the second case where $\tau(\cdot)$ has sharp jumps, we observe something of a phase transition between $n = 2,000$ and $n = 4,000$, as our trees become able to consistently make splits that roughly match the jumps in $\tau(\cdot)$.
In this paper, we proposed an approach to policy learning in the observational study setting that builds on classical ideas for semiparametrically efficient treatment effect estimation. Our main result is that doubly robust estimators of average treatment effects can be adapted for policy evaluation, and that the policy that maximizes the resulting doubly robust value estimate over a pre-specified class $\Pi$ satisfies rate-optimal guarantees for minimax regret. Our approach decouples estimation of nuisance components used for the doubly robust scores from optimization of the doubly robust value function, and thus allows practitioners flexibility in how they implement each step.
Our formal discussion focused on regret bounds for policy learning. A natural follow-up question is to ask for confidence sets guaranteed to contain an optimal policy: For example, if $\Pi$ is the set of depth-$L$ decision trees, can we identify a subset of $\Pi$ guaranteed to contain a value-maximizing policy in $\Pi$ with high probability? Some early results in this direction are reported by rai2018statistical. Meanwhile, armstrong2015inference consider the related task of identifying a subset of the population we are confident will benefit from the policy intervention.
Another natural direction to extend our results is towards dynamic decision making problems, where the policy maker needs to make a sequence of decisions, potentially depending on time-varying covariates. The problem of doubly robust policy evaluation in this setting has been considered by \citet*{thomas2016data} and \citet*{zhang2013robust}. \citet*{nie2019learning} proposed a method for learning observational stopping rules from observational data that is both computationally feasible and robust to confounding. Obtaining a more comprehensive landscape of the problem of dynamic policy learning in observational studies would be of considerable interest.
Finally, all results presented here relied on point-identification of treatment effects, either via a selection on observables assumption or via an instrument that satisfies conditional homogeneity. Some applications, however, do not allow for such clean assumptions, and thus call for methods for policy learning that are robust to failures of identifying assumptions. kallus2018confounding consider the problem of policy learning under an approximate selection-on-observables assumption in the sense of rosenbaum2002observational. It would also be of interest to study what can be done if we only have access to a monotone instrument, as in manski2000monotone.
{0.2pt plus 0.3ex}