Fisher-Schultz Lecture: Generic Machine Learning Inference on Heterogenous Treatment Effects in Randomized Experiments, with an Application to Immunization in India

Fisher-Schultz Lecture: Generic Machine Learning Inference on Heterogenous Treatment Effects in Randomized Experiments, with an Application to Immunization in India

Introduction

Randomized Controlled Trials (RCT) and Machine Learning (ML) are arguably two of the most important developments in data analysis methods for applied researchers. RCTs play an important role in the evaluation of social and economic programs, medical treatments and marketing duflo2007using,imbens2015causal. ML is a name attached to a variety of constantly evolving statistical learning methods including Random Forest, Boosted Trees, Neural Networks, Penalized Regression, Ensembles, and Hybrids; see, e.g., wasserman:ML, ESL, bishop2006pattern, murphy2012machine, hastie2015, goodfellow2016deep and james2021 for prominent textbook treatments. ML has become a key tool for prediction and pattern recognition problems, surpassing classical methods in high dimensional settings.

At first blush, those two sets of methods may seem to have very different applications: in the most basic randomized controlled experiment, there is a sample with a single treatment and a single outcome. Covariates are not necessary and even linear regression is not the best way to analyze the data freeman2008,imbens2015causal. In practice, however, applied researchers are often confronted with more complex experiments. For example, there might be accidental imbalances in the sample, which require selecting control variables in a principled way. ML tools, such as the lasso method proposed in BelloniChernozhukovHansen2011,bcfh17 or the double machine learning method proposed in DML, have proven useful for this purpose. Moreover, some complex RCT designs have so many treatment combinations that ML methods may be useful to select the few treatments that actually work and pool the rest with the control groups for statistical power banerjee2019improving. Finally, researchers and policy makers are often interested in features of the impact of the treatment that go beyond the simple average treatment effect. In particular, very often, they want to know whether the treatment effect depends on covariates, such as gender, age, etc. This heterogeneity is essential to assess if the impact of the program would generalize to a population with different characteristics, and, for economists, to better understand the driving mechanism behind the effects of a particular program. In a review of 189 RCTs published in top economic journals since 2006, we found that 76 (40%) report at least one subgroup analysis, wherein they report treatment effects in subgroups formed by baseline covariates.\footnote{The papers were published in Quarterly Journal of of Economics, American Economic Review, Review of Economics Studies, Econometrica and Journal of Political Economy. We thank Karthik Mularidharan, Mauricio Romero and Kaspar W\"{u}thrich for sharing the list of papers they computed for another project.}

One issue with reporting treatment effects split by subgroups, however, is that there might be a large number of potential ways to form subgroups. Often researchers collect rich baseline surveys, which give them access to a large number of covariates: choosing subgroups ex-post opens the possibility of overfitting. To solve this problem, medical journals and the FDA require pre-registering the sub-sample of interest in medical trials in advance. In economics, this approach has gained some traction with the adoption of pre-analysis plans, which can be filed in the AEA registry for randomized experiments. However, restricting the heterogeneity analysis to pre-registered subgroups amounts to throwing away a large amount of potentially valuable information, especially now that many researchers collect large baseline data sets. It should be possible to use the data to discover ex post whether there is any relevant heterogeneity in treatment effect by covariates.

To do this in a disciplined fashion and avoid the risk of overfitting, scholars have recently proposed using ML tools. Indeed, ML tools seem ideal for exploring heterogeneity of treatment effects when researchers have access to a potentially large array of baseline variables to form subgroups and few guiding principles on which of those are likely to be relevant. Several recent papers, which we review below, develop methods for detecting heterogeneity in treatment effects. Empirical researchers have taken notice.\footnote{In the recent past, several new empirical papers in economics used ML methods to estimate heterogeneous effects. E.g. Rigol showed that villagers outperform the machine learning tools when they predict heterogeneity in returns to capital. Heller predicted who benefits the most from summer internship projects. Deryugina used the methods developed in the present paper to evaluate the heterogeneity in the effect of air pollution on mortality. Creponetal also built on the present paper to develop a methodology to determine if the impact of two different programs can be accounted for by different selection. The methodological papers reviewed later also contain a number of empirical applications.}

This paper develops a generic approach to using any of the available ML tools to predict and make inference on heterogeneous treatment or policy effects. A core difficulty of applying ML tools to the estimation of heterogenous causal effects is that, while they are successful in prediction empirically, it is much more difficult to obtain uniformly valid inference, i.e., inference that remains valid under a large class of data generating processes. In fact, in high dimensional settings, absent strong assumptions, generic ML tools may not even produce consistent estimators of the conditional average treatment effect (CATE), the difference in the expected potential outcomes between treated and control states conditional on covariates. Previous attempts to solve this problem focused either on specific tools (for example, the method proposed by athey:trees, which has become popular with applied researchers, and uses trees), or on situations where those assumptions might be satisfied. Our approach to resolving the fundamental impossibilities in non-parametric inference is different. Motivated by wasserman:adaptive, instead of attempting to get consistent estimation and uniformly valid inference on the CATE itself, we focus on providing valid estimation and inference on features of CATE.

We start by building a ML proxy predictor of CATE, and then target features of the CATE based on this proxy predictor. In particular, we consider three objects, which are likely to be of interest to applied researchers and policy makers: (1) Best Linear Predictor (BLP) of the CATE on the ML proxy predictor; (2) Sorted Group Average Treatment Effects (GATES) or average treatment effect by heterogeneity groups induced by the ML proxy predictor; and (3) Classification Analysis (CLAN) or the average characteristics of the most and least affected units defined in terms of the ML proxy predictor. Thus, we can find out if there is detectable heterogeneity in the treatment effect based on observables, and if there is any, what the treatment effect is for different bins. And finally we can describe which of the covariates are associated with this heterogeneity.

There is a trade-off between more restrictive assumptions or tools and a more ambitious estimation. We address this trade-off by focusing on coarser objects of the function rather than the function itself, but make as little assumptions as possible. This seems to be a worthwhile sacrifice: the objects for which we have developed inference appear to us at this point to be the most relevant, but in the future, one could easily use the same approach to develop methods to estimate other objects of interest. For example, Creponetal used the same technique to construct and estimate a specific form of heterogeneity at the post-processing stage. Even then, as we will see, getting robust and conservative standard errors for heterogeneity requires a larger sample size than just estimating average treatment effects. This reflects a different trade-off: if we do not assume that we can predict ex ante where the heterogeneity might be (in which case we can write it down in a pre-analysis plan), power will be lower, and detecting heterogeneity will require a larger sample. This is a consideration that applied researchers will need to keep in mind when designing and powering their experiments, and when writing pre-analysis plans: if a particular dimension of heterogeneity is deemed important, it should be pre-specified.

Another trade-off between our generic approach and the alternative approaches in the literature is that the parameters of interest are estimated conditional on a particular split of the data, since our estimation and inference methods rely on sample splitting to avoid overfitting and other inferential non-regularities. Conditional on a single data split into a training and a hold-out sample, statistical inference is conceptually straightforward and appealing. Indeed, in this case, statistical inference reduces to the classical inference for linear regression and sample means. Theoretically, if a researcher can credibly pre-commit to a single data split, this gives one clean solution to the inferential problem. However, this introduces additional variability in the results, stemming from the random draw of the particular split of the data. Therefore, researchers often consider multiple sample splits, and aggregate the results. This reduces the probability that two researchers working with the same data will arrive at different conclusions. To formalize and support this approach, we propose quantile-aggregated inference -- which aggregates inferential results by taking medians of estimates and medians and other quantiles of upper and lower confidence intervals obtained from different splits. We show that quantile aggregation formally lowers estimation (reporting) risks over a single-split procedure, and we establish its inferential properties.

The proposed approach is generic in that it can be applied in conjunction with any ML method. To compare and select among ML methods, we develop goodness-of-fit measures for the BLP and GATES. We also take one step backward and use these goodness-of-fit measures to build ML proxies that better target the CATE through causal learning. We show that these causal machines produce provably better proxies of the CATE than generic (predictive) ML methods. Moreover, by designing the ML to target CATE directly, the post-processing methods that we develop can focus on providing valid inference, rather than correcting biases.

We apply our method to a large-scale RCT of nudges to encourage immunization in the state of Haryana, Northern India. This experiment, an important practical application in its own right, is designed and discussed in banerjee2019improving. Immunization is generally recognized as one of the most effective and cost-effective ways to prevent illness, disability, and diseases. Yet, worldwide, close to 20 million children every year do not receive critical immunizations \citep*{unicef2019progress}. While early policy efforts have focused mainly on improving the infrastructure for immunization services, a more recent literature suggests that “nudges” (such as small incentives, leveraging the social network, SMS reminders, social signalling, etc.) may have large effects on the use of those services.\footnote{See, for example, banerjee2010improving,bassani2013financial,wakadha2013feasibility, johri2015strategies,oyo2016interventions,gibson2017mobile,karing2018social,domek2016sms,uddin2016use,regan2017randomized,alatas2019celebrities,banerjee2019leveraging. } This project was a collaboration with the government of Haryana, which was willing to experiment with a combination of nudges, with the goal of choosing the most effective policy and implement it at scale. It built a custom vaccination platform, and ran a large-scale experiment covering seven districts, 140 Primary health centers, 2,360 villages involved in the experiment (including 915 at risk for all the treatments), and 295,038 children in the resulting database. Immunization was very low at baseline: in every single village of the district, the fraction of children whose parents had reported they received the measles vaccine (the last in the sequence) was 39%, and only 19.4% had received the vaccine before the age of 15 months, whereas the full sequence is supposed to be completed in one year. The experiment was a village-level, cross randomized design of three main nudges: providing incentives, sending SMS reminders, and seeding ambassadors. It included several variants for each policy: the level and schedule of the incentives, the number of people receiving reminders, and the mode of selection of the ambassadors, leading to a large number (75) of finely differentiated bundles.

banerjee2019improving developed a methodology to identify the most effective and cost-effective bundle of policies, based on an application of LASSO to a marginal effects specification that imposes some structure on the bundles, and in particular, the idea that policy variants (e.g. level of incentives, or level of coverage of SMS reminders) may be indistinguishable in practice. They found that the most cost-effective policy is to combine “information hubs” (people identified by others as good at diffusing information) and SMS reminders. This is cheap and can be done everywhere. In fact, they showed that this policy is the only one among those tested that would actually save money to the government for each measles shot, while increasing immunization. But the most effective policy, i.e., the policy that increases immunization the most, is the combination of incentives, immunization ambassadors, and SMS reminders, which is much more expensive. Yet, while this policy increases the cost per immunization, the effects are important: the number of monthly measles shots (the last vaccine in the schedule, and thus a marker for full immunization) delivered increases by 3.26, corresponding to 44% of the mean vaccination rate in the control group that got neither SMS nor increasing incentives and information hubs. The government was therefore interested in finding out where the program would be most effective, to implement it only in those places even at the higher cost per immunization.

The pre-analysis plan specified to look for heterogeneity by gender and by “Village-level baseline/national census variables, including assets, beliefs, knowledge, and attitudes towards immunization” but did not identify one or two specific baseline variables to look at. This reflected genuine uncertainty (as is often the case). Many factors can influence policy impact, from attitudes to implementation capabilities to baseline levels, and we did not have a specific theory of where to look. It is precisely the type of context that requires a principled approach to avoid overfitting, and provide a policy-relevant recommendation.\footnote{This approach of finding the best treatment and then looking at where it works the best gets closer to the idea of “personalized medicine”. Using the same data, agarwal2020synthetic go one step further and use a “synthetic intervention” approach to look for the policy that works the best for each kind of village. }

The rest of the paper is organized as follows. Section (ref) formalizes the framework, describes our approach and compares it with the existing literature. Section (ref) presents identification and estimation strategies for the key features of CATE of interest. Section (ref) introduces our inference method that accounts for uncertainty coming from parameter estimation and sample splitting. Section (ref) presents the construction of causal machines that can learn CATE better than purely predictive approaches or some existing proposals for causal approaches. Section (ref) reports the results of the empirical application and provides detailed implementation algorithms. Section (ref) concludes with some remarks. The Online Appendix (OA) gathers proofs of the main theoretical results and additional technical results.

Our Agnostic Approach

This section present our framework and approach. We observe $\mathrm{Data} := (Y_i, Z_i, D_i)_{i=1}^N$, consisting of i.i.d. copies of the random vector $(Y,Z,D)$ having probability law $P$, where $Y$ is the outcome of interest, $D$ is a binary treatment indicator, and $Z$ is a possibly high-dimensional vector of covariates that characterize the observational units. The data is defined on an underlying probability space with measure ${\mathrm{P}}$. The expectation operator is denoted by ${\mathrm{E}}$. When we need to emphasize the dependence of ${\mathrm{P}}$ and ${\mathrm{E}}$ on $P$, we use the notation ${\mathrm{P}}_P$ and ${\mathrm{E}}_P$.

Model and Key Causal Functions.

Let $Y(1)$ and $Y(0)$ be the potential outcomes in the treatment state 1 and the non-treatment state 0 neyman1923applications,rubin74. The main causal functions are the baseline conditional average (BCA):

and the conditional average treatment effect (CATE):

Suppose $D$ is randomly assigned conditional on $Z$, with probability of assignment depending on a subvector of stratifying variables $Z_1 \subseteq Z$, namely

and the propensity score is known and is given by

which we assume is bounded away from zero or one:

This setup is similar to rosenbaum1983.

The observed outcome is $Y = DY(1) + (1-D) Y(0)$. Under the stated assumption, the causal functions are identified by the components of the regression function of $Y$ given $D, Z$:

that is, $\ b_0(Z) = {\mathrm{E}} [Y \mid D=0, Z]$, and

This regression underlies the use of predictive ML methods that learn ${\mathrm{E}}[Y \mid D, Z]$ and then estimate CATE using the formula.

Alternatively one can identify CATE using the following two equivalent “causal" regressions:

where $H$ is the residualized treatment scaled by its variance:

also known as the Horvitz-Thompson transform. We mention these alternative strategies here, because as shown in Section 5, they can lead to better ways of approximating $s_0(Z)$ than through the predictive regression ((ref)), and our inference tools equally apply to ML methods that try to learn $s_0(Z)$ through either of these relations. In fact, in our empirical analysis, the strategies based on (ref) measurably outperform the strategies based on (ref).

Estimation and Inference Challenges

Regardless of the way we try to learn $s_0(Z)$, estimation and inference are challenging in modern high-dimensional settings, because the target function $z \mapsto s_0(z)$ can live in a very complex class. ML methods effectively explore various forms of sparsity to yield “good" approximations to $s_0(z)$. In its simplest form, sparsity reduces the complexity of $z \mapsto s_0(z)$ by assuming that it can be well-approximated by a function that only depends on a low-dimensional subset of $z$, making consistent estimation possible. As a result, these methods can perform much better than classical methods in high-dimensional settings under sparsity. However, sparsity or, more generally, low complexity of the CATE function $s_0$, are untestable assumptions that must be used with caution.

Without some form of sparsity, it is hard, if not impossible, to obtain consistent estimators of $z \mapsto s_0(z)$. There are several fundamental reasons as well as large gaps between theory and practice that are responsible for this. One fundamental reason is that ML methods might not even produce consistent estimators of $z \mapsto s_0(z)$ in high dimensional settings. For example, if $z$ has dimension $d$ and the target function $z \mapsto s_0(z)$ is assumed to have $p$ continuous and bounded derivatives, then the worst case (minimax) lower bound on the rate of learning this function from a random sample of size $N$ cannot be better than $N^{-p/(2p + d)}$ as $N \to \infty$, as shown by stone82. Hence if $p$ is fixed and $d$ is also small, but slowly increasing with $N$, such as $d \geqslant \log N$, then there exists no consistent estimator of $z \mapsto s_0(z)$ generally. Hence, generic ML estimators cannot be regarded as consistent, unless further assumptions are made. Examples of such assumptions include structured forms of linear and non-linear sparsity and super-smoothness.\footnote{The function $z \mapsto s_0(z)$ is super-smooth if it has continuous and bounded derivatives of all orders.} The problem of obtaining uniformly valid inference on $z \mapsto s_0(z)$ using generic ML methods is even more difficult.\footnote{While the previous assumptions make consistent adaptive estimation possible BickelRitovTsybakov2009, confidence sets that adapt to unknown regularity (smoothness or sparsity) do not exist even for low-dimensional nonparametric problems low1997,wasserman:adaptive. Let $z \mapsto s_0(z)$ be a target function that lives in an infinite-dimensional class with unknown regularity $s$ (e.g., smoothness or degree of sparsity). Adaptive consistent estimation (resp. inference) for $z \mapsto s_0(z)$ with respect to $s$ is possible if there exists a consistent estimator (resp. valid confidence set) with a rate of convergence (resp. diameter) that changes with $s$ in a (nearly) rate-optimal way. Construction of adaptive confidence bands then requires making additional untestable assumptions. See, e.g., gine:nickl, where self-similarity conditions are used in low-dimensional nonparametric problems.}

In this paper, we take an agnostic view. We neither rely on any sparsity or low-complexity assumptions to make the ML estimators consistent, nor impose other stronger conditions to make “traditional" confidence intervals valid. We simply treat ML as providing proxy predictors for the objects of interest.

Our Approach

To address the previous challenges, we propose strategies for estimation and inference on key features of $s_0(Z)$ rather than on $s_0(Z)$ itself. Because of this difference in focus, and by relying on sample splitting, we can loosen the restrictions about the properties of the ML estimators.

Let $(M,A)$ denote a random partition of the set of indices $\{1, \ldots, N\}$. The strategies that we consider rely on random splitting of $\mathrm{Data} = (Y_i, D_i, Z_i)_{i=1}^N$ into a main sample, denoted by $\mathrm{Data}_M$ $=$ $(Y_i, D_i, Z_i)_{i \in M}$, and an auxiliary sample, denoted by $\mathrm{Data}_A = (Y_i, D_i, Z_i)_{i \in A}$. We will sometimes refer to these samples as $M$ and $A$. After splitting the sample, we carry out two stages:

Stage 1: From the auxiliary sample $A$, we obtain ML estimators of the baseline functions and treatment effects, which we call the ML proxy predictors, $$ z \mapsto B(z) \text{ and } z \mapsto S(z). $$ Here $S(Z)$ is a possibly biased and noisy predictor of $s_0(z)$ and $B(Z)$ is a possibly biased and noisy predictor of $b_0(Z)$ (or other technical “baseline" functions, as we discuss in Section (ref)). We do not require these predictors to be consistent for the true functions.

Stage 2: We post-process the proxies from Stage 1 to estimate and make inference on features of the CATE function $z \mapsto s_0(z)$ in the main sample $M$. The key features that we target include:

Our approach is generic with respect to the ML method being used, and is agnostic about its formal properties. However, it relies on sample splitting, which introduces a specific source of uncertainty. To account for this, we use many data splits into main and auxiliary samples to produce robust estimators, and we employ quantile aggregation of inference to combine results across splits. Specifically, for point estimation, we report the median of the estimated key features over different random splits of the data. We take medians and other quantiles of many random conditional confidence sets for interval estimation. Finally, we construct p-values by taking medians of many random conditional p-values. We establish the formal inferential properties of this procedure.

Relationship to the Literature.

We focus the review strictly on the literatures about estimation and inference on heterogeneous effects and inference using sample splitting.

This work is related to the literature that uses linear and semiparametric regression methods for estimation and inference on heterogeneous effects. crump2008nonparametric developed tests of treatment effect homogeneity for low-dimensional settings based on traditional series estimators of the CATE. A semiparametric inference method for characterizing heterogeneity, called the sorted effects method, was given in CFL2014. This approach does provide a full set of inference tools, including simultaneous bands for percentiles of the CATE, but is strictly limited to the traditional semiparametric estimators of the regression and causal functions. hansen:kobzur proposed a sparsity-based method called “targeted undersmoothing" to perform inference on heterogeneous effects. This approach does allow for high-dimensional settings, but imposes sparsity as well as additional assumptions that enable the targeted undersmoothing. A related approach, which allows for simultaneous inference on many coefficients (for example, inference on the coefficients corresponding to the interaction of the treatment with other variables) is proposed BCK-LAD using a Z-estimation framework, where the number of interactions can be very large; see also dezeure2016high for a more recent effort in this direction, focusing on de-biased lasso in mean regression problems. This approach, however, still relies on a strong form of sparsity assumptions. zhao2017selective proposed a post-selection inference framework within high-dimensional linear sparse models for the heterogeneous effects. The approach is attractive because it allows for some misspecification of the model.

Another approach is to use tree-based and other methods. Imai:stuff discussed the use of a heuristic support-vector-machine method with lasso penalization for classification of heterogeneous treatments into positive and negative ones. They used the Horvitz-Thompson transformation of the outcome hirano:imbens:ridder,abadieDD such that the new outcome becomes an unbiased, noisy version of CATE.\footnote{Note that using Horvitz-Thompson (HT) transform of outcome, in this and other references, typically gives very noisy signal. One can improve the approach by either including the HT transform interacted with some baseline covariates as regressors in a regression model, as we do in the present paper, or using residualized outcomes in conjunction with HT, as in semenova2020.} athey:trees made use of the Horvitz-Thompson transformation of the outcome to inform the process of building causal trees, with the main goal of predicting CATE. They also provided a valid inference result on average treatment effects for groups defined by the tree leaves, conditional on the data split into two subsamples: one used to build the tree leaves and the one to estimate the predicted values given the leaves. Like our methods, this approach is essentially assumption-free. Our paper is a complement, in that our approach can be used with any ML method. wager:athey proposed a subsampling-based construction of a causal random forest, providing valid pointwise inference for CATE (see also the review in wager:athey on prior uses of random forests in causal settings) for the case when covariates are very low-dimensional (and essentially uniformly distributed).\footnote{The dimension $d$ is fixed in wager:athey; the analysis relies on the Stone's model with smoothness index $\beta=1$, in which no consistent estimator exists once $d \geqslant \log n$ in the minimax sense.} This condition rules out the typical high-dimensional settings that arise in many empirical problems, especially in current RCTs, where the number of baseline covariates is potentially very large.

Several other studies look at model-based strategies for performing inference on CATE. semenova2020 used ML to perform inference on the “partial" CATE, ${\mathrm{E}}[ s_0(Z) \mid X]$, where $X$ is a prespecified low-dimensional set of covariates. Specifically, they constructed an estimator of a denoised HT transform of the outcome and projected it using a nonparametric series estimator on the set of low-dimensional prespecified covariates of interest $X$, whose dimension is much lower than the dimension of $Z$.\footnote{Specifically, the denoised HT transform of outcome is $\tilde Y = g(1,Z) - g(0,Z) + H ( Y - g(D,Z) )$, where $H = (D-p(Z))/[p(Z)(1-p(Z)]$ and $g(D,Z) = {\mathrm{E}} (Y \mid D, Z)$. semenova2020 used ML to estimate $g(D,Z)$ and the propensity score $p(Z)$, in case the latter is unknown.} The main advantage of this approach is that it delivers familiar nonparametric inference on partial CATE (even though inference on the full CATE remains intractable). fan2022estimation, zimmert, CNS:AutoLocal, and CNS:simple developed kernel versions of this procedure. Related ideas but based on partialling-out (using residualized outcomes and treatment) appear in semenova:panel, semenova:ROML, nie:20, foster2019orthogonal, and kennedy2020optimal. Relative to the approach taken here, the assumptions made in these papers are more restrictive, but deliver stronger results. For example, in the approach of semenova2020, one has to specify the baseline covariates $X$ for the partial CATE analysis, (which is exactly what we are trying to avoid in our approach). Second, the methods critically rely on the consistency of ML to estimate the nuisance components well, which can be restrictive in high dimensional settings, as discussed above.

The idea of using a “hold out” sample to validate the result of a ML procedure to discover heterogeneity was suggested in Heller, who used the method proposed in wager:athey and compared their results to the heterogeneity in a holdout sample. Our inference approach is different because it calls for multiple splits. This procedure itself is also of independent interest and could be applied to many problems, where sample splitting is used to produce ML predictions abadie2013endogenous. Related references include wasserman2009high, and meinshausen2009p, where the ideas are related, but the details are quite different, as we shall explain below. The premise is the same; however, as in meinshausen2009p and rinaldo2016bootstrapping -- we should not rely on a single random split of the data and should adjust inference in some way. Our construction of p-values builds upon ideas in meinshausen2009p, though what we propose is simpler, and our confidence intervals appear to be new. Of course, sample splitting ideas are classical, going back to hartigan1969using,kish1974inference,barnard,cox1975note,mosteller:tukey, though having been mostly underdeveloped and overlooked for inference, as characterized by rinaldo2016bootstrapping. Finally, our inference method shares with the literature on post-selection inference in statistics that the target estimands are random functions depending on a ML proxy fithian2014optimal,lee2016exact.

Main Identification Results and Estimation Strategies

In this section we condition on $\mathrm{Data}_A$ (the auxiliary sample in a random split) and therefore consider the functions $$ z \mapsto B_A(z) := B(z; \mathrm{Data}_A) \text{ and } z \mapsto S_A(z) := S(z; \mathrm{Data}_A). $$ as fixed functions. To lighten the notation, we keep the dependence on $\mathrm{Data}_A$ implicit in all the expectations, objects and parameters, that is, e.g., we use ${\mathrm{E}}[\cdot]$ and $S(Z)$ instead of ${\mathrm{E}}[\cdot \mid \mathrm{Data}_A]$ and $S_A(Z)$. We shall make the dependence on $\mathrm{Data}_A$ explicit when we discuss estimation and inference in Section (ref).

Best Linear Predictor of CATE

The first inferential target is the best linear predictor of the CATE using the proxy $S(Z)$.

By construction, $\mathsf{BLP}_{}[s_0(Z) \mid S(Z)]$ is an unbiased predictor of $s_0(Z)$, which improves over $S(Z)$ in the mean-squared error sense, that is $$ {\mathrm{E}} \{ s_0(Z) - \mathsf{BLP}_{}[s_0(Z) \mid S(Z)]\}^2 \leqslant {\mathrm{E}} [ s_0(Z) - S(Z)]^2. $$ Indeed, we can quantify the improvement by $$ {\mathrm{E}} [ s_0(Z) - S(Z)]^2 - {\mathrm{E}} \{ s_0(Z) - \mathsf{BLP}_{}[s_0(Z) \mid S(Z)]\}^2 = (1-\beta_2)^2 \operatorname{Var}[S(Z)] + [{\mathrm{E}} S(Z) - {\mathrm{E}} s_0(Z)]^2, $$ which is positive unless $S(Z)$ is an unbiased predictor and, either $\beta_2 = 1$ or $\operatorname{Var}[S(Z)] = 0$.\footnote{The previous expression follows from the decompositions $ {\mathrm{E}} [ s_0(Z) - S(Z)]^2 = {\mathrm{E}} [ \{s_0(Z) - {\mathrm{E}} s_0(Z)\} - \{S(Z) - {\mathrm{E}} S(Z)\} + \{{\mathrm{E}} s_0(Z) - {\mathrm{E}} S(Z)\}]^2 $ and ${\mathrm{E}} \{ s_0(Z) - \mathsf{BLP}_{}[s_0(Z) \mid S(Z)]\}^2 = {\mathrm{E}} \{ [s_0(Z) - {\mathrm{E}} s_0(Z)] - [\mathsf{BLP}_{}[s_0(Z) \mid S(Z)] - {\mathrm{E}} s_0(Z)]\}^2,$ using that ${\mathrm{E}} \mathsf{BLP}_{}[s_0(Z) \mid S(Z)] = {\mathrm{E}} s_0(Z)$ and $\beta_2 = \mathrm{Cov}_{}(s_0(Z), S(Z))/\operatorname{Var}_{}(S(Z))$. } Accordingly, compared to the ML proxy, the BLP can be seen as a refined predictor of the individual CATE, $s_0(Z)$. If $S(Z)$ is a perfect proxy for $s_0(Z)$, then $\beta_2 =1.$ In general, $\beta_2 \neq 1$, correcting for noise in $S(Z)$. If $S(Z)$ is complete noise, uncorrelated to $s_0(Z)$, then $\beta_2 =0$. Furthermore, if there is no heterogeneity, that is, $s_0(Z) = s$, then $ \beta_2 = 0. $ Rejecting the hypothesis $\beta_2 = 0$ therefore means both that there is heterogeneity in $s_0(Z)$ and $S(Z)$ is a relevant predictor.

We provide two strategies for identifying and estimating $\mathsf{BLP}_{}[s_0(Z) \mid S(Z)]$.

Strategy A: Weighted Residual BLP

Consider the weighted linear projection:

where $S:=S(Z)$, $w(Z) := \{p(Z)(1-p(Z))\}^{-1},$ $X :=(X_1', X_2')',$ $$X_1 = [1, B(Z), p(Z), p(Z) S(Z)]', \quad X_2 :=[D- p(Z), (D- p(Z))(S - {\mathrm{E}}_{} S) ]'.$$ The term $B(Z)$ could be replaced by any “noise-reducing" proxy function. For example, the algorithms of semenova:panel and nie:20, targeting the delivery of $S(Z)$, also construct $B(Z)$ that are meant to approximate ${\mathrm{E}}[Y \mid Z]$ but not $b_0(Z)$; see also Section (ref) for other examples of such algorithms. We include $X_1$ to reduce finite sample noise of the estimators of the BLP parameter based on this strategy.\footnote{Note that $X_1$ can include other functions of $Z$. In our experiments, the use of $B(Z)$ strongly improves the precision of estimating BLP (and other quantities such as GATEs and CLAN introduced below).}

Note that $\alpha_1$ and $\alpha_2$ are identified under weak assumptions. Further, we note that the interaction $(D- p(Z)) (S - {\mathrm{E}}_{} S)$ is {orthogonal} to $D- p(Z)$ under the weight $w(Z)$, and to all functions of $Z$ such as $X_1$. Consequently, we obtain the following result that shows that the linear projection (ref) identifies the BLP.

The identification result in Theorem (ref) is constructive. We can base a corresponding estimation strategy on the empirical analog:

where $X_i = [X_{1i}',X_{2i}']'$, $X_{1i} := [1, B(Z_i), p(Z_i), p(Z_i) S(Z_i)]'$, $X_{2i} :=[D_i- p(Z_i), (D_i- p(Z_i))(S_i - \Bbb{E}_{N,M} S_i) ]'$, and $\Bbb{E}_{N,M}$ denotes the empirical expectation with respect to the main sample, i.e. $$\Bbb{E}_{N,M} h(Y_i, D_i, Z_i) := |M|^{-1}\sum_{i \in M} h(Y_i, D_i, Z_i) .$$

Figure (ref) provides two examples. The left panel shows a case where $s_0(Z) = 0$ with zero effect and zero heterogeneity in the CATE, whereas the right panel shows a case where $s_0(Z) =Z$ with strong heterogeneity in the CATE. In both cases, we evenly split 1,000 observations between the auxiliary and main samples, $Z$ follows uniform distribution on $(-1,1)$, $b_0(Z) = 3Z$, $U$ is standard normal, independently of $Z$, and $Y$ is generated by (ref). We obtain the proxy predictor $S(Z)$ by Breiman's random forest, using the ranger implementation in R ranger17.

In the first example, the ML proxy is pure noise by construction, and the BLP post-processor correctly eliminates the noise, producing a CATE prediction that is roughly a 68% better approximation to the CATE under the RMSE metric. Furthermore, using our inferential methods of Section 4, we cannot reject the null hypothesis that the BLP is zero. In the second case, under strong heterogeneity, the signal in the ML proxy dominates the noise component. As a result, the BLP does not change the ML proxy drastically, but still gives a meaningful improvement to the CATE under the RMSE metric. These improvements agree with the theoretical arguments given above.\footnote{To show that the risk reductions are not a fluke, we repeated the calculations in $1,000$ simulations. We found average risk reductions of $65\%$ and $18\%$ in the no heterogeneity and strong heterogeneity examples, respectively.}

Strategy B: HT BLP

This strategy makes use of the Horvitz-Thompson transform $H$ defined in (ref). It is well known that the transformed response $YH$ provides an unbiased signal about CATE:

and it follows by the properties of the best linear predictor that $$ \mathsf{BLP}_{}[s_0(Z) \mid S(Z)] = \mathsf{BLP}_{} [YH \mid S(Z)].$$ Note that $\mathsf{BLP}_{}[s_0(Z) \mid S(Z)]$ is a more precise unbiased predictor of $s_0(Z)$ than $YH$ because, by construction $$ \operatorname{Var}(\mathsf{BLP}_{}[s_0(Z) \mid S(Z)]) = \operatorname{Var}(\mathsf{BLP}_{}[YH \mid S(Z)]) \leqslant \operatorname{Var}(YH). $$ The R-squared of $\mathsf{BLP}_{} [YH \mid S(Z)]$ quantifies the percent reduction in variance of the BLP relative to $YH$.

The simple linear projection $\mathsf{BLP}_{} [YH \mid S(Z)]$ is completely fine for identification purposes, but can severely underperform in estimation and inference due to lack of precision. We can repair the deficiencies by considering, instead, the linear projection:

where $\tilde X := (X_1' H, \tilde X_2 ')'$, $ \tilde X_2 := (1, S- {\mathrm{E}}_{} S)'$, and $X_1 := [1, B(Z), p(Z), p(Z) S(Z)]'$ as before. The term $X_1$ could contain other functions of $Z$. We include $X_1H$ in order to reduce noise.

The following theorem shows that the linear projection (ref) also identifies the BLP.

Theorem (ref) leads to an estimator defined through the empirical analog:

where $\tilde X_i := (X_{1i}' H_i, \tilde X_{2i} ')'$, $ \tilde X_{2i} := (1, S_i- \Bbb{E}_{N,M} S_i)'$, and $X_{1i} := [1, B(Z_i), p(Z_i), p(Z_i) S(Z_i)]'$

Sorted Group Average Treatment Effects

The second inferential target is the set of group average treatment effects, where the groups are induced by $S(Z)$.

We provide two strategies for identifying and estimating the GATES.

Strategy A: Weighted Residual GATES

Consider the weighted linear projection equation:

where $ W := (X_1', W_2')'$, $X_1$ contains a vector of functions of $Z$, e.g., $X_1 = (B(Z), p(Z) \{1(G_k)\}_{k=1}^K)'$ and $W_2 := (\{ [D-p(Z)] \cdot \{1(G_k)\}_{k=1}^K )'.$ The presence of $D-p(Z)$ in the interaction $[D- p(Z)] \cdot 1(G_k) $ orthogonalizes this regressor relative to all other regressors that are functions of $Z$, such as $X_1$. The controls in $X_1$, as in the BLP, are included to reduce noise in estimation.

Theorem (ref) below shows that the linear projection (ref) identifies the GATES. We can therefore base an estimation strategy on the empirical analog:

where $\widehat \alpha = (\widehat \alpha_1, \ldots, \widehat \alpha_K)'$.

Strategy B: HT GATES

Here we employ a linear projection on Horvitz-Thompson transformed variables:

where $ \tilde W := ( X_1' H, \tilde W_2')'$, $X_1$ includes functions of $Z$, e.g. $X_1$ the same as above, and $\tilde W_2 := [\{1(G_k)\}_{k=1}^K ]'.$

Theorem (ref) shows that the linear projection (ref) also identifies the GATES. We can therefore base an estimation strategy on the empirical analog:

where $\widehat \mu = (\widehat \mu_1, \ldots, \widehat \mu_K)'$. The resulting estimator has similar performance to the estimator in (ref), and under some conditions their first-order properties coincide.

We now provide a formal statement of the identification results.

Classification Analysis

When the BLP and GATES analyses reveal substantial heterogeneity, it is interesting to know the properties of the subpopulations that are the most and least affected. Here we focus on the “least affected group" $G_1$ and “most affected group" $G_K$, where the labels “most" and “least" can be swapped depending on the context.

The parameters $\delta_1$ and $\delta_K$ are identified because they are averages of variables that are directly observed. The CLAN quantifies the differences between the most and least affected groups and singles out the covariates that are associated with the heterogeneity in the CATE. The CLAN can be extended to comparisons of features other than averages, such as variances, covariances or distributions. In the empirical analysis, we estimate the CLAN parameters by taking averages in $M$:

using $G_{k,i} = 1\{S(Z_i) \in I_k\}$, where $I_k = [\ell_{k-1},\ell_{k})$ and $\ell_k$ is the $(k/K)$-quantile of $\{S_{i}\}_{i\in M}$.

Goodness of Fit Measures for Fitting CATE

In practical applications it is useful to have goodness-of-fit measures to guide the selection of ML proxies.

For the analysis based on the BLP of CATE, we propose to use:

Maximizing $\Lambda$ is equivalent to maximizing the correlation between the ML proxy predictor $S(Z)$ and the CATE $s_0(Z)$, or equivalent to maximizing the $R^2$ in the regression of $s_0(Z)$ on $S(Z)$. Therefore, an ML method that attains a higher $\Lambda$ is a preferred method.

Analogously, for the GATES analysis, we propose to use:

This is the part of variation of $s_0(z)$, ${\mathrm{E}} s_0(Z)^2$, explained by $\bar S(Z) = \sum_{k=1}^K \gamma_k 1(S(Z) \in I_k)$. Hence choosing the ML proxy $S(Z)$ to maximize $\bar \Lambda$ is equivalent to maximizing the $R^2$ in the regression of $s_0(Z)$ on $\bar S(Z)$ (without a constant). If the groups $G_k = \{S \in I_k\}$ have equal size, namely ${\mathrm{P}}(S(Z) \in I_k) = 1/K$ for each $k=1,..., K$, then $$ \bar \Lambda = \frac{1}{K}\sum_{k=1}^K \gamma^2_k. $$ Therefore, a ML method that attains a higher $\bar \Lambda$ is a preferred method. The empirical versions of the parameters above are:

The choice of the ML method using goodness-of-fit measures does not pose any additional inferential challenge, when there is clearly an ML method that dominates the others -- so that we select the best method with probability approaching one. This means that the inferential methods of the next section would not need any further adjustment. When this is not the case, there are several possibilities depending on the scientific reporting objectives. For example, suppose there are two (near) winners and we want to construct a $(1-\alpha)$-confidence set, as in our empirical analysis. Then in the spirit of sensitivity analysis, we can report the union of the $(1-\alpha)$-confidence sets. This ensures the inferential coverage guarantee of $1-\alpha$ continues to apply, if the reader of the empirical report chooses one or the other winner at random. On the other hand, the inferential guarantee needs to be discounted to $1-2\alpha$ via Bonferroni adjustment, if the reader of the empirical report chooses one or the other report depending on the empirical results themselves. We use the first approach because the readers of our empirical analysis are not likely to follow the latter approach.

Finally, if the data sets are big, we could use additional splitting to choose the best-performing ML method, before taking the resulting ML proxies to the main sample.\footnote{ We also refer to Section (ref) which discusses other, more exploratory ideas for building the best ML algorithms for targeting CATE already in the first stage. }

Split-Sample Robust Estimation and Inference Methods

Estimation and Inference: The Generic Targets

Let $\theta$ denote a generic target parameter or functional. For example,

Let $(a,m)$ denote a fixed partition of $\{1,\ldots,N\}$. In this section we make explicit when the probabilities and expectations are conditional on the auxiliary sample $a$, $$\mathrm{Data}_a := \{ (Y_i, D_i, X_i)\}_{i\in a},$$ and also index the ML proxies $B= B_a$ and $S= S_a$ and estimands $\theta=\theta_a$ by $a$ to make explicit the dependence on $\mathrm{Data}_a$. As we mention in Comment (ref), this dependence vanishes as the size of the auxiliary sample becomes large under suitable conditions.

Single Split

We begin the discussion of inference conditional on a single split of data induced by the partition $\{(a,m)\}$ of $\{1,...,N\}$ into sets of cardinality $(N-n, n)$. All of the examples admit an estimator $\widehat \theta_a$ that is approximately Gaussian, conditionally on $\mathrm{Data}_a$, namely as $(N-n, n) \to \infty$ and for any $z$,

We provide sufficient regularity conditions for ((ref)) in Lemma (ref), which may be of independent interest (in single-split inference context).\footnote{The conventional asymptotic normality of the OLS estimator does not imply validity of the single-split inference procedure automatically, because the distribution of the estimator is random. We need to guarantee that this random distribution converges in probability to a normal (non-random), which requires conditions for splits being "regular" with high probability. We refer to Lemma (ref) in the OA for the formal result.}

As a consequence of (ref), the standard p-values $$ p^+_a := 1-\Phi( \widehat \sigma^{-1}_a (\widehat \theta_a - \theta_0)), \ \ p^-_a := \Phi( \widehat \sigma^{-1}_a (\widehat \theta_a - \theta_0)), $$ for testing the null hypothesis $ \theta_a = \theta_0$ against the alternatives $\theta_a>\theta_0$ and $\theta_a<\theta_0$, respectively, are approximately uniform under the null, namely $$ {\mathrm{P}} (p^{\pm}_a < \alpha \mid \mathrm{Data}_a) = \alpha + o_P(1). $$ As another consequence of (ref), the standard confidence interval (CI) $$ [L_a, U_a] := [\widehat \theta_a \pm \Phi^{-1}(1-\alpha/2) \widehat \sigma_a ] $$ covers $\theta_a$ with approximate probability $1-\alpha$ conditional on $\mathrm{Data}_a$: $$ {\mathrm{P}} [\theta_a \in [L_a, U_a ] \mid \mathrm{Data}_a] = 1-\alpha -o_P(1). $$ Thus, we have straightforward inference conditional on a single data split.

Multiple Splits

In practice, researchers often prefer using multiple splits $(a,m)'s$ to reduce estimation risk and demonstrate that the results are robust to how they split the data. Therefore, we need a way to aggregate the results across different splits, and propose quantile aggregation methods and analyze their properties.

Different partitions of $[N]$ yield different estimands and estimators. To formalize this randomness, we consider $A$ as a uniform random variable taking values $a \in \mathcal{A}$, that is, $ A \sim U(\mathcal{A})$. Therefore, conditional on $\mathrm{Data}$, the estimand $\theta_{A}$ is a random variable. In what follows, we will mainly target our inference on the median value of $\theta_A$, instead of $\theta_a$, the estimand for a specific partition. Furthermore, different partitions yield different estimators $\widehat \theta_A$ and approximate distributions for these estimators. Therefore, conditional on $\mathrm{Data}$, estimators $\widehat \theta_A$, p-values $p_A$, and intervals $[L_A, U_A]$ are all random variables. We will use quantile aggregation methods to summarize them.

It is useful to recall some definitions of quantiles for discrete variables. For a random variable $X$ with law ${\mathrm{P}}_X$ and $k$ points of support, and index $u \in (0,1)$, the lower and upper $u$-quantiles are $ \mathrm{\underline{Q}}_u(X) := \inf\{ x \in \Bbb{R}: {\mathrm{P}}_X(X \leqslant x) \geqslant u \},$ and $ \mathrm{\overline{Q}}_u(X) : = \sup\{ x \in \Bbb{R}: {\mathrm{P}}_X(X \geqslant x) \geqslant 1-u)$, respectively. For $w_u := \lfloor u k \rfloor/(\lfloor u k \rfloor+ \lceil u k \rceil)$, we define $$\mathrm{{Q}}_u(X) : = w_u \mathrm{\underline{Q}}_u(X) + (1-w_u)\mathrm{\overline{Q}}_u(X)$$ as the central quantile.\footnote{For example, the quantile function in R uses this definition Rcite.} If $X$ is continuous, all three definitions coincide. To define upper, lower, and central medians, we use $u=1/2$ in the definitions above and $\mathrm{M}$ instead of $\mathrm{Q}_u$.

We now formally define the median-aggregated p-value.

Aggregation using lower median p-values was first proposed by meinshausen2009p in the context of split-sample hypothesis testing in linear regression with selection. Here we take the central medians, since they are more likely to behave like regular p-values.\footnote{For example, the sample lower median of $\{U,1-U\}$, $U \sim U(0,1)$, obeys $P(\underline{\mathrm{M}} < \alpha) = 2\alpha$ for $\alpha< 1/2$. In contrast, the central median obeys ${\mathrm{P}}({\mathrm{M}} < \alpha) < \alpha$ for $\alpha < 1/2$.}

We next define the quantile-aggregated point and interval estimators.

We can interpret these definitions as risk-reducing inferential summaries.

Analogous risk contraction properties hold for the mean aggregation, but we focus on medians for robustness reasons.

In what follows, we study the formal inferential guarantees of $[L,U]$. The default choice of $\beta$ is $1/2$, but we obtain useful theoretical guarantees for other choices $\beta< 1/2$ as well.

Principal Regularity Condition

As the main regularity condition, we assume approximate normality of the split-sample t-statistics:\footnote{Here and below we use the standard error $\widehat \sigma_A$ instead of the theoretical standard deviation $\sigma_A$ in all statements, but we can exchange the two if $\widehat \sigma_A/\sigma_A \to_P 1$, which holds under typical conditions, e.g. Lemma (ref).}

Suppose that the data $\{(Y_i, Z_i, D_i)\}_{i=1}^N$ are generated as i.i.d. copies of $(Y, Z, D)$. In this case, for any $a \in \mathcal{A}$: $$ {\mathrm{P}}\{ \widehat \sigma^{-1}_A (\widehat \theta_A - \theta_A) < z \} = {\mathrm{E}} \left [\frac{1}{|\mathcal{A}|} \sum_{a \in \mathcal{A}} 1(\widehat \sigma^{-1}_a (\widehat \theta_a - \theta_a) < z) \right] = {\mathrm{P}}\{ \widehat \sigma^{-1}_a (\widehat \theta_a - \theta_a) < z \}, $$ because the expression on the right does not depend on $a \in \mathcal{A}$ under the i.i.d. sampling. This observation simplifies the verification of (R1) for least squares type estimators; see Lemma (ref). While the i.i.d. case is our main focus, the main results in this section rely only on the conditions labeled as R, which are likely to hold more generally.

Below we give various theoretical guarantees for our inferential summaries using this condition and adding more conditions to get stronger results.

Hypothesis Testing with Multiple Splits

We start the analysis by testing homogeneous hypotheses $\theta_A = \theta_0$, which imply that $\theta_a$ does not vary with $a$. Suppose, for example, that we want to test that the slope of the BLP is zero with probability one, $\beta_{2A} = 0, $ against the alternative $\beta_{2A}>0$ with positive probability. This problem amounts to both testing the heterogeneity in CATE and the relevance of the ML score $S_A$ as a predictor. Another interesting hypothesis is whether $\beta_{2A} = 1,$ with probability one, that is, whether $S_A$ is well-calibrated and needs no post-processing.

More generally, suppose we are testing the hypothesis

with probability one, against $H^+_1: \theta_A>\theta_0$ with positive probability. Testing using the median p-value $p^+$ will have power against the null when the majority of $\theta_a$'s violate the null, so we can interpret the rejection accordingly. Similarly, we can test against $H^-_1: \theta_A<\theta_0$ or $H_1: \theta_A \neq \theta_0$ with positive probability using the median p-values $p^-$ or $\bar p$.

Below we establish the properties of the median $p$-values under (R1). To get the sharpest results, we can invoke a concentration condition for approximate medians:

This condition states that the approximate median over-the-splits t-statistics tend to concentrate more than any single-split t-statistic. This condition is intuitive, but it is hard to give general primitive conditions for it.\footnote{When the $t$-stats are independent, then their median concentrates in a fixed interval around 1/2 with probability approaching 1 exponentially fast; see romano:SS. Therefore, the approximate median concentration condition holds. This happens when $m$'s are non-overlapping and some further homogeneity conditions hold. On the other hand, suppose the $p$-values are the same asymptotically, then the inequality in the concentration condition binds, but does not fail. This situation is not common in our context, though.} We believe (ref) is quite plausible. We verified it for typical values of $\alpha< 1/4$ numerically in various experiments that mimic empirical applications, and were unable to find any counterexample.

Therefore under the median concentration condition, the median p-values have the standard property. Without the median concentration condition, the median p-values need to be multiplied by 2. However, based on our computational experiments, median p-values are conservative even for the nominal level $\alpha$ (once $\alpha< 1/4$), mainly due to the concentration property holding with $\gamma''_N <0$. We therefore do not recommend multiplying by 2; see also romano:SS for a similar point. The exact form of $o(1)$, given in the proof, allows us to convert the results into those holding uniformly in a set of probability measures ${\mathrm{P}}$. The proof of the first result partly relies on the idea of meinshausen2009p to use Markov inequality to bound quantiles of an arbitrary collection of marginally uniform random variables.

Prediction Intervals with Multiple Splits

Outside of the settings with homogeneity, the estimand $\theta_A$ is a random variable, and we might be interested in characterizing its typical values. Our first approach serves this purpose, and is connected to conformal/permutation inference that is commonly used for predicting unobserved outcomes hoeffding1952large,vovketal,candes2022.

Here our goal is to have a prediction interval for $\theta_A$, and the challenge we face is that $\theta_A$ is not observed directly, which places us outside the standard conformal setting. However, for each $a \in \mathcal{A}$, we have a (random) confidence interval $[L_a, U_a]$ that has the covering property:

This condition is implied by the basic regularity condition (R1) in our context.

We take the quantile-aggregated confidence interval $[L,U]$ as our prediction interval for $\theta_A$.

We can use the prediction interval $[L,U]$ to characterize the “majority" of the central values of the random target $\theta_A$ that one could get from sample splitting. For this we set $\beta=.25$ and "small" $\alpha=o(1)$, then $MP=[L,U]$ has the property:

That is, $MP$ contains majority of central values of $\theta_A$. On the other hand, if we set $\beta =1/2$, we get a median aggregated interval $$MI=[L,U].$$ In this setting, we can think of $MI$ as predicting the median of $\theta_A$ with small margin of error, if $\alpha=o(1)$:

The latter gives us a useful interpretation of median confidence intervals, and notably this interpretation applies under the weakest possible regularity condition (R1) in our setting.

Confidence Intervals for Median Parameter with Multiple Splits

Instead of predicting the “majority" of $\theta_A$, we may focus the inference on a single target.

The choice of the target has the intuitive appeal of representing a typical value. Moreover, in many cases, $\theta_A$ will concentrate around its median value $\theta^*$, making it an even more natural target. What follows is the principal regularity condition that covers this concentration scenario.

The concentration condition above is a convenient property for the interpretability of the inference. Here the rate of concentration of $\theta_A$ around $\theta^*$ should be faster than the rate $\widehat \sigma_A$ of $\widehat \theta_A$ estimating $\theta_A$. Thus, the concentration condition implicitly requires the auxiliary set $a$ to be large and the main set $m$ to be small compared to $a$; so that randomness in the inferential target is small compared to the size of the estimation error $\widehat \sigma_a$ in the main sample, which typically is proportional to $n^{-1/2}$. Otherwise, if $n$ is large relative to $N-n$, then each $\theta_a$ would be estimated accurately, but there would be a lot of variation of $\theta_a$ across $a \in \mathcal{A}$.\footnote{As a practical diagnostic, we recommend the researchers to report the variation of $\widehat \theta_a$ across sample splits, in addition to the confidence interval for $\theta^*$. We are grateful to Guido Imbens for this suggestion.} Condition (R3) is high-level; we demonstrate the plausibility of this condition for the BLP parameter in Appendix (ref) of the OA using notions of estimation and algorithmic stability. \footnote{Estimation stability implies that $\theta_A$ concentrates around a fixed value $\theta_\bullet$, in which case the median also concentrates around $\theta_\bullet$. Estimation stability follows from the ML proxy $S_A$ being consistent for some fixed proxy function $s_\bullet$, but not necessarily consistent for the true CATE. This condition can be readily verified using statistical learning theory, as we do in Section (ref) for causal learners of CATE. The algorithmic stability condition is strictly weaker than estimation stability, though not as readily available. Our use of these stability criteria is inspired by similar ideas in wager:PNAS, CWZ:JASA, and vasilis:stable, applied to a different context. Appendix (ref) discusses all of this further.}

The following results summarize the properties of the proposed median confidence interval under various conditions.

Under the strongest assumptions, the target $\theta^*$ is covered with a probability of at least $1-\alpha-o(1)$. Under the minimal set of assumptions, the coverage probability is $1- 2 \alpha - o(1)$. In our numerical results, the confidence intervals tend to be conservative even under the minimal condition, with coverage exceeding $1- \alpha$. Therefore, using $1-\alpha$ as the nominal level is our recommended choice based on this evidence.

Further Consideration: Causal Machines that Learn CATE Better

Our main proposal so far is to take proxies from any first stage black-box machine and post-process them to better target CATE and perform inference on functionals of CATE, such as the BLP and GATES. But, can we design the machines to target CATE directly in the first stage? If we can, then the post-processing methods of the previous section would mostly focus on providing inference, and less on correcting biases of the first stage inputs. Building on athey:trees, we propose two types of such causal machines, and connect them to the emerging literature on orthogonal machine learning, such as nie:20, semenova:panel, and foster2019orthogonal, among others.

Focusing ML Methods on CATE in Stage 1

We propose two options, taking ideas from our stage 2 analysis to stage 1. Specifically, we can train ML proxies in the auxiliary sample based on either:

where $B(Z)$ is now a technical “baseline" function of covariates $Z$, as described below, whose role is to reduce noise in the learning problem.

We can refer to the first causal learner as the weighted residual (WR) learner, and the second causal learner as the HT learner. Examples of parameter spaces include spaces of linear functions generated by a set of dictionary transformations of $Z$, reproducing kernels, linear combinations of decision trees, neural networks, and others. In (A) the parameter spaces are meant, but not required, to contain the functions $z \mapsto \tilde b_0(Z) := b_0(z) + p(z) s_0(z)$ and $z \mapsto s_0(Z)$. In (B) the parameter spaces are meant, but not required, to contain the functions $z \mapsto \bar b_0(Z) := b_0(Z) + (1-p(z)) s_0(z)$ and $z \mapsto s_0(Z)$.

Both (A) and (B) improve over the standard predictive learners that predict $Y$ using the best approximation to $ {\mathrm{E}}[Y \mid D, Z]$ in a given class, but not necessarily the best approximation to the CATE $s_0(Z)$ itself, and may be of independent interest. Moreover, the loss functions in (A) and (B) are also helpful for validation purposes, and choosing the best or aggregating classes of ML methods for targeting the CATE function.

The proposal (B) generalizes and refines the strategy of athey:trees of predicting $YH$ using (a tree form of) $S$ by introducing denoising by $B$. semenova2020 developed a related HT strategy that applies series/sieve learners to the denoised HT-transformed outcome, but it explicitly relies on consistent estimation of the regression function, unlike our approach. We further discuss connections of proposal (A) to the unweighted residual learners of semenova:panel, nie:20, and others below.\footnote{Both of our proposals appeared to be new around the first circulation of this paper as ArXiv:1712.04802. See links to the recent literature below, which proposed related, but different ideas. Relative to the initial version of the paper, the current version contributes with several formal learning properties of (A) and (B). We are grateful to the referees for suggesting that we develop these properties.}

Theorem (ref) shows that the minimizers of the two loss functions provide the best approximation in the mean-square sense to the actual CATE function $s_0(Z)$ in the class $\mathcal{S}$. This property occurs even though we do not observe $s_0(Z)$, and such performance is usually qualified as “oracle." A sufficient condition for the existence of minimizers is that the set $\mathcal{S}$ be convex and closed in the $L^2(P)$ norm.

We illustrate the benefits of using the causal learning objectives (A) and (B) in Figure (ref). In the two panels, we compare the CATE learners derived from the standard predictive Random Forest and Neural Network with the Causal Learners from Random Forests and Neural Network that solve the objective functions (ref) and (ref). We find that the causal learners are better at approximating the CATE function, thereby providing better proxies for CATE. The improvements in the RMSE of approximating the CATE provided by the causal learners range from $21\%$ to $44\%$.

Likewise, the left panel of Figure (ref) shows that we can improve the standard predictive OLS by the Causal OLS that solves the objective functions (A) and (B). Here we estimate linear models in $Z$ for the baseline function and CATE. The improvement in the RMSE of approximating the CATE provided by causal OLS is about $35\%$. This finding might be of interest to researchers using OLS in empirical work.

Finally, the right panel of Figure (ref) shows that one can improve the Causal Forest by a causal boosting step that solves the objective functions (A) and (B) by looking for a shallow forest deviation away from the cross-fitted Causal Forest proxy. The improvements in the RMSE of approximating the true CATE provided by this step are $54-63\%$. The explanation for this improvement is that the Causal Forest, while explicitly targeting CATE, actually solves a different objective than (A) or (B).\footnote{In our understanding it performs a residual learning approach, but not the weighted residual learning approach, which makes the method under-perform relative to the Forest Causal Learner based upon (A) or (B).} We also verified that this improvement only applies when the propensity score is not constant, like in this example.

Learning Guarantees

The learning guarantees of the causal learners (ref) and (ref) follow from the state-of-art statistical learning theory \citep*{rakhlin:offset}, in particular through the use of the expected off-set Rademacher complexity (ORC).

rakhlin:offset define the expected ORC of the function class $\mathcal{H}$ as: $$ \mathcal{R}^o(A, \mathcal{H},c) := {\mathrm{E}} \sup_{h \in \mathcal{H}} \frac{1}{|A|} \sum_{i \in A} \left[ e_i h(Z_i) - c h(Z_i)^2 \right], $$ where $\{e_i\}$ are i.i.d. Rademacher variables taking values $-1$ and $1$ with probability $1/2$, that are generated independently of the data $\{Z_i\}_{i \in A}$, and $c>0$ is a positive constant. ORC is a statistical measure of complexity that captures the ability of the functional class to fit i.i.d. Rademacher noise. The more complex the class is, the higher the ORC. As shown in rakhlin:offset, the expected ORC naturally scales like $ d/|A|$, where $d$ is the effective dimension of the function class and $|A|$ is the sample size. For example, for linear classes, $d$ is the actual dimension of the linear class; and for VC classes, $d$ scales like the VC dimension. rakhlin:offset show that ORC is the sharpest characterization of complexity: in particular, the ORC is upper bounded by the standard critical radii of function classes defined in terms of local Rademacher complexity wainwright:book or in terms of the standard uniform covering entropy Dudley2000. This makes the ORC bounds readily available for all function classes used in modern ML.

The following result establishes formal learning guarantees for causal learners in randomized experiments under general settings that do not assume we can learn $s_0$ consistently.

The result shows that if the ORC of the functional parameter spaces is small, the excess risk of this estimator relative to the oracle approximation to the CATE is small. Note that the lower bound also bounds the distance of $S$ to the oracle (best-in-class) $s_\bullet$ predictor of the CATE. Therefore, the bounds on the excess risk and distance readily follow from the existing characterization of the ORC. For example, if $\mathcal{H}$ has VC type covering entropy with VC index $d$, then the ORC is of order $d/|A|$.

Since the “base" functions $B$'s play only a noise-reducing rule, we can always select $\mathcal{B}$ to be no more complex than $\mathcal{S}$. For example, we can use $\mathcal{B} \subseteq \mathcal{S}$ or pre-train $B$ using a separate auxiliary sample, in which case $\mathcal{B}$ is a singleton.\footnote{This also applies to cross-fitting, which is a better form of sample-splitting for practice.} In either case, learning the technical baseline function does not affect the rate of learning the oracle prediction $s_\bullet$; and the ORC is determined solely by the complexity of $\mathcal{S}$.

Further, we can use losses (ref) and (ref) for choosing the best ML method. We provide detailed discussion in Appendix (ref) of the OA.

Application: Where are nudges for immunization the most effective?

We apply our methods to an RCT in India that was conducted to improve immunization and provide detailed implementation algorithms. We first describe the setting and results, and then provide the implementation details. Our main specification reports median intervals (MI) from causal learners via Boosting as described in Algorithm (ref). We also estimate results using predictive ML methods and report them in Appendix (ref) of the OA as they perform worse than the causal learners. Finally, inferential results are robust to using prediction intervals for the majority values (MP, reported in Appendix (ref) of the OA). For the sample splitting, we allocate 1/3 of the sample to the main sample.\footnote{We find similar results using 1/2 splits. These results are available upon request.}

Setting

Immunization is widely believed to be one of the most cost-effective ways to save children's lives. Much progress has been made in increasing immunization coverage since the 1990s. For example, according to the World Health Organization (WHO), global measles deaths have decreased by 73% from 536,000 estimated deaths in 2000 to 142,000 in 2018. In the last few years, however, global vaccination coverage has remained stuck at around 85% (until the COVID-19 epidemics, when they plummeted). In 2018, 19.7 million children under the age of one year did not receive basic vaccines. Around 60% of these children lived in ten countries: Angola, Brazil, the Democratic Republic of the Congo, Ethiopia, India, Indonesia, Nigeria, Pakistan, the Philippines, and Vietnam. The WHO estimates that immunization saves 2-3 million deaths every year and that an additional 1.5 million deaths could be averted every year if global vaccination coverage improves (this is comparable to 689,000 deaths from COVID-19 between January and August 2020).\footnote{See WHO “10 facts on immunization”, \url{https://www.who.int/features/factfiles/immunization/facts/en/index1.html}}

While most of the early efforts have been devoted to building an immunization infrastructure and ensuring that immunization is available close to people's homes, there is a growing recognition that it is important to also address the demand for immunization. Part of the low demand reflects deep-seated mistrust, but in many cases, parents seem to be perfectly willing to immunize their children. For example, in our data for Haryana, India, among the sample's older siblings who should all have completed their immunization course, 99% had received polio drops, and about 90% had an immunization card. 90% of the parents claimed to believe immunization is beneficial, and 3% claimed to believe it is harmful. However, only 37% of the older children had completed the course and received the measles vaccine, according to their parents (which is likely to be an overestimate), and only 19.4% had done so before the fifteen month of life, when it is supposed to be done between the 10th and the 12th month. It seems that parents lose steam over the course of the immunization sequence, and nudges could be helpful to boost demand. Indeed, recent literature cited in the introduction suggests that “nudges,” such as small incentives, leveraging the social network, SMS, etc., may have a large effect on those services.

In 2017, Esther Duflo, one of the authors of this paper, led a team that conducted a large-scale experiment with the government of Haryana in North India to test various strategies to increase the takeup of immunization services. The government health system rolled out an e-health platform designed by a research team and programmed by an MIT group (SANA health), in which nurses collected data on which child was given which shot at each immunization camp. The platform was implemented in over 2,000 villages in seven districts and provides excellent administrative data on immunization coverage.\footnote{banerjee2019improving discuss validation data from random checks conducted by independent surveyors.} From the individual data, we constructed the monthly sum of the number of children eligible for the program (i.e., age 12 months or younger at their first vaccines) who received each particular immunization at a program location. These children were aged between 0 and 15 months. This paper focuses on the number of children who received the measles shot, as it is the last vaccine in the sequence and thus a reliable marker for full immunization.

Before the launch of the interventions, survey data were collected in 912 of those villages using a sample of 15 households with children aged 1-3 per village. The baseline data covers demographic and socio-economic variables and the immunization history of these children, who were too old to be included in the intervention. In these 912 villages, three different interventions (and their variants) were cross-randomized at the village level:

For each of these interventions, there were several possible variants: incentives were either low or high and either flat or increasing with each shot; the immunization ambassadors were either randomly selected or chosen to be information hubs, using the “gossip” methodology developed by banerjee2019leveraging, a trusted person, or both; and reminders were sent to either 33% or 66% of the people concerned. Moreover, each intervention was cross-cut, generating 75 possible treatment combinations.

banerjee2019leveraging developed and implemented a two-step methodology to identify the most cost-effective and the most effective policy to increase the number of children completing the full course of immunization at the village level and estimate its effects (correcting for bias due to the fact that the policy is found to be the best). First, they used a specific version of LASSO to determine which policies are irrelevant and which policy variants can be pooled together. Second, they obtained consistent estimates of this restricted set of pooled policies using post-LASSO CHS:PnP. They found that the most cost-effective policy (and the only one to reduce the cost of each immunization compared to the status quo) is to combine information hub ambassadors (trusted or not) and SMS reminders. But the policy that increases immunization the most is the combination of information-hub ambassador, the presence of reminders, and increasing incentives (regardless of levels). This is also the most expensive package, so the government was interested in prioritizing villages: where should they scale up the full package? This is an excellent application of this methodology because there was no strong prior.

Results

We compare 25 treated villages where this particular policy bundle was implemented with 78 control villages that received neither sloped incentives and social network intervention nor reminder. Our data constitute an approximately balanced monthly panel of the 103 treated and control villages for 12 months (the duration of the intervention). The outcome variable, $Y$, is the number of children 15 months or younger in a given month in a given village who receive the measles shot. The treatment variable, $D$, is an indicator of the household being in a village that receives the policy. The covariates, $Z$, include 36 baseline village-level characteristics such as religion, caste, financial status, marriage and family status, education, and baseline immunization. The propensity score is constant.

Table (ref) shows sample averages in the control and treated groups for some of the variables used in the analysis weighted by village population, as the rest of the analysis. Treatment and control villages have similar baseline characteristics (in particular, the immunization status of the older cohort was similar). The combined treatment was very effective on average. During the course of the intervention, on average, 7.30 children per month aged 15 months or less got the measles shot that completes the immunization sequence in control villages, and 10.08 did so in treatment villages. This is a raw difference of 2.77 or 38% of the control mean. Note that while these effects are not insignificant, we are far from reaching full immunization: The baseline survey suggests that about $38\%$ of children aged 1-3 had received the measles shot at baseline, and $19.4\%$ had received it before they turned 15 months. These estimates imply that the fraction getting their measles shot before 15 months would only go up to $26.7\% (19.4+0.38*19.4).$

The implementation details for the heterogeneity analysis follow Algorithm (ref) below, with three characteristics due to the design: we weight village-level estimations by village population, include district--time fixed effects, and cluster standard errors at the village level. Table (ref) compares the four ML methods for producing proxy predictors $S(Z_i)$ using the criteria in (ref) and (ref). We find that Elastic Net and Neural Network outperform the other methods, with Elastic Net beating Neural Network by a smaller margin than the other methods. Accordingly, we shall focus on these two methods for the rest of the analysis.

Table (ref) presents the results of the BLP of CATE on the ML proxies. We report estimates of the coefficients $\beta_1$ and $\beta_2$, which correspond to the ATE and heterogeneity loading (HET) parameters in the BLP. The ATE estimates in Columns 1 and 3 indicate that the package treatment increases the number of immunized children by $2.81$ based on elastic net estimates and by $2.44$ based on neural network estimates. Reassuringly, these estimates are on either side of the raw difference in means ($2.77$). Focusing on the HET estimates, we find strong heterogeneity in treatment effects, as indicated by the statistically significant estimates. Moreover, the estimates are close to 1, suggesting that the ML proxies are good predictors of the CATE.

Next, we estimate the GATES by quintiles of the ML proxies. Figure (ref) presents the estimated GATES coefficients $\gamma_{1}-\gamma_{5}$ along with joint confidence bands and the ATE estimates. In Table (ref) we present the result from the hypothesis test that the difference of the ATE for the most and least affected groups is statistically significant. We find that this difference is $21.60$ and $18.13$ based on elastic net and neural network methods, respectively, and is statistically significant. Given that the ATE estimates in the whole population are about $2.5$, these results suggest a large and potentially policy-relevant heterogeneity.

The analysis so far reveals very large heterogeneity, with two striking results. First, the results are very large for the most affected villages. In these villages, an average of $13.23$ extra children eligible for baseline incentives get the measles vaccines every month (starting from a mean of $2.19$ in the elastic net estimation). Second, the impact is negative and significant in the least affected villages (an average decline of $8.00$ immunization per month, starting from $12.68$ in the elastic net estimation). It looks like in some contexts, the combined package of small incentives, reminders, and persuasion by members of the social network put people off immunization.

Given these large differences, it is important to determine whether this heterogeneity seems to be associated with pre-existing characteristics. To answer this question, we ask what variables are associated with the heterogeneity detected in BLP and GATES via CLAN. Table (ref) reports the CLAN estimates for a selected set of covariates and Tables (ref)--(ref) in Appendix (ref) of the OA for the rest of covariates. Regardless of the method used, the estimated differences in the means of most and least affected groups for the number of vaccines to child since birth, number of polio drops to child, the fraction of children receiving measles vaccines by 15 months of age, and fraction of children receiving measles vaccine at credible locations, are negative and statistically significant. Those are various measures of pretreatment immunization levels, all survey-based, that have nothing to do with our measure of impact. These results suggest that the villages with low levels of pretreatment immunization are the most affected by the incentives. These are, in fact, the only variables that consistently pop up from the CLAN. Thus, in this instance, the policy preferred ex-ante by the government (since it is equality-enhancing) also happens to be the most effective.

While the heterogeneity associated with the baseline immunization rates cannot be causally interpreted (it could always be proxying for other things), it still sheds light on the negative effect we find for the least affected group. Note that this effect is not mechanical. Even in the least affected villages, there was a good number of children who did not receive the measles shot, and they were not close to reaching full immunization, where they could not have experienced an increase. It may be that it would have been difficult to vaccinate $13.23$ extra children every month, but there was scope to experience an increase in immunization. We had no prior on that the effect would be larger in the villages with the lowest immunization rate. On the contrary, immunization rates could have been low precisely because parents were more doubtful about immunization. For example, immunization is particularly low in Muslim-majority villages, which is believed to reflect their lack of trust in the health system. There were, therefore, reasons to be genuinely uncertain about where immunization would have had the largest effect.

One possible interpretation of the negative impact in some villages is that villagers were intrinsically motivated to get immunized. The nudging with small incentives and mild social pressure may have backfired, crowding out intrinsic motivation without providing a strong enough extrinsic motivation to act as in gneezy2000fine. A point estimate of $13.23$ extra immunization per month in the most affected group might seem high: a multiplication by $6.0$ of the baseline level (based on the elastic net specification). This increase in immunizations is not inconsistent with the literature: in a set of villages with a very low immunization rate in Rajasthan, banerjee2010improving find that small incentives increase immunization from $18\%$ to $39\%$ (relative to a treatment that just improves infrastructures, and $6\%$ relative to the control group) in a low immunization region (in the entire sample, not in the places where it is most effective), which was also a very large increase. Given the restrictions imposed on the data set (only children 1 year or less at their first immunization were included), the data cover children who were 15 months or younger when getting the measles shot. Among the older cohort in the most affected group, only $12.7\%$ of children were vaccinated before 15 months. Taking this as a benchmark for the control group, the estimate would still imply that only $64\%$ of the treatment group was immunized before 15 months: a big improvement but not implausible.

Our last exercise is to compute the cost-effectiveness of the program in various groups. To do so, we compute in each village the average number of immunizations delivered per dollar spent in a month in each group. The dollar spent is the fixed cost to run an immunization program per month (nurse salaries, administrative overheads, etc.) plus the marginal cost of each vaccine multiplied by the number of vaccines administered (incentives distributed to local health workers, vaccines doses, syringes, etc.) in both treatment and control villages, plus the extra cost of running each particular treatment (the cost of the tablets used for recording in all the treatment villages, the cost of contacting and enrolling the ambassadors, and the cost of the incentives). We then estimate the cost-effectiveness in each GATES group as ${\mathrm{E}}[X(1) - X(0) \mid G_k]$, where $X$ is the immunizations per dollar. ${\mathrm{E}}[X(1) - X(0) \mid G_k] = {\mathrm{E}}[X \mid D = 1 , G_k] - {\mathrm{E}}[X \mid D = 0, G_k]$ by the randomization assumption, and we can estimate each of ${\mathrm{E}}[X \mid D = 1 , G_k]$ and ${\mathrm{E}}[X \mid D = 0 , G_k]$ analogously to CLAN, that is by taking sample averages within treatment groups for each sample split and aggregating over sample splits.

The results are presented in Table (ref). They highlight the crucial importance of treatment effect heterogeneity for policy decisions in this context. Overall, as shown in banerjee2019leveraging the treatment is not cost-effective compared to the control (the immunization per dollar spent goes down). This analysis reveals that this is driven (not surprisingly) by negative impacts on cost-effectiveness in the groups where it is least effective. However, in the fourth and fifth quintile of cost-effectiveness, we cannot reject that the immunization per dollar spent is the same in the control group and in the treatment group, despite the added marginal cost of the incentives and the vaccines: this is because the fixed cost of running the program is now spread over a larger number of immunizations.

We performed the following additional analysis. First, since the main result in banerjee2019leveraging is that the most cost-effective option on average is the combination of SMS plus Information hubs, an alternative policy question may therefore be whether there are places where it may be more cost-effective to add the incentives to this cheaper treatment. We replicated the heterogeneity analysis comparing these two treatments and looked at the cost-effectiveness of GATES in these two options. There is also considerable heterogeneity in this comparison (see Figure (ref) of the OA). The results for cost-effectiveness are shown in Table (ref) of the OA. There, again, we find that in the two quintiles where adding incentives is most effective, it would be cost-effective, even compared to an alternative status quo of just having SMS and information hubs.

Second, we reproduce the analysis with the first vaccine in the series to be given at birth (Penta 1). The results are very similar to the results for the measles vaccine, with regions with the lowest immunization experiencing the largest effects. Third, still for the penta 1 vaccination, we apply the method to the comparison of villages with or without increasing incentives for vaccination. Once again, we find significant heterogeneity with the most affected regions being those with low immunization to start with (these results are available upon request).

Implementation Details

We describe two general algorithms and provide some specific implementation details for the empirical example.

We implemented our causal learners via a boosting approach that looks for relatively simple deviations from the initial predictive learner to improve CATE predictions. The reason is that this approach performed better in our simulation experiments than directly solving the objective functions (A) and (B) over large parameter spaces. We also observed that this boosting implementation performed better in the empirical example. The performance improvement occurs because the objective functions (A) or (B) tend to be much noisier than the objective functions in predictive learning and harder to tune. The following algorithm summarizes the implementation.

Ideally, step 1 of the algorithm should make use of cross-fitting, but we present a simplified version for clarity.

Conclusion and Extensions

We propose to focus inference on key features of heterogeneous effects in randomized experiments, and develop the corresponding methods. These key features include best linear predictors of the effects and average effects sorted by groups, as well as average characteristics of most and least affected units. Our approach is valid in high dimensional settings, where the effects are estimated by machine learning methods. The main advantage of our approach is its agnostic nature; it avoids making strong assumptions. Estimation and inference relies on data splitting, where the latter allows us to avoid overfitting and all kinds of non-regularities. Our inference aggregates the results across many splits, reducing the replication risks, and could be of independent interest. An empirical application illustrates the practical use of the approach.

Our hope is that applied researchers use the method to discover whether there is heterogeneity in their data in a disciplined way. A researcher might be concerned about the application of our method due to the possible power loss induced by sample splitting. This power loss is the price to pay when the researcher is not certain or willing to fully specify the form of the heterogeneity prior to conducting the experiment. Thus, if the researcher has a well-defined pre-analysis plan that spells out a small number of heterogeneity groups in advance, then there is no need of splitting the sample.\footnote{More generally, the plan needs to specify a parametric form for the heterogeneity as a low dimensional function of pre-specified covariates CFL2014. In this case, ML tools can still be used to efficiently estimate the CATEs in the presence of control variables but are not required to detect heterogeneity bcfh17, DML.} However, this situation is not common. In general, the researchers might not be able to fully specify the form of the heterogeneity due to lack of information, economic theory, or willingness to take a stand at the early stages of the analysis. They might also face data limitations that preclude the availability of the desired covariates. Here we recommend the use of our method to avoid overfitting and p-hacking, and impose discipline to the heterogeneity analysis at the cost of some power loss due to sample splitting.\footnote{ It is not clear whether this loss is real though, as we are not aware of any alternative method that avoids sample splitting and that works at the same level of agnosticism as ours. In a previous version of the paper we provided a numerical example using a simple parametric model where standard methods without sample splitting are available. We find that the extent of the power loss for not using the parametric form of the heterogeneity roughly corresponds to reducing the sample size by half in a test for the presence of heterogeneity, although the exact comparison depends on features of the data generating process.} If discovering and exploiting heterogeneity in treatment effect is a key goal of the research, the researcher should indeed plan for larger sample sizes (relative to just testing whether the treatment has an effect), but the required sample size remains within the realm of what is feasible in the field. In many applications we are aware of, there was apparent heterogeneity according to some covariates of interest, but the disciplined ML heterogeneity exercise found no systematic difference. This could be because this heterogeneity was a fluke, or because the method does not have the power to detect it in a small sample. In any case, what this experience suggests is that one should not rely on ex-post heterogenous effects in such cases.

The application to immunization in India is of substantive interest. Our findings suggest that a combination of small incentives, relay by information hub, and SMS reminders can have very large effect on vaccine take up in some villages where immunization was low at baseline, and even be as cost-effective than the status quo, but can also backfire in other places. This suggests that these types of strategy need to be piloted in the relevant context before being rolled out, and that heterogeneity needs to be taken into account.

Our inference approach generalizes to any problem of the following sort, studied in semenova2020 using more conventional inference approaches. Suppose we can construct an unbiased signal $\tilde Y$ such that $${\mathrm{E}} [\tilde Y \mid Z] = s_0(Z),$$ where $s_0(Z)$ is now a generic target function. Let $S(Z)$ denote an ML proxy for $s_0(Z)$. This setting generalizes the CATE inference problem of this paper.\footnote{The unbiased signals arise from multiplying an outcome with a Riesz representer for the effect of interest. In CATE setting, the representer is the Horwitz-Thomposon transform.} We can then apply our inferential framework to this setting. The setting covers a variety of causal effects of interest, for example, inference on conditional average causal derivatives, conditional average effects from transporting and distributional shifts of covariates, as we further explain in Appendix (ref) of the OA.