GENERIC MACHINE LEARNING INFERENCE ON HETEROGENOUS TREATMENT
EFFECTS IN RANDOMIZED EXPERIMENTS

arXiv:1712.04802v3 [stat.ML] 23 Feb 2018

V. CHERNOZHUKOV, M. DEMIRER, E. DUFLO, I. FERNANDEZ-VAL

Abstract. We propose strategies to estimate and make inference on key features of heterogeneous effects in randomized experiments. These key features include best linear predictors of the effects using machine learning proxies, average
effects sorted by impact groups, and average characteristics of most and least impacted units. The approach is valid in high
dimensional settings, where the effects are proxied by machine learning methods. We post-process these proxies into
the estimates of the key features. Our approach is generic, it can be used in conjunction with penalized methods, deep
and shallow neural networks, canonical and new random forests, boosted trees, and ensemble methods. Our approach
is agnostic and does not make unrealistic or hard-to-check assumptions; we don’t require conditions for consistency of
the ML methods. Estimation and inference relies on repeated data splitting to avoid overfitting and achieve validity.
For inference, we take medians of p-values and medians of confidence intervals, resulting from many different data
splits, and then adjust their nominal level to guarantee uniform validity. This variational inference method is shown
to be uniformly valid and quantifies the uncertainty coming from both parameter estimation and data splitting. The
inference method could be of substantial independent interest in many machine learning applications. An empirical
application to the impact of micro-credit on economic development illustrates the use of the approach in randomized
experiments. An additional application to the impact of the gender discrimination on wages illustrates the potential use
of the approach in observational studies, where machine learning methods can be used to condition flexibly on very
high-dimensional controls.
Key words: Agnostic Inference, Machine Learning, Confidence Intervals, Causal Effects, Variational P-values and
Confidence Intervals, Uniformly Valid Inference, Quantification of Uncertainty, Sample Splitting, Multiple Splitting,
Assumption-Freeness

1. Introduction
Randomized experiments play an important role in the evaluation of social and economic programs and medical treatments (e.g., [35, 27]). Researchers and policy makers are often interested
in features of the impact of the treatment that go beyond the simple average treatment effects. In
particular, very often, they want to know whether treatment effect depends on basic covariates,
such as gender, age, etc. It is essential to assess if the impact of the program would generalize to
a different population with different characteristics, and for economists, to better understand the
driving mechanism behind a particular program effects.
One issue with reporting treatment effects split by subgroups, however, is that there are often a large number of potential sample splits: choosing subgroups ex-post opens the possibility
Date: December, 2017.
We thank Susan Athey, Moshe Buchinsky, Denis Chetverikov, Siyi Luo, Max Kasy, Susan Murphy, Whitney Newey,
and seminar participants at UCLA and AEA for valuable comments.
1

2

V. CHERNOZHUKOV, M. DEMIRER, E. DUFLO, I. FERNANDEZ-VAL

of overfitting. To solve this problem, medical journals and the FDA require pre-registering subsample of interest in medical trials in advance. In economics, this approach has gained some traction, with the adoption of pre-analysis plan (which can be filed in the AEA registry for randomized
experiments). Restricting heterogeneity analysis to pre-registered subgroups, however, 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 co-variates.
To do this in a disciplined fashion and avoid the risk of overfitting, scholars have recently proposed using machine learning tools (see e.g. [6] and below for a review). Indeed, ML tools seem
to be ideal to explore heterogeneity of treatment effect, when researchers have access to a potentially large array of baseline variables to form subgroup, and little 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.1
This paper develops a generic approach to use any of the tools of machine learning 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. In fact, in high
dimensional settings, absent strong assumptions, generic ML tools may not even produce consistent estimates of the conditional average treatment effect (or CATE) (the difference in the expected
potential outcomes between treated and control groups conditional on the covariates Z).
Previous attempts to solve this problem focus on specific tools and situations where those assumptions are satisfied. Our approach to resolve the fundamental impossibilities in non-parametric
inference is different. Motivated by [29], instead of attempting to get consistent estimate and uniformly valid inference on CATE itself, we focus on providing valid estimation and inference for
features of CATE. We start by building a ML proxy predictor of CATE, and then develop valid inference for features of the CATE based on this proxy predictor. In particular, we develop valid
inference for the following three objects, which are likely to be of interest to applied researchers
and policy makers: First, the Best Linear Predictor (BLP) of the CATE based on the ML proxy predictor, second, the Sorted Group Average Treatment Effects (GATES), average treatment effect by
heterogeneity groups induced by the ML proxy predictor S(Z); third, the Classification Analysis
(CLAN): 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 is the treatment effect for different bins. And finally
we can describe which of the covariates is correlated with this heterogeneity.
1

In the last few months alone, several new empirical papers in economics used ML methods to estimate heterogenous

effects. E.g. [44], which shows that villagers outperform the machine learning tools when they predict heterogeneity in
returns to capital. [25] predict who benefits the most from a summer internship projects. The methodological papers
reviewed later also contain a number of empirical applications.

GENERIC ML FOR FEATURES OF HETEROGENOUS TREATMENT EFFECTS

3

Thus, in the trade off between more restrictive assumptions and a more ambitious estimation,
we chose a different point than previous papers: focus 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 that for which we have developed inference appear to us at this point to be
the most relevant, but in the future, one could easily develop methods for other objects of interest.
The 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; see [46]. Let Z be a vector of of covariates that characterize the observational units. The main causal functions are the baseline conditional average
(BCA):
b0 (Z) := E[Y (0) | Z],

(1.1)

and the conditional average treatment effect (CATE):
s0 (Z) := E[Y (1) | Z] − E[Y (0) | Z].

(1.2)

Suppose the binary treatment variable D is randomly assigned conditional on Z, with probability of assignment depending only on a subvector of stratifying variables Z1 ⊆ Z, namely
D ⊥⊥ (Y (1), Y (0)) | Z,

(1.3)

and the propensity score is known and is given by
p(Z) := P[D = 1 | Z] = P[D = 1 | Z1 ],

(1.4)

which we assume is bounded away from zero or one:
p(Z) ∈ [p0 , p1 ] ⊂ (0, 1).

(1.5)

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:
Y = b0 (Z) + Ds0 (Z) + U, E[U | Z, D] = 0,
that is,
b0 (Z) = E[Y | D = 0, Z],

(1.6)

s0 (Z) = E[Y | D = 1, Z] − E[Y | D = 0, Z].

(1.7)

and

We observe Data = (Yi , Zi , Di )N
i=1 , consisting of i.i.d. copies of the random vector (Y, Z, D)
having probability law P . Expectation with respect to P is denoted by E = EP . The probability
law of the entire data is denoted by P = PP and the corresponding expectation is denoted by
E = EP .

4

V. CHERNOZHUKOV, M. DEMIRER, E. DUFLO, I. FERNANDEZ-VAL

Properties of Machine Learning Estimators of s0 (Z) Motivating the Agnostic Approach. Machine learning (ML) is a name attached to a variety of new, constantly evolving statistical learning
methods: Random Forest, Boosted Trees, Neural Networks, Penalized Regression, Ensembles, and
Hybrids (see, e.g., [50] for a recent review, and [28] for a prominent textbook treatment). In modern
high-dimensional settings, machine learning methods effectively explore the various forms of nonlinear structured sparsity to yield “good” approximations to s0 (z) whenever such assumptions are
valid, based on equations (1.6) and (1.7). As a result these methods often work much better than
classical methods in high-dimensional settings, and have found widespread uses in industrial and
academic applications.
Motivated by the practical predictive success of the methods, it is really tempting to apply the
ML methods directly to try to learn the CATE function z 7→ s0 (z) (by learning the two regression
functions for treated and untreated and taking the difference). However, it is hard, if not impossible, to obtain uniformly valid inference on z 7→ s0 (z) using generic ML methods, under credible
assumptions and practical tuning parameter choices. There are several fundamental reasons as
well as huge gaps between theory and practice that are responsible for this.
One fundamental reason is that the ML methods might not even produce consistent estimators
of z 7→ s0 (z) in high dimensional settings. For example, if z has dimension d and the target function
z 7→ s0 (z) is assumed to have p continuous and bounded derivatives, then the worst case (minimax)
lower bound on the rate of learning this function cannot be better than N −p/(2p+d) as N → ∞, as
shown by Stone [47]. Hence if p is fixed and d is also small, slowly increasing with N , such as
d > log N , then there exists no consistent estimator of z 7→ s0 (z) generally.
Hence, ML estimators cannot be regarded as consistent, unless further very strong assumptions
are made. Examples of such assumptions include structured forms of linear and non-linear sparsity and super-smoothness. While these (sometime believable and yet untestable) assumptions
make consistent adaptive estimation possible (e.g.,[17]), inference remains a more difficult problem, as adaptive confidence sets do not exist even for low-dimensional nonparametric problems
([41, 29]). Indeed, adaptive estimators (including modern ML methods) have biases of comparable
or dominating order as compared to sampling error. Further assumptions such as ”self-similarity”
are needed to bound the biases and expand the confidence bands by the size of bias (see [30, 22])
to produce partly adaptive confidence bands. For more traditional statistical methods there are
constructions in this vein that make use of either undersmoothing or bias-bounding arguments
([30, 22]). These methods, however, are not yet available for ML methods in high dimensions (see,
however, [31] for a promising approach called ”targeted undersmoothing” in sparse linear models).
Suppose we did decide to be optimistic (or panglossian) and imposed the strong assumptions,
that made the theoretical versions of the ML methods provide us with high-quality consistent estimators of z 7→ s0 (z) and valid confidence bands based on them. This would still not give as a practical construction we’d want for our applications. The reason is that there is a huge gap between

GENERIC ML FOR FEATURES OF HETEROGENOUS TREATMENT EFFECTS

5

theoretical versions of the ML methods appearing in various theoretical papers and the practical
versions (with the actual, data-driven tuning parameters) coded up in statistical computing packages used by practitioners2.3 The use of ML, for example, involves many tuning parameters with
practical rules for choosing them, while theoretical works provide little guidance or backing for
such practical rules; see e.g., the influential book [28] for many examples of such rules.4
In this paper we shall take an agnostic view and we shall not rely on any structured assumptions,
which might be difficult to verify or believe in practice, and we don’t impose conditions that make
the ML estimators consistent. We simply treat ML as providing proxy predictors for the objects of
interest.
Our Agnostic Approach. Our approach will be radically different – we propose strategies for estimation and inference on
key features of s0 (Z) rather than s0 (Z) itself.
Because of this difference in focus we can avoid making strong assumptions about the properties
of the ML estimators.
Let (M, A) denote a random partition of the set of indices {1, . . . , N }. The strategies that we
consider rely on random splitting of
Data = (Yi , Di , Zi )N
i=1
into a main sample, denoted by DataM = (Yi , Di , Zi )i∈M , and an auxiliary sample, denoted by
DataA = (Yi , Di , Zi )i∈A . We will sometimes refer to these samples as M and A. We assume that
the main and auxiliary samples are approximately equal in size, though this is not required theoretically.
From the auxiliary sample A, we obtain ML estimators of the baseline and treatment effects,
which we call the proxy predictors,
z 7→ B(z) = B(z; DataA ) and z 7→ S(z) = S(z; DataA ).
These are possibly biased and noisy predictors of b0 (z) and s0 (z), and in principle, we do not even
require that they are consistent for b0 (z) and s0 (z). We simply treat these estimates as proxies,
which we post-process to estimate and make inference on the features of the CATE z 7→ s0 (z). We
condition on the auxiliary sample DataA , so we consider these maps as frozen, when working with
the main sample.
2There are cases where such gap does not exist, e.g., see [13, 15] for the lasso.
3For example, even the wide use of K-fold cross-validation in high-dimensional settings for machine learning remains

theoretically unjustified (there do exist, however, related subsample-based methods [52, 40] that do achieve near-oracle
performance for tuning selection).
4Unfortunately, theoretical work often only provides existence results: there exist theoretical ranges of the tuning
parameters that make the simple versions of the methods work for predictive purposes (under very strong assumptions),
leaving no satisfactory guide to practice.

6

V. CHERNOZHUKOV, M. DEMIRER, E. DUFLO, I. FERNANDEZ-VAL

Using the main sample and the proxies, we shall target and develop valid inference about key
features of s0 (Z) rather than s0 (Z), which include
(1) Best Linear Predictor (BLP) of the CATE s0 (Z) based on the ML proxy predictor S(Z);
(2) Sorted Group Average Treatment Effects (GATES): average of s0 (Z) (ATE) by heterogeneity groups induced by the ML proxy predictor S(Z);
(3) Classification Analysis (CLAN): average characteristics of the most and least affected units
defined in terms of the ML proxy predictor S(Z).
Our approach is generic with respect to the ML method being used, and is agnostic about its formal
properties.
Our point estimation will make use of many splits of the data into main and auxiliary samples
to produce robust estimates. Our estimation and inference will systematically account for two
sources of uncertainty:
(I) Estimation uncertainty conditional on the auxiliary sample.
(II) Splitting uncertainty induced by random partitioning of the data into the main and auxiliary samples.
Because we account for the second source, we call the resulting collection of methods as variational
estimation and inference methods (VEINs). For point estimates we report the median of the estimated key features over different random splits of the data. For the confidence intervals we take
the medians of many random conditional confidence sets and we adjust their nominal confidence
level to reflect the splitting uncertainty. We construct p-values by taking medians of many random
conditional p-values and adjust the nominal levels to reflect the splitting uncertainty. Note that
considering many different splits and accounting for variability caused by splitting is very important. Indeed, with a single splitting practice, empiricists may unintentionally look for a ”good” data
split, which supports their prior beliefs about the likely results, thereby invalidating inference.5
Relationship to the Literature. We focus the review strictly on the literatures about estimation
and inference on heterogeneous effects and inference using sample splitting.
We first mention work that uses linear and semiparametric regression methods. A semiparametric inference method for characterizing heterogeneity, called the sorted effects method, was
given in [20]. This approach does provide a full set of inference tools, including simultaneous
5This problem is “solved” by fixing the Monte-Carlo seed and the entire data analysis algorithm before the empirical

study. Even if such a huge commitment is really made and followed, there is a considerable risk that the resulting
data-split may be non-typical. Our approach allows one to avoid taking this risk.

GENERIC ML FOR FEATURES OF HETEROGENOUS TREATMENT EFFECTS

7

bands for percentiles of the CATE, but is strictly limited to the traditional semiparametric estimators for the regression and causal functions. [31] proposed a sparsity based method called ”targeted undersmoothing” to perform inference on heterogeneous effects. This approach does allow
for high-dimensional settings, but makes strong assumptions on 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) was first given in [14] using a Z-estimation
framework, where the number of interactions can be very large; see also [26] 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. [53] proposed a post-selection inference framework within the high-dimensional linear sparse models for the heterogeneous effects.
The approach is attractive because it allows for some misspecification of the model.
Next we discuss the use of tree-based and other methods. [34] 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 (e.g., as in [33, 1]) such that that the new outcome becomes an unbiased, noisy version of
CATE. [5] made use of the Horvitz-Thompson transformation of the outcome variable to inform
the process of building causal trees, with the main goal of predicting CATE. They also provide a
valid inference result on average treatment effects for groups defined by the tree leaves, conditional
on the data split in two subsamples: one used to build the tree leaves and the one to estimate the
predicted values given the leaves. Unlike our generic method introduced below, this inference approach is limited to trees and does not account for splitting uncertainty.6 Like our methods, this approach is essentially assumption-free. [49] provided a subsampling-based construction of a causal
random forest, providing valid pointwise inference for CATE (see also the review in [49] on prior
uses of random forests in causal settings) for the case when covariates are very low-dimensional
(and essentially uniformly distributed).7 Unfortunately, this condition rules out the typical highdimensional settings that arise in many empirical problems, including the ones considered in this
paper.
Our approach is radically different from these existing approaches, in that we are changing the
target, and instead of hunting for CATE z 7→ s0 (z), we focus on key features of z 7→ s0 (z). Our
approach is generic, it can be used in conjunction with penalized methods, deep and shallow neural
6In principle, our (generic) inference methods described below applies to the [5]’s problem to provide valid inference

that does account for sample splitting uncertainty. Also, going beyond trees is crucial, since the best ML algorithms are
typically not trees.
7The dimension d is fixed in [49]; the analysis relies on the Stone’s model with smoothness index β = 1, in which no
consistent estimator exists once d > log n. It’d be interesting to establish consistency properties and find valid inferential
procedures for the random forest in high-dimensional (d ∝ n or d  n) approximately sparse cases, with continuous and
categorical covariates, but we are not aware of any studies that cover such settings, which are of central importance to
us.

8

V. CHERNOZHUKOV, M. DEMIRER, E. DUFLO, I. FERNANDEZ-VAL

networks, canonical and new random forests, boosted trees, and ensemble methods. Our approach
is agnostic and does not make unrealistic or hard-to-check assumptions, for example, we don’t
even require conditions for consistency of the ML methods. We simply treat the ML methods as
providing a proxy predictor z 7→ S(z), which we post-process to estimate and make inference
on the key features of the CATE z 7→ s0 (z). Some of our strategies rely on Horvitz-Thompson
transformations of outcome and some do not. The inspiration for our approach draws upon an
observation in [29], namely that some fundamental impossibilities in non-parametric inference
could be avoided if we focus inference on coarser features of the non-parametric functions rather
than the functions themselves.
Our inference approach is also of independent interest, and could be applied to many problems,
where sample splitting is used to produce ML predictions, e.g. [2]. Related references include
[51, 42], where the ideas are related but quite different in details, which we shall explain below.
The premise is the same, however, as in [42, 45] – we should not rely on a single random split
of the data and should adjust inference in some way. Our approach takes the medians of many
conditional confidence intervals as the confidence interval and the median of many conditional
p-values as the p-value, and adjusts their nominal levels to account for the splitting uncertainty.
Our construction of p-values builds upon ideas in [16, 42], though what we propose is radically
simpler, and our confidence intervals appear to be brand new. Of course sample splitting ideas
are classical, going back to [32, 38, 12, 23, 43], though having been mostly underdeveloped and
overlooked for inference, as characterized by [45]. The overlooked status perhaps allowed us to
innovate on the inference side.
2. Main Identification Results and Estimation Strategies
2.1. BLP of CATE. We consider two strategies for identifying and estimating the best linear predictor of s0 (Z) using S(Z):
BLP[s0 (Z) | S(Z)] := arg

min
f (Z)∈Span(1,S(Z))

E[s0 (Z) − f (Z)]2 ,

which, if exists, is defined by projecting s0 (Z) on the linear span of 1 and S(Z) in the space L2 (P ).
BLP of CATE: The First Strategy. Here we shall identify the coefficients of the BLP from the
weighted linear projection:
Y = α0 X1 + β1 (D − p(Z)) + β2 (D − p(Z))(S − ES) + , E[w(Z)X] = 0,
where S := S(Z),
w(Z) = {p(Z)(1 − p(Z))}−1 ,
X1 := X1 (Z),

e.g.,

X := (X1 , X2 )

X1 = [1, B(Z)],

X2 := [D − p(Z), (D − p(Z))(S(Z) − (S − ES)].
Note that the above equation uniquely pins down β1 and β2 under weak assumptions.

(2.1)

GENERIC ML FOR FEATURES OF HETEROGENOUS TREATMENT EFFECTS

9

The interaction (D − p(Z))(S − ES) is orthogonal to D − p(Z) under the weight w(Z) and to all
other regressors that are functions of Z under any Z-dependent weight.8
A consequence is our first main identification result, namely that
β1 + β2 (S(Z) − ES) = BLP[s0 (Z) | S(Z)],
in particular β1 = Es0 (Z) and β2 = Cov(s0 (Z), S(Z))/ Var(S(Z)).
Theorem 2.1 (BLP 1). Consider z 7→ S(z) and z 7→ B(z) as fixed maps. Assume that Y and X have
finite second moments and that EXX 0 is full rank. Then, (β1 , β2 ) defined in (2.1) also solves the best linear
predictor/approximation problem for the target s0 (Z):
(β1 , β2 )0 = arg min E[s0 (Z) − b1 − b2 S(Z)]2 ,
b1 ,b2

in particular β1 = ES0 (Z) and β2 = Cov(s0 (Z), S(Z))/ Var(S(Z)).
The identification result is constructive. We can base the corresponding estimation strategy on
the empirical analog:
Yi = α
b0 X1i + βb1 (Di − p(Zi )) + βb2 (Di − p(Zi ))(Si − EN,M Si ) + b
i ,

i ∈ M,

EN,M [w(Zi )b
i Xi ] = 0,
where EN,M denotes the empirical expectation with respect to the main sample, i.e.
X
EN,M g(Yi , Di , Zi ) := |M |−1
g(Yi , Di , Zi ).
i∈M

The properties of this estimator, conditional on the auxilliary data, are well known and follow as
a special case of Lemma B.1 in the Appendix.
Comment 2.1 (Main Implications of the result). If S(Z) is a perfect proxy for s0 (Z), then β2 = 1.
In general, β2 6= 1, correcting for noise in S(Z). If S(Z) is complete noise, uncorrelated to s0 (Z),
then β2 = 0. Furthermore, if there is no heterogeneity, that is s0 (Z) = s, then β2 = 0. Rejecting
the hypothesis β2 = 0 therefore means that there is both heterogeneity and S(Z) is its relevant
predictor.

Figure 1 provides two examples. The left panel shows a case without heterogeneity in the CATE
where s0 (Z) = 0, whereas there right panel shows a case with strong heterogeneity in the CATE
where s0 (Z) = Z. In both cases we evenly split 1000 observations between the auxiliary and main
8The orthogonalization ideas embedded in this strategy do have classical roots in econometrics (going back to at

least Frisch and Waugh in the 30s), and similar strategies underlie the orthogonal or double machine learning approach
(DML) in [21]. Our paper has different goals than DML, attacking the problem of inference on heterogeneous effects
without rate and even consistency assumptions. The strategy here is more nuanced in that we are making it work under
misspecification or inconsistent learning, which is likely to be true in very high-dimensional problems.

V. CHERNOZHUKOV, M. DEMIRER, E. DUFLO, I. FERNANDEZ-VAL

0.5
−0.5

0.0

S(Z)

0
−4

−1.0

−2

S(Z)

2

1.0

4

10

−1.0

−0.5

0.0

0.5

1.0

−1.0

Z

−0.5

0.0

0.5

1.0

Z

Figure 1. Example. In the left panel we have a homogeneous CATE s0 (Z) = 0; in
the right panel we have heterogeneous CATE s0 (Z) = Z. The proxy predictor S(Z)
is produced by the Random Forest, shown by green line, the true BLP of CATE is
shown by black line, and the estimated BLP of CATE is shown by blue line. The
true and estimated BLP of CATE are more attenuated towards zero than the proxy
predictor.
samples, Z is uniformly distributed in (−1, 1), and the proxy predictor S(Z) is estimated by random
forest in the auxiliary sample following the standard implementation, see e.g. [28]. When there
is no heterogeneity, post-processing the ML estimates helps reducing sampling noise bringing the
estimated BLP close to the true BLP; whereas under strong heterogeneity the signal in the ML
estimates dominates the sampling noise and the post-processing has little effect.
Comment 2.2 (Digression: Naive Strategy that is not Quite Right). It is tempting and “more natural” to estimate
Y = α̃1 + α̃2 B + β̃1 D + β̃2 D(S − ES) + , E[X̃] = 0,
where X̃ = (1, B, D, D(S − ES)). This is a good strategy for predicting the conditional expectation
of Y given Z and D. But, β̃2 6= β2 , and β̃1 + β̃2 (S − ES) is not the best linear predictor of s0 (Z). 
BLP of CATE: The Second Strategy. The second strategy makes use of the Horvitz-Thompson
transformation:
D − p(Z)
.
H = H(D, Z) =
p(Z)(1 − p(Z))
It is well known that the transformed response Y H provides an unbiased signal about CATE:
E[Y H | Z] = s0 (Z)

GENERIC ML FOR FEATURES OF HETEROGENOUS TREATMENT EFFECTS

11

and it follows that
BLP[s0 (Z) | S(Z)] = BLP[Y H | S(Z)].
This simple strategy is completely fine for identification purposes, but can severely underperform
in estimation and inference due to lack of precision. We can repair the deficiencies as follows.
We consider, instead, the linear projection:
Y H = µ0 X1 H + β1 + β2 (S − ES) + ˜,

E˜
X̃ = 0,

(2.2)

where B := B(Z), S := S(Z), and X̃ := (X10 H, X̃20 )0 , X̃2 = (1, (S − ES)0 )0 , where X1 = X1 (Z), e.g.
B(Z) or (B(Z), S(Z), p(Z))0 . The terms X1 are present in order to reduce noise.
We show that, as a complementary main identification result,
β1 + β2 (S − ES) = BLP[s0 (Z) | S(Z)].

Theorem 2.2 (BLP 2). Consider z 7→ S(z) and z 7→ B(z) as fixed maps. Assume that Y has finite second
moments and X̃ = (X1 H, 1, (S − ES)) is such that EX̃ X̃ 0 is finite and full rank. Then, (β1 , β2 ) defined
in (2.2) solves the best linear predictor/approximation problem for the target s0 (Z):
(β1 , β2 )0 = arg min E[s0 (Z) − b1 − b2 S(Z)]2 ,
b1 ,b2

in particular β1 = Es0 (Z) and β2 = Cov(s0 (Z), S(Z))/ Var(S(Z)).
The corresponding estimator is defined through the empirical analog:
Yi Hi = µ
b0 X1i Hi + βb1 + βb2 (Si − EN,M Si ) + b
i ,

EN,M b
i X̃i = 0,

and the properties of this estimator, conditional on the auxiliary data, are well known and given
in Lemma B.1.
2.2. The Sorted Group ATE. The target parameters are
E[s0 (Z) | G],
where G is an indicator of group membership.
Comment 2.3. There are many possibilities for creating groups based upon ML tools applied to
the auxiliary data. For example, one can group or cluster based upon predicted baseline response
(as in the “endogenous stratification” analysis, [2]) or based upon actual predicted treatment effect
S. We focus on the latter approach for defining groups, although our identification and inference
ideas immediately apply to other ways of defining groups, and could be helpful in these contexts.
We build the groups to explain as much variation in s0 (Z) as possible
Gk := {S ∈ Ik },

k = 1, ..., K,

12

V. CHERNOZHUKOV, M. DEMIRER, E. DUFLO, I. FERNANDEZ-VAL

where Ik = [`k−1 , `k ) are non-overlaping intervals that divide the support of S into regions [`k−1 , `k )
with equal or unequal masses:
−∞ = `0 < `1 < . . . < `K = +∞.
The parameters of interest are the Sorted Group Average Treatment Effects (GATES):
E[s0 (Z) | Gk ],

k = 1, . . . , K.

Given the definition of groups, it is natural for us to impose the monotonicity restriction
E[s0 (Z) | G1 ] 6 ... 6 E[s0 (Z) | GK ],
which holds asymptotically if S(Z) is consistent for s0 (Z) and the latter has an absolutely continuous distribution. Under the monotonicity condition, the estimates could be rearranged to obey the
weak monotonicity condition, improving the precision of the estimator. The joint confidence intervals could also be improved by intersecting them with the set of monotone functions. Furthermore,
as before, we can test for homogeneous effects, s0 (Z) = s, by testing whether,
E[s0 (Z) | G1 ] = ... = E[s0 (Z) | GK ].
GATES: The First Strategy. Here we shall recover GATES parameters from the weighted linear
projection equation:
0

Y = α X1 +

K
X

γk · (D − p(Z)) · 1(Gk ) + ν,

E[w(Z)νW ] = 0,

(2.3)

k=1

for B := B(Z), S := S(Z), W = (X10 , W20 )0 ,
0
W2 = ({(D − p(Z))1(Gk )}K
k=1 ) .

The presence of D − p(Z) in the interaction (D − p(Z))1(Gk ) orthogonalizes this regressor relative
to all other regressors that are functions of Z. The controls X1 , e.g. B, can be included to improve
precision.
The second main identification result is that the projection coefficients γk are the GATES
parameters:
K
γ = (γk )K
k=1 = (E[s0 (Z) | Gk ])k=1 .

Given the identification strategy, we can base the corresponding estimation strategy on the following empirical analogs:
Yi = α
b0 X1i + γ
b0 W2i + νbi ,

i ∈ M, EN,M [w(Zi )b
νi Wi ] = 0.

(2.4)

The properties of this estimator, conditional on the auxilliary data, are well known and stated as a
special case of Lemma B.1.
A formal statement appears below, together with a complementary result.

13

1.0

GENERIC ML FOR FEATURES OF HETEROGENOUS TREATMENT EFFECTS

4

●

●
●

●
●

●
●

●

●

●
●

0.0

●
●

E s(Z)|Group

●

0

●
●
●

●

●

●

●
●

●

−0.5

−2

E s(Z)|Group

2

0.5

●
●

−4

●
●

−1.0

●

1

2

3

4

5

1

2

Heterogeneity Group

3

4

5

Heterogeneity Group

Figure 2. In the left panel we have the homogeneous CATE s0 (Z) = 0; in the right
panel we have heterogeneous CATE s0 (Z) = Z. The proxy predictor S(Z) for CATE
is produced by the random forest, whose sorted averages by groups are shown as
red dots, exhibiting large biases. These are the naive estimates. The true sorted
group average treatment effects (GATES) E[s0 (Z) | Gk ] are shown by black dots, and
estimated GATES are shown by blue dots. The true and estimated GATES correct
for the biases relative to the naive strategy shown in red. The estimated GATES
shown by blue dots are always closer to the true GATEs shown by black dots than
the naive estimates shown in red.
Figure 2 provides two examples using the same designs as in fig. 1. Post-processing the ML
estimates again has stronger effect when there is no heterogeneity, but in both cases help bring the
estimated GATES close to the true GATES.
GATES: The Second Strategy. Here we employ linear projections on Horvitz-Thompson transformed variables:
K
X
Y H = µ0 X1 H +
γk · 1(Gk ) + ν, E[ν W̃ ] = 0,
(2.5)
k=1

for B := B(Z), S := S(Z), W̃ = (X10 H, W̃20 ), W̃20 = ({1(Gk )}K
k=1 ).
Again, we show that the projection parameters are GATES:
K
γ = (γk )K
k=1 = (E[s0 (Z) | Gk ])k=1 .

14

V. CHERNOZHUKOV, M. DEMIRER, E. DUFLO, I. FERNANDEZ-VAL

Given the identification strategy, we can base the corresponding estimation strategy on the following empirical analogs:
Yi Hi = µ
b0 X1i Hi + γ
b0 W̃2i + νbi ,

i ∈ M, EN,M [b
νi W̃i ] = 0.

(2.6)

The properties of this estimator, conditional on the auxiliary data, are well known and given in
Lemma B.1. The resulting estimator has similar performance to the previous estimator, and under
some conditions their first-order properties coincide.
The following is the formal statement of the identification result.
Theorem 2.3 (GATES). Consider z 7→ S(z) and z 7→ B(z) as fixed maps. Assume that Y has finite
second moments and the W ’s and W̃ defined above are such that EW W 0 and EW̃ W̃ 0 are finite and have full
rank. Consider γ = (γk )K
k=1 defined by the weighted regression equation (2.3) or by the regression equation
(2.5). These parameters defined in two different ways are equivalent and are equal to the expectation of s0 (Z)
conditional on the proxy group {S ∈ Ik }:
γk = E[s0 (Z) | Gk ].
2.3. Classification Analysis (CLAN). When the BLP and GATES analyses reveal substantial heterogeneity, it is interesting to know the properties of the subpopulations that are most and least
affected. Here we focus on the “least affected group” G1 and “most affect group” GK . Under the
monotonicity assumption, it is reasonable that the first and the last groups are the most and least
affected, where the labels “most” and “least” can be swapped depending on the context.
Let g(Y, Z) be a vector of characteristics of an observational unit. The parameters of interest are
the average characteristics of the most and least affected groups:
δ1 = E[g(Y, Z) | G1 ]

and

δK = E[g(Y, Z) | GK ].

The parameters δK and δ1 are identified because they are averages of variables that are directly
observed. We can compare δK and δ1 to quantify differences between the most and least affected
groups. We call this type of comparisons as classification analysis or CLAN.
3. ”Variational” Estimation and Inference Methods
3.1. Estimation and Inference: The Generic Targets. Let θ denote a generic target parameter or
functional, for example,
• θ = β2 is the heterogeneity predictor loading parameter;
• θ = β1 + β2 (S(z) − ES) is the “personalized” prediction of s0 (z);
• θ = γk is the expectation of s0 (Z) for the group Gk ;
• θ = γK − γ1 is the difference in the expectation of s0 (Z) between the most and least affected
groups;

GENERIC ML FOR FEATURES OF HETEROGENOUS TREATMENT EFFECTS

15

• θ = δK − δ1 is the difference in the expectation of the characteristics of the most and least
impacted groups.
3.2. Quantification of Uncertainty: Two Sources. There are two principal sources of sampling
uncertainty:
(I) Estimation uncertainty regarding the parameter θ, conditional on the data split;
(II) Uncertainty or ”variation” induced by the data splitting.
Conditional on the data split, quantification of estimation uncertainty is standard. To account
for uncertainty with respect to the data splitting, it makes sense to examine the robustness and
variability of the estimates/confidence intervals with respect to different random splits. One of our
goals is to develop methods, which we call ”variational estimation and inference” (VEIN) methods,
for quantifying this uncertainty, which can be of independent interest in many settings where the
sample splitting is used.
Quantifying Source (I): Conditional Inference. We first recognize that the parameters implicitly
depend on
DataA := {(Yi , Di , Xi )}i∈A ,
the auxiliary sample, used to create the ML proxies B = BA and S = SA . Here we make the
dependence explicit: θ = θA .
All of the examples admit an estimator θbA such that under mild assumptions,
2
θbA | DataA ∼a N (θA , σ
bA
),

in the sense that, as |M | → ∞,
−1 b
P(b
σA
(θA − θA ) 6 z | DataA ) →P Φ(z).

Implicitly this requires the auxiliary data DataA to be ”sufficiently regular”, and this should happen with high probability.
As a consequence, the confidence interval (CI)
[LA , UA ] := [θbA ± Φ−1 (1 − α/2)b
σA ]
covers θA with approximate probability 1 − α:
P[θA ∈ [LA , UA ] | DataA ] = 1 − α − oP (1).
This leads to straighforward conditional inference, which does not account for the sample splitting
uncertainty.

16

V. CHERNOZHUKOV, M. DEMIRER, E. DUFLO, I. FERNANDEZ-VAL

Quantifying Source (II): “Variational” Inference. Different partitions (A, M ) of {1, ..., N } yield
different targets θA . Conditional on the data, we treat θA as a random variable, since (A, M ) are
random sets that form random partitions of {1, . . . , N } into samples of size n and N − n. Different partitions also yield different estimators θbA and approximate distributions for these estimators.
Hence we need a systematic way of treating the randomness in these estimators and their distributions.
Comment 3.1. In cases where the data sets are not large, it may be desirable to restrict attention to
balanced partitions (A, M ), where the proportion of treated units is equal to the designed propensity score.
We want to quantify the uncertainty induced by the random partitioning. Conditional on Data,
the estimated θbA is still a random variable, and the confidence band [LA , UA ] is a random set. For
reporting purposes, we instead would like to report an estimator and confidence set, which are
non-random conditional on the data.
Adjusted Point and Interval Estimators. Our proposal is as follows. As a point estimator,
we shall report the median of θbA as (A, M ) vary (as random partitions):
θb := Med[θbA | Data].
This estimator is more robust than the estimator based on a single split. To account for partition uncertainty, we propose to report the following confidence interval (CI) with the nominal
confidence level 1 − 2α:
[l, u] := [Med[LA | Data], Med[UA | Data]].
Note that the price of splitting uncertainty is reflected in the discounting of the confidence level
from 1 − α to 1 − 2α. Alternatively, we can report the confidence interval based on inversion
of a test based upon p-values, constructed below.
The above estimator and confidence set are non-random conditional on the data. The confidence set reflects the uncertainty created by the random partitioning of the data into the main and
auxilliary data.
Comment 3.2. For a random variable X with law PX we define
Med(X) := inf{x ∈ R : PX (X 6 x) > 1/2},
Med(X) := sup{x ∈ R : PX (X > x) > 1/2},
Med(X) := (Med(X) + Med(X))/2.
Note that the lower median is the usual definition of the median. The upper median is the next
distinct quantile of the random variable (or it is the usual median after reversing the order on R).
For example, when X is uniform on {1, 2, 3, 4}, then Med(X) = 2 and Med(X) = 3; and if X is

GENERIC ML FOR FEATURES OF HETEROGENOUS TREATMENT EFFECTS

17

uniform on {1, 2, 3}, then Med(X) = Med(X) = 2. For continuous random variables the upper
and lower medians defined above coincide. For discrete random variables they can differ, but the
differences will be small for variables that are close to being continuous.

Suppose we are testing H0 : θA = θ0 against H1 : θA < θ0 , conditional on the auxiliary data,
then the p-value is given by
pA = Φ(b
σ −1 (θbA − θ0 )).
A

−1 b
The p-value for testing H0 : θA = θ0 against H1 : θA > θ0 , is given by pA = 1 − Φ(b
σA
(θA − θ0 )).

Under the null hypothesis pA is approximately distributed as the uniform variable, pA ∼ U (0, 1),
conditional on DataA . Note that, conditional on Data, pA still has randomness induced by random
partitioning of the data, which we need to address.
Adjusted P-values. We say that testing the null hypothesis, based on the p-values pA , that
are random conditional on data, has significance level α if
P(pA 6 α/2 | Data) > 1/2

or p.5 = Med(pA | Data) 6 α/2.

That is, for at least 50% of the random data splits, the realized p-value pA falls below the level
α/2. Hence we can call p = 2p.5 the sample splitting-adjusted p-value, and consider its small
values as providing evidence against the null hypothesis.
Comment 3.3. Our construction of p-values builds upon the false-discovery-rate type adjustment
ideas in [16, 42], though what we propose is radically simpler, and is minimalistic for our problem,
whereas the idea of our confidence intervals below appears to be completely new.

The main idea behind this construction is simple: the p-values are distributed as marginal uniform variables {Uj }j∈J , and hence obey the following property.
Lemma 3.1 (A Property of Uniform Variables). Consider M , the (usual, lower) median of a sequence
{Uj }j∈J of uniformly distributed variables, Uj ∼ U (0, 1) for each j ∈ J, where variables are not necessarily
independent. Then,
P(M 6 α/2) 6 α.
Proof. Let M denote the median of {Uj }j∈J . Then M 6 α/2 is equivalent to |J|−1
α/2)] − 1/2 > 0. So
X
P[M 6 α/2] = E1{|J|−1
[1(Uj 6 α/2)] > 1/2}.

P

j∈J [1(Uj

6

j∈J

By Markov inequality this is bounded by
X
2E|J|−1
[1(Uj 6 α/2)] 6 2E[1(Uj 6 α/2)] 6 2α/2 = α.
j∈J

where the last inequality holds by the marginal uniformity.



18

V. CHERNOZHUKOV, M. DEMIRER, E. DUFLO, I. FERNANDEZ-VAL

Main Inference Result: Variational P-values and Confidence Intervals. We present a formal result on adjusted p-values using this condition:
PV. Suppose that A is a set of regular auxiliary data configurations such that for all x ∈ [0, 1],
under the null hypothesis:
sup |PP [pA 6 x | DataA ∈ A] − x| 6 δ = o(1),

P ∈P

and inf P ∈P PP [DataA ∈ A] =: 1 − γ = 1 − o(1). In particular, suppose that this holds for
the p-values
−1 b
−1 b
pA = Φ(b
σA
(θA − θA )) and pA = 1 − Φ(b
σA
(θA − θA )).

Lemma B.1 shows that this condition is plausible for the least squares estimators defined in the
previous section under mild conditions.
Theorem 3.1 (Uniform Validity of Variational P-Value). Under condition PV and the null hypothesis
holding,
PP (p.5 6 α/2) 6 α + 2(δ + γ) = α + o(1),
uniformly in P ∈ P.
In order to establish the properties of the confidence interval [l, u], we first consider the properties of the related confidence interval, which is based on the inversion of the p-value based tests:
CI := {θ ∈ R : pu (θ) > α/2, pl (θ) > α/2},

(3.1)

for α < .25 , where, for σ
bA > 0,
−1 b
σA
(θA − θ)] | Data),
pl (θ) := Med(1 − Φ[b

(3.2)

−1 b
pu (θ) := Med(Φ[b
(θA − θ)] | Data).
σA

(3.3)

The confidence interval CI has the following representation in terms of the medians of t-statistics
implied by the proof Theorem 3.2 stated below:


h b
i
θA
−1 (1 − α/2) | Data < 0 

Med θ−
−
Φ
h σbAb
i
CI = θ ∈ R :
.
(3.4)

Med θ−θA + Φ−1 (1 − α/2) | Data > 0 
σ
bA

This CI can be (slightly) tighter than [l, u], while the latter is much simpler to construct.
The following theorem establishes that both confidence sets maintain the approximate confidence level 1 − 2α.
Theorem 3.2 (Uniform Validity of Variational Confidence Intervals). CI can be represented as (3.4)
and CI ⊆ [l, u], and under condition PV,
PP (θA ∈ CI) > 1 − 2α − 2(δ + γ) = 1 − 2α − o(1),
uniformly in P ∈ P.

GENERIC ML FOR FEATURES OF HETEROGENOUS TREATMENT EFFECTS

19

4. Other Considerations and Extensions
1. Choosing the Best ML Method Targeting CATE in Stage 1. There are several options. The
best ML method can be chosen using the auxiliary sample, based on either (A) the ability to predict
Y H using BH and S or (B) the ability to predict Y using B and (D − p(Z))(S − E(S)) under the
weight w(Z) (as in the first type of strategies we developed earlier). To be specific, we can solve
either of the following problems:
(A) minimize the errors in the prediction of Y H on BH and S:
X
(B, S) = arg min
[Yi Hi − B(Zi )Hi − S(Zi )]2 ,
B∈B,S∈S

i∈A

where B and S are parameter spaces for z 7→ B(z) and z 7→ S(z); or
(B) minimize the errors in the weighted prediction of Y on B and (D − p(Z))(S − E(S)):
X
(B, S) = arg min
w(Zi )[Yi − B(Zi ) − (Di − p(Zi )){S(Zi ) − S̄(Zi )}]2 ,
B∈B,S∈S

where S̄(Zi ) = |A|−1
z 7→ S(z).

i∈A

P

i∈A S(Zi )

and B and S are parameter spaces for z 7→ B(z) and

This idea improves over simple but inefficient strategy of predicting Y H just using S, which have
been suggested before for causal inference. It also improves over the simple strategy that predicts
Y using B and DS (which chooses the best predictor for E[Y | D, Z] in a given class but not necessarily the best predictor for CATE s0 (Z)). Note that this idea is new and is of major independent
interest.
2. Choosing the Best ML Method BLP Targeting CATE in Stage 2. The best ML method can
also be chosen in the main sample by maximizing
Λ := |β2 |2 Var(S(Z)) = Corr2 (s0 (Z), S(Z))Var(s0 (Z)).

(4.1)

Maximizing Λ is equivalent to maximizing the correlation between the ML proxy predictor S(Z)
and the true score s0 (Z), or equivalent to maximizing the R2 in the regression of s0 (Z) on S(Z).
3. Choosing the Best ML Method GATES Targeting CATE in Stage 2. Analogously, for GATES
the best ML method can also be chosen in the main sample by maximizing
Λ̄ = E

K
X
k=1

!2
γk 1(S ∈ Ik )

=

K
X

γk2 P(S ∈ Ik ).

(4.2)

k=1

P
This is the part of variation Es20 (Z) of s0 (z) explained by S̄(Z) = K
k=1 γk 1(S(Z) ∈ Ik ). Hence
choosing the ML proxy S(Z) to maximize Λ̄ is equivalent to maximizing the R2 in the regression

20

V. CHERNOZHUKOV, M. DEMIRER, E. DUFLO, I. FERNANDEZ-VAL

of s0 (Z) on S̄(Z) (without a constant). If the groups Gk = {S ∈ Ik } have equal size, namely
P(S(Z) ∈ Ik ) = 1/K for each k = 1, ..., K, then
Λ̄ =

K
1 X 2
γk .
K
k=1

4. Stratified Splitting. The idea is to balance the proportions of treated and untreated in both
A and M samples, so that the proportion of treated is equal to the experiment’s propensity scores
across strata. This formally requires us to replace the i.i.d. assumption by the i.n.i.d. assumption
(independent but not identically distributed observations) when accounting for estimation uncertainty, conditional on the auxiliary sample. This makes the notation more complicated, but the
results in Lemma B.1 still go through with notational modifications.
5. When Proxies have Little Variation. The analysis may generate proxy predictors S that have
little variation, so we can think of them as “weak”, which makes the parameter β2 weakly identified.
We can either add small noise to the proxies, which is called jittering, so that inference results go
through, or we may switch to testing rather than estimation. For practical reasons, we prefer the
jittering approach.

5. Further Potential Applications to Prediction and Causal Inference Problems
Our inference approach generalizes to any problem of the following sort.
Generalization. Suppose we can construct an unbiased signal Ỹ such that
E[Ỹ | Z] = s0 (Z),
where s0 (Z) is now a generic target function. Let S(Z) denote an ML proxy for s0 (Z). Then,
using previous arguments, we immediately can generate the following conclusions:
(1)
(2)
(3)
(4)
(5)

The projection of Ỹ on the ML proxy S(Z) identifies the BLP of s0 (Z) using S(Z).
The grouped averages of the target (GAT) E[s0 (Z) | Gk ] are identified by E[Ỹ | Gk ].
Using ML tools we can train proxy predictors S(Z) to predict Ỹ in auxiliary samples.
We post-process S(Z) in the main sample, by estimating the BLP and GATs.
We apply variational inference on functionals of the BLP and GATs.

The noise reduction strategies, like the ones we used on the context of H-transformed outcomes,
can be useful in these cases, but their construction could depend on the context.
Example 1. Forecasting or Predicting Regression Functions using ML proxies. This is the most
common type of the problem arising in forecasting. Here the target is the best predictor of Y using
Z, namely s0 (Z) = E[Y | Z], and Ỹ = Y trivially serves as the unbiased signal. The interesting part

GENERIC ML FOR FEATURES OF HETEROGENOUS TREATMENT EFFECTS

21

here is the use of variational inference tools developed in this paper for constructing confidence
intervals for the predicted values produced by the estimated BLP of s0 (Z) using S(Z).
Example 2. Predicting Structural Derivatives using ML proxies. Suppose we are interested in
best predicting the partial derivative s0 (z), where s0 (z) = ∂x g(x, z), where g(x, z) = E[Y | X =
x, Z = z]. In the context of demand analysis, Y is the log of individual demand, X is the log-price
of a product, and Z includes prices of other products and characteristics of individuals. Then the
unbiased signal is given by Ỹ = Y [∂x log p(X | Z)], where p(x | z) is the conditional density
function of x given z. That is, E[Ỹ | Z] = s0 (Z) under mild conditions on the density using the
integration by parts formula.

6. Empirical Applications and Implementation Algorithms
We consider two empirical applications. The first application is an RCT conducted in Morocco,
which investigates the effects of microfinance access on several outcomes. The second application
is an observational study, which analyzes the gender wage gap using the U.S. Current Population
Survey (CPS). While the results of the second study should not be interpreted as causal, we aim to
obtain the best approximation to the causal effects of gender-based discrimination on wages using
observational data, by flexibly conditioning on many controls using ML methods. We conclude
this section by providing the implementation algorithms for both examples.
6.1. Heterogeneity in the Effect of Microcredit Availability. We analyze a randomized experiment designed to evaluate the impact of microcredit availability on borrowing and self-employment
activities, which was previously studied in [24]. The experiment was conducted in 162 villages
in Morocco, divided into 81 pairs of villages with similar observable characteristics (number of
households, accessibility to the center of the community, existing infrastructure, type of activities
carried out by the households, and type of agriculture activities). One of the villages in each pair
was randomly assigned to treatment and the other to control. Between 2006 and 2007 a microfinance institution started operating in the treated villages. Two years after the intervention an
endline household survey was conducted with 5,551 households, which constitute our sample.
There was no other microcredit penetration in these villages, before and for the duration of the
study. Therefore, we interpret the treatment as the availability of microcredit.
Recent randomized evaluations of access to microcredit at the community level have found limited impacts of microcredit9. Despite evidence that access to microfinance leads to an increase in
borrowing ([4], [10], [48]) and business creation or expansion ([4], [7], [10], [48]), most studies have
found that this does not translate into an increase in economic outcomes such as profit, income,
labor supply and consumption ([4], [10], [24]). Moreover, there is also no evidence of substantial
gains in human development outcomes, such as education and health ([10],[48]). Studies which
9See [11] for a summary of the recent literature

22

V. CHERNOZHUKOV, M. DEMIRER, E. DUFLO, I. FERNANDEZ-VAL

estimate the impact of microfinance by randomizing microcredit at the individual level confirm
these findings ([8], [36], [37]).
One question that remains elusive is whether the lack of evidence on the average effects masks
heterogeneity, in which there are potential winners and losers of the microcredit expansion. Understanding this heterogeneity can have important implications for evaluating the welfare effects
of microcredit, designing policies and targeting the groups that would benefit from microfinance.
Indeed, the idea that there might be heterogeneity in the impact of microcredit has been a common theme among RCTs evaluating microfinance programs. Having found mostly positive but
insignificant coefficients, the papers cited above attempt to explore heterogeneous treatment effects, mostly using quantile treatment effects. For profits, most studies seem to find positive impact at the higher quantiles (and in the data set we study here, [24] actually find negative impacts at
lower level). Using Bayesian hierarchical methods to aggregate the evidence across studies, Meager (2017) cautions that these results on quantiles may not be generalizable: the profit variables
seems to have too much noise to lend itself to quantile estimation.
A number of recent papers also consider heterogeneous treatment effects by studying the effect
of microfinance on a subpopulation. In a follow-up study of [10], [9] investigates whether the heterogeneity is persistent six years after the microfinance was introduced. They find that credit has a
much bigger impact on the business outcomes of those who started a business before microfinance
entered than of those without prior businesses. Using the same dataset as in this application [24]
classify households into three categories in term of their probability to borrow before the intervention and find that microcredit access has a significant impact on investment and profit, but still no
impact on income and consumption among those who are most likely to borrow. It is worth noting
that the original strategy for this study was to construct groups which, ex ante had different probability to borrow, in order to separately estimate direct effect of microcredit on those most likely
to borrow, and indirect on those very unlikely to borrow. The researchers had initially collected a
short survey on the entire village, and had then tried to predict the probability to borrow using a
model fitted in a first group of villages. But the model proved to have low predictive power in the
other villages, and in their paper then end up predicting the probability to borrow ex-post. In the
paper, they worry that this ex-post classification leads them to over-fit.
The identification strategy developed in this paper provides several advantages in studying heterogeneity in the treatment effects of microfinance. First, contrary to the literature, which relies on
ad hoc subgroup analysis across a few baseline characteristics, we are agnostic about the source of
heterogeneity. While the variable “had a prior business” has proven to be a robust and generalizable predictors of differences in treatment effect (Meager, 2017) and could therefore be pre-specified
in future pre-analysis plans, we have little idea about what else predict heterogeneity. Second, our
approach is valid in high dimensional settings, allowing us to include a rich set of characteristics
in an unspecified functional form. Finally, using the CLAN estimation we are able to identify the

GENERIC ML FOR FEATURES OF HETEROGENOUS TREATMENT EFFECTS

23

characteristics of the most and least affected subpopulations, which could be an important input
for a welfare analysis or targeting households who are likely to benefit from access to microfinance.
We focus on heterogeneity in treatment effects on four household outcome variables, Y : the
amount of money borrowed, the output from self-employment activities, profit from self-employment
activities, and monthly consumption. The treatment variable, D, is an indicator for the household
residing in a treated village. The covariates, Z, include some baseline household characteristics
such as number of members, number of adults, head age, indicators for households doing animal husbandry, doing other non-agricultural activity, having an outstanding loan over the past 12
months, household spouse responded to the survey, another household member (excluding the
HH head) responded to the survey, and 81 village pair fixed effects (these are the variables that
are available for all households). We also include indicators for missing observation at baseline as
controls. Table 6 shows some descriptive statistics for the variables used in the analysis (all monetary variables are expressed in Moroccan Dirams, or MAD). Treated and control households have
similar characteristics and the unconditional average treatment effect on loans, output, profit and
consumption are respectively 1,128, 5,237, 1,844 and -31.

Table 1. Descriptive Statistics of Households
All

Treated

Control

Outcome Variables
Total Amount of Loans

2,359

2,930

1,802

Total output from self-employment activities (past 12 months)

32,499

35,148

29,911

Total profit from self-employment activities (past 12 months)

10,102

11,035

9,191

Total monthly consumption

3,012

2,996

3,027

Number of Household Members

3.879

3.872

3.886

Number of Members 16 Years Old or Older

2.604

2.601

2.607

Head Age

35.976

35.937

36.014

Declared Animal Husbandry Self-employment Activity

0.415

0.426

0.404

Declared Non-agricultural Self-employment Activity

0.146

0.129

0.164

Borrowed from Any Source

0.210

0.224

0.196

Spouse of Head Responded to Self-employment Section

0.067

0.074

0.061

Member Responded to Self-employment Section

0.044

0.048

0.041

Baseline Covariates

We implement our methods using the algorithm described in Section 6.3. By design the propensity score p(Zi ) = 1/2 for all the households. Table 2 compares the four ML methods for producing
the proxy predictors S(Zi ) considered in Stage 1. We find that the Random Forest and Elastic Net

24

V. CHERNOZHUKOV, M. DEMIRER, E. DUFLO, I. FERNANDEZ-VAL

Table 2. Comparison of ML Methods: Microfinance Availability

Elastic Net

Boosting

Neural Network

Random Forest

2,808,960

2,086,256

2,464,276

2,706,767

875

421

630

1253

143,620,159

78,127,786

76,945,520

133,825,537

8,697

2,823

4,219

5,229

32,307,828

17,105,855

20,404,000

39,286,050

4,595

2,319

1,847

4,478

44,583

26,387

34,019

36,927

101

101

96

109

Amount of Loans
Best BLP (Λ)
Best GATES (Λ̄)
Output
Best BLP (Λ)
Best GATES (Λ̄)
Profit
Best BLP (Λ)
Best GATES (Λ̄)
Consumption
Best BLP (Λ)
Best GATES (Λ̄)

Notes: Medians over 100 splits in half.

outperform the Boosted Tree and Neural Network across all outcome variables for both metrics.
Accordingly, we focus on these two methods for the rest of the analysis.10
Table 3 presents results of the BLP of CATE using the ML proxies S(Z) for the four outcome
variables. We report estimates of the coefficients β1 and β2 , which correspond to the ATE and
heterogeneity loading (HET) parameters in the BLP, respectively. In parentheses, we report confidence intervals adjusted for variability across the sample splits using the median method; and
in brackets, we report adjusted p-values. The estimated ATEs of microfinance availability are consistent with the findings of [24] and are similar to the unconditional ATE, as expected by virtue
of the randomization. The ATE on the amount of loans and output are positive and statistically
significant at least at the 10% level with both ML methods. Microfinance availability does not have
a significant impact on profit and consumption.
Turning to the heterogeneity results, we reject the hypothesis that HET is zero at the 10% level
for the amount of loans, output and profit with the elastic net method, suggesting the presence
of heterogeneity in the effect of microfinance availability. The results are consistent across both
ML methods except for output, for which HET coefficient on the Random Forest proxy is not significantly different from zero at the 10%. Finally, the BLP analysis does not reveal any significant
heterogeneity in the effect on consumption. Overall, these results suggest that microfinance availability has heterogenous impacts on business-related outcomes that does not seem to translate into
10The results obtained using Boosted Tree and Neural Network are similar to the results reported, but they are slightly

less precise. These results are not reported but are available from the authors upon requests.

GENERIC ML FOR FEATURES OF HETEROGENOUS TREATMENT EFFECTS

25

Table 3. BLP of Microfinance Availability

Elastic Net

Amount of Loans

Output

Profit

Consumption

Random Forest

ATE (β1 )

HET (β2 )

ATE (β1 )

HET (β2 )

1,163

0.238

1,185

0.375

(544,1736)

(0.021,0.448)

(561,1771)

(0.028,0.774)

[0.000]

[0.060]

[0.000]

[0.069]

5,095

0.262

5,027

0.192

(232,10033)

(0.085,0.433)

(-89,10194)

(-0.100,0.508)

[0.079]

[0.008]

[0.109]

[0.391]

1,553

0.244

1,603

0.279

(-1344,4389)

(0.079,0.416)

(-1276,4536)

(0.046,0.518)

[0.584]

[0.008]

[0.521]

[0.039]

-59.1

0.157

-58.6

0.196

(-161.5,44.2)

(-0.058,0.385)

(-166.6, 43.3)

(-0.160,0.574)

[0.514]

[0.278]

[0.508]

[0.553]

Notes: Medians over 100 splits. 90% confidence interval in parenthesis.
P-values for the hypothesis that the parameter is equal to zero in brackets.

a detectable contemporaneous effect on the standard of living as represented by consumption, even
for the most positively affected households.
We next estimate the GATES. We divide the households into K = 5 groups based on the quintiles
of the ML proxy predictor S(Z) and estimate the average effect for each group. Figures 3-6 present
the estimated GATES coefficients γ1 − γ5 along with joint confidence bands. We also report the
ATE and its confidence interval that were obtained in the BLP analysis for comparison. The GATES
provide a richer understanding of the heterogeneity. In particular, the figures reveal that there are
groups of winners, the most affected groups, for which the GATES on amount of loans, output
and profit are significantly different from zero. These groups are likely to drive the heterogeneity
in the treatment effect that we find in the BLP analysis. We further investigate the GATES by
comparing the most and least affected groups in Table 4. Here we find that the difference of GATES
of these two groups is significantly different from zero at least at the 10% level on amount of loans,
level on output and profit, whereas we fail to reject the hypothesis that this difference is zero at
conventional levels on consumption. Looking at the least affected group, it is reassuring to see that
we have no evidence of negative impact on profit and income, mitigating the concerns that there are
adversely affected households. However, there is negative and insignificant effect on consumption
for the same group. A possible explanation for this result is that investment is lumpy and some
households cut back consumption to increase investment. All the results for the GATES are fairly
robust to the ML method.

26

V. CHERNOZHUKOV, M. DEMIRER, E. DUFLO, I. FERNANDEZ-VAL

Table 4. GATES of 20% Most and Least Affected Groups
Elastic Net

Amount of Loans

Output

Profit

Consumption

Random Forest

20% Most
(γ5 )

20% Least
(γ1 )

Difference
(γ5 − γ1 )

20% Most
(γ5

20% Least
(γ1 )

Difference
(γ5 − γ1 )

2,677

-197

2,995

2,870

94.707

2,814

(1298,4076)

(-1835,1307)

(945,5103)

(1149,4587)

(-1663,1723)

(503,5193)

[0.000]

[1.000]

[0.008]

[0.002]

[1.000]

[0.032]

22,367

-3,039

25,088

21,606

626

21,035

(7678,36920)

(-12546,6535)

(7028,42698)

(5862,38022)

(-11871,13529)

(125,43170)

[0.007]

[1.000]

[0.015]

[0.015]

[1.000]

[0.097]

10,644

-1,152.242

11,768

11,540

-2,031

14,037

(2146,19096)

(-7250,4952)

(1077,22422)

(2965,20955.576)

(-8721,4796)

(2459,25833)

[0.028]

[1.000]

[0.061]

[0.014]

[1.000]

[0.037]

66.4

-333

383

62

-300

332

(-166.2,289.8)

(-695.6,23.2)

(-38.0,805.6)

(-271,346)

(-683,66)

(-196,835)

[1.000]

[0.140]

[0.152]

[1.000]

[0.228]

[0.429]

Notes: Medians over 100 splits. 90% confidence interval in parenthesis.
P-values for the hypothesis that the parameter is equal to zero in brackets.

We conclude by looking at the average characteristics of the most and least affected groups to
understand what generates heterogeneity in the treatment effects. We omit the results for consumption as we do not detect heterogeneity for this outcome. We focus on three characteristics in
this analysis: Age of household head, non-agricultural self-employment activity and whether the
household borrowed from any source. Table 5 reports the CLAN for the 10% least and most affected
groups defined by the deciles of the CATE proxy S(Z) as well as the difference between the two.
We find households with young heads, non-agriculture self-employment and that borrowed less
from any source are likely to borrow more from the microfinance institution, suggesting that formal loans are substitutes rather than complements. For output only the average self-employment
sector indicator is found to be significantly different between the most and least affected groups. Finally, the estimates for profit are similar to the estimates of amount of loans. It is important to note
that the significance of some of these differences is sensitive to the ML method used to generate
the proxy predictors, while the sign of the differences is robust.
6.2. Gender Wage Gap in the U.S.. We consider an application to the gender wage gap using
data from the U.S. March Supplement of the Current Population Survey (CPS) in 2015. We select
white, non-hispanic individuals who are aged 25 to 64 years and work more than 35 hours per week
during at least 50 weeks of the year. We exclude self-employed workers; individuals living in group
quarters; individuals in the military, agricultural or private household sectors; individuals with
inconsistent reports on earnings and employment status; individuals with allocated or missing

GENERIC ML FOR FEATURES OF HETEROGENOUS TREATMENT EFFECTS

27

Table 5. CLAN of Microfinance Availability

Elastic Net

Random Forest

10% Most
(δ10 )

10% Least
(δ1 )

Difference
(δ10 − δ1 )

10% Most
(δ10 )

10% Least
(δ1 )

Difference
(δ10 − δ1 )

29.3
(26.3,32.4)
0.199
(0.159,0.238)
0.144
(0.099,0.189)
-

35.2
(32.2,38.2)
0.068
(0.030,0.108)
0.169
(0.124,0.212)
-

-6.6
(-10.9,-2.4)
[0.004]
0.123
(0.069,0.178)
[0.000]
-0.038
(-0.101,0.025)
[0.448]

23.0
(19.8,26.2)
0.134
(0.096,0.173)
0.109
(0.064,0.153)
-

33.8
(30.5,36.9)
0.118
(0.076,0.156)
0.217
(0.175,0.262)
-

-10.4
(-14.9,-6.0)
[0.000]
0.022
(-0.033,0.075)
[0.875]
-0.107
(-0.164,-0.050)
[0.001]

36.280
(33.4,39.1)
0.275
(0.233,0.315)
0.193
(0.142,0.241)
-

36.708
(33.6,39.6)
0.050
(0.007,0.093)
0.215
(0.167,0.262)
-

-0.896
(-5.242,3.432)
[1.000]
0.226
(0.169,0.285)
[0.000]
-0.033
(-0.102,0.034)
[0.687]

29.090
(25.8,32.3)
0.215
(0.172,0.257)
0.165
(0.121,0.208)
-

30.831
(27.5,34.1)
0.088
(0.045,0.129)
0.189
(0.146,0.234)
-

-1.925
(-6.648,2.799)
[0.849]
0.130
(0.070,0.190)
[0.000]
-0.024
(-0.086,0.039)
[0.895]

34.1
(31.2,37.0)
0.181
(0.140,0.222)
0.180
(0.130,0.230)
-

40.4
(37.5,43.4)
0.108
(0.068,0.149)
0.257
(0.207,0.307)
-

-6.5
(-10.7,-2.5)
[0.003]
0.082
(0.022,0.138)
[0.014]
-0.091
(-0.160,-0.022)
[0.020]

29.2
(25.7,32.6)
0.153
(0.113,0.192)
0.144
(0.098,0.190)
-

33.7
(30.390,37.108)
0.099
(0.058,0.139)
0.162
(0.122,0.206)
-

-5.8
(-10.566,-1.217)
[0.029]
0.051
(-0.003,0.105)
[0.129]
-0.032
(-0.095,0.029)
[0.578]

Amount of Loans
Head Age

Non-agricultural self-emp.

Borrowed from Any Source

Output
Head Age

Non-agricultural self-emp.

Borrowed from Any Source

Profit
Head Age

Non-agricultural self-emp.

Borrowed from Any Source

Notes: Medians over 100 splits. 90% confidence interval in parenthesis.
P-values for the hypothesis that the parameter is equal to zero in brackets.

information in any of the variables used in the analysis; and individuals with hourly wage rate
below $3. The resulting sample consists of 32, 523 workers including 18, 137 men and 14, 382 of
women.
We will estimate the BLP, GATES and CLAN of the conditional gender wage gap, i.e., the difference in the average wage between female and male workers with the same observable characteristics; see [18] for a recent survey on this topic. The outcome variable Y is the logarithm of
the hourly wage rate constructed as the ratio of the annual earnings to the total number of hours
worked, which is constructed in turn as the product of number of weeks worked and the usual

28

V. CHERNOZHUKOV, M. DEMIRER, E. DUFLO, I. FERNANDEZ-VAL

number of hours worked per week. The treatment D is an indicator for female worker, and the observable characteristics Z include 5 marital status indicators (widowed, divorced, separated, never
married, and married); 5 educational attainment indicators (less than high school graduate, high
school graduate, some college, college graduate, and advanced degree); 4 region indicators (midwest, south, west, and northeast); potential experience constructed as the maximum of age minus
years of schooling minus 7 and zero, i.e., experience = max(age − education − 7, 0); 5 occupation
indicators (management, professional and related; service; sales and office; natural resources, construction and maintenance; and production, transportation and material moving); and 12 industry
indicators (mining, quarrying, and oil and gas extraction; construction; manufacturing; wholesale
and retail trade; transportation and utilities; information; financial services; professional and business services; education and health services; leisure and hospitality; other services; and public
administration).11
Table 6. Descriptive Statistics of Workers
All

Women

Men

All

Women

Men

Log wage

3.16

3.02

3.26

O.manager

0.49

0.56

0.43

Female

0.44

1.00

0.00

O.service

0.09

0.10

0.09

MS.married

0.70

0.65

0.73

O.sales

0.22

0.30

0.16

MS.widowed

0.01

0.02

0.01

O.construction

0.09

0.01

0.15

MS.separated

0.02

0.02

0.01

O.production

0.11

0.04

0.17

MS.divorced

0.12

0.15

0.09

I.minery

0.03

0.01

0.04

MS.Nevermarried

0.16

0.16

0.16

I.construction

0.06

0.01

0.09

E.lhs

0.02

0.01

0.03

I.manufacture

0.13

0.08

0.18

E.hsg

0.25

0.21

0.28

I.retail

0.13

0.11

0.14

E.sc

0.28

0.30

0.27

I.transport

0.04

0.02

0.06

E.cg

0.28

0.29

0.27

I.information

0.02

0.02

0.03

E.ad

0.17

0.19

0.15

I.finance

0.08

0.10

0.07

R.northeast

0.26

0.26

0.26

I.professional

0.11

0.10

0.12

R.midwest

0.27

0.28

0.27

I.education

0.25

0.41

0.12

R.south

0.33

0.34

0.33

I.leisure

0.04

0.05

0.04

R.west

0.22

0.21

0.23

I.services

0.03

0.03

0.04

Experience

21.21

21.03

21.35

I.public

0.07

0.07

0.08

Source: March Supplement CPS 2015

Table 6 reports sample means of the variables used in the analysis. Working women are more
highly educated than working men, have about the same potential experience, and are less likely
to be married and more likely to be divorced. They work relatively more often in managerial and
sales occupations and in the industries providing education and health services. Working men
are relatively more likely to work in construction and production occupations within non-service
industries. The unconditional gender wage gap is 24%.
11The sample selection and variable construction is the same as in [20].

GENERIC ML FOR FEATURES OF HETEROGENOUS TREATMENT EFFECTS

29

In this application we do not observe the propensity scores. We estimate them from a large
sample of 129,983 workers from the U.S. March Supplement of the CPS pooling the years 2012 to
2017, excluding 2015, as described in Section 6.3. We ignore the sampling error in the estimated
propensity scores for the rest of the analysis because the sample size to estimate these scores is
much larger than the sample size to estimate the parameters of interest.
Table 7. Comparison of ML Methods

Elastic Net Boosting Neural Network Random Forest
Best BLP (Λ)

0.084

0.079

0.076

0.082

Best GATES (Λ̄)

0.060

0.058

0.058

0.059

Notes: Medians over 1,000 splits in half.

Table 7 compares the ML methods. In this case the elastic net comes out as the winner both
based on BLP and GATES targeting of CATE, followed closely by the random forest. As in the
previous section, we focus on these two ML methods for the rest of the analysis. Table 8 shows the
results for the BLP coefficients. The average conditional gender wage gap of −22% comes close to
the unconditional gap, but the slope of the BLP uncovers substantial heterogeneity across workers.
The estimates and joint confidence intervals for the GATES in Figure 7 confirm this finding. Here
we divide the sample in K = 5 groups defined by the quintiles of the CATE proxy S(Z) and show
that average conditional gender wage gap ranges from about −33% to −9% across the groups.
Table 9 compares the averages of the gender wage gap in the first (most affected) and fifth (least
affected) groups. We reject that the average gaps for the 20% most and least affected are equal at
any conventional level. The results for the BLP and GATES are robust to the ML method used in
Stage 1.
Table 8. BLP of Gender Wage Gap

Elastic Net
ATE (β1 )

HET (β2 )

Random Forest
ATE (β1 )

-0.229
0.726
-0.228
(-0.251,-0.207) (0.540,0.912) (-0.249,-0.207)
[0.000]
[0.000]
[0.000]

HET (β2 )
0.552
(0.409,0.694)
[0.000]

Notes: Medians over 1,000 splits. 90% confidence interval in parenthesis.
P-values for the hypothesis that the parameter is equal to zero in brackets.

Table 10 reports the results of a classification analysis where we divide the workers in K = 10
groups defined by the deciles of the CATE proxy S(Z) and compare the average characteristics of

30

V. CHERNOZHUKOV, M. DEMIRER, E. DUFLO, I. FERNANDEZ-VAL

Table 9. GATES of Gender Wage Gap

20% Least Affected 20% Most Affected Difference

Elastic Net

Random Forest

(γ5 )

(γ1 )

(γ5 − γ1 )

-0.088

-0.328

0.240

(-0.136,-0.040)

(-0.377,-0.279)

(0.171,0.308)

[0.001]

[0.000]

[0.000]

-0.095

-0.329

0.234

(-0.142,-0.048)

(-0.377,-0.280)

(0.166,0.301)

[0.000]

[0.000]

[0.000]

Notes: Medians over 1,000 splits. 90% confidence interval in parenthesis.
P-values for the hypothesis that the parameter is equal to zero in brackets.

the most and least affected group. This analysis reveals another dimension of the heterogeneity in
the gender wage gap by showing that the 10% most affected workers are more likely to have either
of the following characteristics relative to the 10% less affected workers on average: high wages,
high experience, not advanced degree, married, not divorced or never married. These differences
are robust to the two ML method. The elastic net and random forest, however, disagree in the
comparison of the averages of some of the education variables such as high school graduate and
some college. Due to this sensitivity, we recommend reporting results from at least the two best
ML methods. Overall, the heterogeneity patterns that we find are consistent with both glass ceiling
theories for the gender wage gap [3], and with heterogenous individual preferences where young,
single, and highly educated women are more career-oriented. [20] found similar patterns applying
semiparametric regression methods to the same data.
6.3. Implementation Algorithm. In this section we describe an algorithm based on the first identification strategy and provide some specific implementation details for the two empirical examples.
Algorithm 1 (Inference Algorithm). The inputs are given by the data on units i ∈ [N ] = {1, ..., N }.
Step 0. Fix the number of splits S and the significance level α, e.g. S = 100 and α = 0.05.
Step 1. Compute the propensity scores p(Zi ) for i ∈ [N ].
Step 2. Consider S splits in half of the indices i ∈ {1, ..., N } into the main sample, M , and the
auxiliary sample, A. Over each split s = 1, .., S, apply the following steps:
a. Tune and train each ML method separately to learn B(·) and S(·) using A. For each i ∈ M ,
compute the predicted baseline effect B(Zi ) and predicted treatment effect S(Zi ). If there is
zero variation in B(Zi ) and S(Zi ) add Gaussian noise with a variance of 0.1 to the proxies.
b. Estimate the BLP parameters by weighted OLS in M , i.e.,
Yi = α
b0 X1i + βb1 (Di − p(Zi )) + βb2 (Di − p(Zi ))(Si − EN,M Si ) + b
i , i ∈ M

GENERIC ML FOR FEATURES OF HETEROGENOUS TREATMENT EFFECTS

31

Table 10. CLAN of Gender Wage Gap

Elastic Net

Log Wage

Experience

Random Forest

10% Least

10% Most

Difference

10% Least

10% Most

Difference

(δ10 )

(δ1 )

(δ10 − δ1 )

(δ10 )

( δ1 )

(δ10 − δ1 )

2.949
(2.918,2.981)
-

3.436
(3.404,3.467)
-

-0.484
(-0.528,-0.440)
[0.000]

3.074
(3.043,3.106)
-

3.233
(3.201,3.264)
-

-0.157
(-0.201,-0.113)
[0.000]

12.902
25.664
-12.778
14.176
25.109
-10.893
(12.386,13.419) (25.14,26.187) (-13.513,-12.055) (13.666,14.687) (24.602,25.622) (-11.615,-10.176)
[0.000]
[0.000]

Less than H. School

0.039
(0.030,0.048)
-

0.023
(0.014,0.032)
-

0.017
(0.004,0.028)
[0.016]

0.025
(0.016,0.034)
-

0.039
(0.030,0.047)
-

-0.013
(-0.025,-0.001)
[0.078]

High School

0.246
(0.226,0.266)
-

0.167
(0.147,0.187)
-

0.080
(0.052,0.109)
[0.000]

0.162
(0.141,0.182)
-

0.268
(0.248,0.288)
-

-0.106
(-0.134,-0.078)
[0.000]

Some College

0.184
(0.166,0.201)
-

0.094
(0.076,0.110)
-

0.089
(0.066,0.112)
[0.000]

0.186
(0.166,0.206)
-

0.238
(0.217,0.258)
-

-0.050
(-0.079,-0.021)
[0.001]

College Graduate

0.326
(0.303,0.35)
-

0.570
(0.546,0.594)
-

-0.235
(-0.269,-0.202)
[0.000]

0.361
(0.337,0.384)
-

0.354
(0.33,0.377)
-

0.010
(-0.023,0.043)
[1.000]

Advanced Degree

0.201
(0.182,0.22)
-

0.135
(0.116,0.154)
-

0.066
(0.04,0.091)
[0.000]

0.259
(0.241,0.278)
-

0.096
(0.078,0.115)
-

0.164
(0.137,0.19)
[0.000]

Married

0.101
(0.088,0.114)
-

0.953
(0.94,0.966)
-

-0.851
(-0.869,-0.832)
[0.000]

0.359
(0.338,0.38)
-

0.861
(0.84,0.881)
-

-0.499
(-0.529,-0.47)
[0.000]

Widowed

0.010
(0.005,0.015)
-

0.008
(0.003,0.012)
-

0.003
(-0.004,0.009)
[0.728]

0.008
(0.004,0.013)
-

0.011
(0.006,0.016)
-

-0.002
(-0.009,0.005)
[0.981]

Divorced

0.092
(0.08,0.103)
-

0.019
(0.008,0.03)
-

0.072
(0.056,0.088)
[0.000]

0.095
(0.081,0.109)
-

0.083
(0.069,0.097)
-

0.012
(-0.008,0.031)
[0.478]

Separated

0.028
(0.021,0.035)
-

0.010
(0.003,0.017)
-

0.017
(0.008,0.027)
[0.001]

0.013
(0.006,0.019)
-

0.021
(0.014,0.027)
-

-0.008
(-0.017,0.001)
[0.170]

Nevermarried

0.764
(0.749,0.779)
-

0.003
(-0.011,0.018)
-

0.76
(0.738,0.781)
[0.000]

0.523
(0.505,0.541)
-

0.022
(0.004,0.04)
-

0.500
(0.474,0.526)
[0.000]

Notes: Medians over 1,000 splits. 90% confidence interval in parenthesis.
P-values for the hypothesis that the parameter is equal to zero in brackets.

0 , D − p(Z ), (D − p(Z ))(S − E
0
such that EN,M [w(Zi )b
i Xi ] = 0 for Xi = [X1i
i
i
i
i
i
N,M Si )] , where
w(Zi ) = {p(Zi )(1 − p(Zi ))}−1 and X1i includes a constant, B(Zi ) and S(Zi ).

32

V. CHERNOZHUKOV, M. DEMIRER, E. DUFLO, I. FERNANDEZ-VAL

c. Estimate the GATES parameters by weighted OLS in M , i.e.,
0

Yi = α
b X1i +

K
X

γ
bk · (Di − p(Zi )) · 1(Si ∈ Ik ) + νbi , i ∈ M,

k=1
0 , {(D − p(Z ))1(S ∈ I )}K ]0 , where w(Z ) =
such that EN,M [w(Zi )b
νi Wi ] = 0 for Wi = [Xi1
i
i
i
i
k k=1
{p(Zi )(1 − p(Zi ))}−1 , X1i includes a constant, B(Zi ) and S(Zi ), Ik = [`k−1 , `k ), and `k is the
(k/K)-quantile of {Si }i∈M .
d. Estimate the CLAN parameters in M by

δb1 = EN,M [g(Yi , Zi ) | Si ∈ I1 ]

and

δbK = EN,M [g(Yi , Zi ) | Si ∈ IK ],

where Ik = [`k−1 , `k ) and `k is the (k/K)-quantile of {Si }i∈M .
e. Compute the two performance measures for the ML methods
b = |βb2 |2 Var(S(Z))
d
Λ

K
X
b̄ = 1
γ
bk2 .
Λ
K
k=1

b̄ over the splits.
b and Λ
Step 3: Choose the best ML methods based on the medians of Λ
Step 4: Compute the estimates, (1 − α)-level conditional confidence intervals and conditional
p-values for all the parameters of interest. Monotonize the confidence intervals if needed. For
example, construct a (1 − α) joint confidence interval for the GATES as
{b
γk ± b
c(1 − α)b
σk , k = 1, . . . , K},

(6.1)

where b
c(1 − α) is a consistent estimator of the (1 − α)-quantile of maxk∈1,...,K |b
γk − γk |/b
σk and σ
bk
is the standard error of γ
bk conditional on the data split. Monotonize the band (6.1) with respect to
k using the rearrangement method of [19].
Step 5: Compute the adjusted (1 − 2α)-confidence intervals and adjusted p-values using the
VEIN methods described in Section 3.
Comment 6.1 (ML Methods). We consider four ML methods to estimate the proxy predictors:
elastic net, boosted trees, neural network with feature extraction, and random forest. The ML
methods are implemented in R using the package caret [39]. The names of the elastic net, boosted
tree, neural network with feature extraction, and random forest methods in caret are glmnet, gbm,
pcaNNet and rf, respectively. For each split of the data, we choose the tuning parameters separately
for B(z) and S(z) based on mean squared error estimates of repeated 2-fold cross-validation, except
for random forest, for which we use the default tuning parameters to reduce the computational
time.12 In tuning and training the ML methods we use only the auxiliary sample. In all the methods
we rescale the outcomes and covariates to be between 0 and 1 before training.
12We have the following tuning parameters for each method: Elastic Net: alpha (Mixing Percentage), lambda (Reg-

ularization Parameter), Boosted trees: n.trees (Number of Boosting Iterations), interaction.depth (Max Tree Depth),
shrinkage (Shrinkage), n.minobsinnode (Min. Terminal Node Size), size (Number of Hidden Units) , decay (Weight
Decay), mtry (Number of Randomly Selected Predictors).

GENERIC ML FOR FEATURES OF HETEROGENOUS TREATMENT EFFECTS

33

Comment 6.2 (Microfinance Application). We adopt two strategies to improve precision. First,
the linear projections of the BLP and GATES control for village pair fixed effects along with the
predicted baseline effect, B(z) and predicted treatment effect, S(z). Second, as suggested in Section
4, we use stratified sample splitting where the strata are village pairs. We cluster the standard errors
at the village level to account for potential correlated shocks within each village. All reported
results are medians over S = 100 splits and α = 0.05.
Comment 6.3 (Gender Wage Gap Application). We estimate the propensity scores using a neural
network with feature extraction from a large sample of 129,983 workers from the U.S. March Supplement of the CPS pooling the years 2012 to 2017, excluding 2015. We select the neural network
among elastic net, boosted tree, neural network, and random forest, and the tuning parameters by
repeated 2-fold cross validation on the mean squared error. To avoid tuning ML methods for each
data split separately, we tune ML methods by using 20% random extract from the 2012-2017 pooled
sample and use those tuning parameters for all splits. This allows us to use S = 1, 000 splits in this
empirical application. Finally, to reduce the disproportionate impact of extreme propensity score
weights we trimmed the predicted propensity scores between 0.01 and 0.99 in the main sample.
These dropped observations make up approximately 4% of the sample.

7. Concluding Remarks
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 new approach is valid in high dimensional settings, where the effects are
estimated by machine learning methods. The main advantage of our approach is its credibility: the
approach is agnostic about the properties of the machine learning estimators, and does not rely on
incredible or hard-to-verify assumptions. Estimation and inference relies on data splitting, where
the latter allows us to avoid overfitting and all kinds of non-regularities. Our inference quantifies uncertainty coming from both parameter estimation and the data splitting, and could be of
independent interest. Empirical applications illustrate the practical uses of the approach.

Appendix A. Proofs
Proof of Theorem 2.1. The subset of the the normal equations, which correspond to β1 and β2 , are
given by E[w(Z)(Y − α0 X1 − β 0 X2 )X2 ] = 0. Substituting Y = b0 (Z) + s0 (Z)D + U , and using the
definition X2 = X2 (Z, D) = [D − p(Z), (D − p(Z)(S − ES)]0 , that X1 = X1 (Z), and the law of

34

V. CHERNOZHUKOV, M. DEMIRER, E. DUFLO, I. FERNANDEZ-VAL

iterated expectations, we notice that:
E[w(Z)b0 (Z)X2 ] = E[w(Z)b0 (Z) E[X2 | Z]] = 0,
=0

E[w(Z)U X2 ] = E[w(Z) E[U | Z, D] X2 (Z, D)] = 0,
0

E[w(Z)X1 X2 ] = E[w(Z)X1 (Z) E[X2 (Z, D) | Z]] = 0.
=0

Hence the normal equations simplify to: E[w(Z)(s0 (Z)D − β 0 X2 )X2 ] = 0. Since
E[{D − p(Z)}{D − p(Z)} | Z] = p(Z)(1 − p(Z)) = w−1 (Z),
and S = S(Z), the components of X2 are orthogonal by the law of iterated expectations:
Ew(Z)(D − p(Z))(D − p(Z))(S − ES) = E(S − ES) = 0.
Hence the normal equations above further simplify to
E[w(Z){s0 (Z)D − β1 (D − p(Z))}(D − p(Z))] = 0,
E[w(Z){s0 (Z)D − β2 (D − p(Z))(S − ES)}(D − p(Z))(S − ES)] = 0.
Solving these equations and using the law of iterated expectations, we obtain
β1 =

Ew(Z)s0 (Z)w−1 (Z)
Ew(Z){s0 (Z)D(D − p(Z))}
=
= Es0 (Z),
Ew(Z)(D − p(Z))2
Ew(Z)w−1 (Z)

Ew(Z){s0 (Z)D(D − p(Z))(S − ES)}
Ew(Z)(D − p(Z))2 (S − ES)2
Ew(Z)s0 (Z)w−1 (Z)(S − ES)
=
= Cov(s0 (Z), S)/Var(S).
Ew(Z)w−1 (Z)(S − ES)2
The conclusion follows by noting that these coefficients also solve the normal equations
β2 =

E{[s0 (Z) − β1 − β2 (S − ES)][1, (S − ES)]0 } = 0,
which characterize the optimum in the problem of best linear approximation/prediction of s0 (Z)
using S.

Proof of Theorem 2.2. The normal equations defining β = (β10 , β20 ) are given by E[(Y H − µ0 X1 H −
β 0 X̃2 )X̃2 ] = 0. Substituting Y = b0 (Z) + s0 (Z)D + U , and using the definition X̃2 = X̃2 (Z) =
[1, (S(Z) − ES(Z))]0 , that X1 = X1 (Z), and the law of iterated expectations, we notice that:
E[b0 (Z)H X̃2 (Z)] = E[b0 (Z) E[H | Z] X̃2 (Z)] = 0,
=0

E[U H X̃2 (Z)] = E[E[U | Z, D] H(D, Z)X̃2 (Z)] = 0,
0

E[X1 (Z)H X̃2 (Z)] = E[X1 (Z) E[H | Z] X̃2 (Z)] = 0.
=0

Hence the normal equations simplify to:
E[(s0 (Z)DH − β 0 X̃2 )X̃2 ] = 0.

GENERIC ML FOR FEATURES OF HETEROGENOUS TREATMENT EFFECTS

35

Since 1 and S − ES are orthogonal, the normal equations above further simplify to
E{s0 (Z)DH − β1 } = 0,
E[{s0 (Z)DH − β2 (S − ES)}(S − ES)] = 0.
Using that
E[DH | Z] = [p(Z)(1 − p(Z))]/[p(Z)(1 − p(Z))] = 1,
S = S(Z), and the law of iterated expectations, the equations simplify to
E{s0 (Z) − β1 } = 0,
E{s0 (Z) − β2 (S − ES)}(S − ES) = 0.
These are normal equations that characterize the optimum in the problem of best linear approximation/prediction of s0 (Z) using S. Solving these equations gives the expressions for β1 and β2
stated in the theorem.

Proof of Theorem 2.3. The proof is similar to the proof of Theorem 2.1- 2.2. Moreover, since the
proofs for the two strategies are similar, we will only demonstrate the proof for the second strategy.
The subset of the normal equations, which correspond to γ = (γk )K
k=1 , are given by E[(Y H −
0
0
µ W̃1 − β W̃2 )W̃2 ] = 0. Substituting Y = b0 (Z) + s0 (Z)D + U , and using the definition W̃2 =
0
W̃2 (Z) = [1(S ∈ Ik )K
k=1 ] , that W̃1 = X1 (Z)H, and the law of iterated expectations, we notice that:
E[b0 (Z)H W̃2 (Z)] = E[b0 (Z) E[H | Z] W̃2 (Z)] = 0,
=0

E[U H W̃2 (Z)] = E[E[U | Z, D] H(D, Z)W̃2 (Z)] = 0,
0

E[W̃1 W̃2 (Z)] = E[X1 (Z) E[H | Z] W̃2 (Z)] = 0.
=0

Hence the normal equations simplify to:
E[{s0 (Z)DH − β 0 W̃2 }W̃2 ] = 0.
0
Since components of W̃2 = W̃2 (Z) = [1(Gk )K
k=1 ] are orthogonal, the normal equations above further simplify to

E[{s0 (Z)DH − γk 1(Gk )}1(Gk )] = 0.
Using that
E[DH | Z] = [p(Z){1 − p(Z)}]/[p(Z){1 − p(Z)}] = 1,
that S = S(Z), and the law of iterated expectations, the equations simplify to
E[{s0 (Z) − γk 1(Gk )}1(Gk )] = 0 ⇐⇒
γk = Es0 (Z)1(Gk )/E[1(Gk )] = E[s0 (Z) | Gk ].


36

V. CHERNOZHUKOV, M. DEMIRER, E. DUFLO, I. FERNANDEZ-VAL

Proof of Theorem 3.1. We have that p.5 6 α/2 is equivalent to EP [1(pA 6 α/2) | Data] > 1/2. So
PP [p.5 6 α/2] = EP 1{EP [1(pA 6 α/2) | Data] > 1/2}.
By Markov inequality,
EP 1{EP [1(pA 6 α/2) | Data] > 1/2} 6 2PP [pA 6 α/2]
Moreover,
PP (pA 6 α/2) 6 EP [PP [pA 6 α/2 | DataA ∈ A] + γ] 6 α/2 + δ + γ.


Proof of Theorem 3.2. To show the second claim, we note that
PP (θA 6∈ CI) = PP (pl (θA ) 6 α/2) + PP (pu (θA ) 6 α/2)
6 α + δ + γ + α + δ + γ,
where the inequality holds by Theorem 3.1 on the p-values. The last bound is upper bounded by
2α + o(1) by the regularity condition PV for the p-values, uniformly in P ∈ P.
To show the first claim, we need to show the following inequalities:
sup{θ ∈ R : pu (θ) > α/2} 6 u,

inf{θ ∈ R : pl (θ) > α/2} > l.

We demonstrate the first inequality, and the second follows similarly.
We have that
−1 b
(θA − θ)} | Data] > α/2}
σA
{θ ∈ R : pu (θ) > α/2} = {θ ∈ R : Med[Φ{b
−1 b
= {θ ∈ R : Φ{Med[b
(θA − θ) | Data]} > α/2}
σA
−1 b
= {θ ∈ R : Med[b
(θA − θ) | Data] > Φ−1 (α/2)}
σA
−1
(θ − θbA ) | Data] < Φ−1 (1 − α/2)}
= {θ ∈ R : Med[b
σA
(
"
#
)
θ − θbA
−1
=
θ ∈ R : Med
− Φ (1 − α/2) | Data < 0 ,
σ
bA

where we have used the equivariance of Med and Med to monotone transformations, implied from
their definition. We claim that by the definition of
u := Med[θbA + σ
bA Φ−1 (1 − α/2) | Data],
we have

"

#
u − θbA
−1
Med
− Φ (1 − α/2) | Data > 0.
σ
bA

Indeed, by the definition of u we have that


E 1(u − θbA − σ
bA Φ−1 (1 − α/2) > 0) | Data > 1/2.

GENERIC ML FOR FEATURES OF HETEROGENOUS TREATMENT EFFECTS

37

Since σ
bA > 0 by assumption we have
1(u − θbA − σ
bA Φ

−1

(1 − α/2) > 0) = 1

!
u − θbA
−1
− Φ (1 − α/2) > 0 ,
σ
bA

and it follows that
P

 u − θb
A
σ
bA


− Φ−1 (1 − α/2) > 0 | Data > 1/2.

The claimed inequality sup{θ ∈ R : pu (θ) > α/2} 6 u follows.



Appendix B. A Lemma on Uniform in P Conditional Inference
Lemma B.1. Fix positive constants c and C, and a small constant δ > 0. Let Ỹ and X denote generic
outcome and generic d-vector of regressors, whose use and definition may differ in different places of the
paper. Assume that for each P ∈ P, EP |Ỹ |4+δ < C and let 0 < w 6 w(Z) 6 w < ∞ denote the generic
weights, and that {(Ỹi , Zi , Di )}N
i=1 are i.i.d. copies of (Ỹ , Z, D). Let {DataA ∈ AN } be the event such
that the ML algorithm, operating only on DataA , produces a vector XA = X(Z, D; DataA ) that obeys, for
A = Ỹ − X 0 βA defined by: EP [A w(Z)XA | DataA ] = 0, the following inequalities, uniformly in P ∈ P
0
0
EP [kXA k4+δ | DataA ] 6 C, mineig EP [XA XA
| DataA ] > c, mineig EP [2A XA XA
| DataA ] > c.

Suppose that PP {DataA ∈ AN } > 1 − γ → 1 uniformly in P ∈ P, as N → ∞. Let βbA be defined by:
EN,M [w(Z)XAb
A ] = 0,

b
A = YA − X 0 βbA .

0 −1 be an estimator of
0 (E
0 )−1 E
2A XA XA
Let VbN,A := (EN,M XA XA
N,M XA XA )
N,M b
0
0
0
VN,A = (EP [XA XA
| DataA ])−1 EP [2A XA XA
| DataA ](EP [XA XA
| DataA ])−1 .

Then for any convex set R in Rd , we have that uniformly in P ∈ P:
−1/2

PP [VbN,A (βbA − βA ) ∈ R | DataA ] →P P(N (0, I) ∈ R),
−1/2
PP [VbN,A (βbA − βA ) ∈ R | {DataA ∈ AN }] → P(N (0, I) ∈ R),

and the same results hold with VbN,A replaced by VN,A .
Proof. It suffices to demonstrate the argument for an arbitrary sequence {Pn } in P. Let z 7→
X̃A,N (z) be a deterministic map such that the following inequalities hold, for ẽA defined by
EPn [ẽA w(Z)X̃A,N (Z)] = 0
and X̃A = X̃A (Z):
0
0
EPn [kX̃A,N k4 ] < C, mineig EPn [X̃A,N X̃A,N
] > c, mineig EPn [ẽ2A X̃A,N X̃A,N
] > c.

Let d := dim βA . Then we have that (abusing notation):
BN := sup

sup

X̃A,N h∈BL1

(Rd )

−1/2

|EPn h(ṼN,A (βbA − βA ) | X̃A,N ) − Eh(N (0, I))| → 0,

38

V. CHERNOZHUKOV, M. DEMIRER, E. DUFLO, I. FERNANDEZ-VAL

by the standard argument for asymptotic normality of the least squares estimator, which utilizes
the Lindeberg-Feller Central Limit Theorem. Here
0 −1
0
0 −1
) ,
(EX̃A X̃A
) E˜
2A X̃A X̃A
ṼN,A := (EX̃A X̃A

and BL1 (Rd ) denotes the set of Lipschitz maps h : Rd → [0, 1] with the Lipschitz coefficient
bounded by 1.
Then, for the stochastic sequence XA,N = XA,N (DataA ), we have
sup
h∈BL1

(Rd )

−1/2

|EPn [h(VN,A (βbA − βA )) | XA,N ] − E[h(N (0, I))]| 6 BN + 2(1 − 1{DataA ∈ AN }) →Pn 0.

1/2 −1/2
Since under the stated bounds on moments, VbN,A VN,A →Pn I by the standard argument for consistency of the Eicker-Huber-White sandwich, we further notice that

sup
h∈BL1 (Rd )

−1/2
−1/2
|EPn [h(VbN,A (βbA − βA )) | XA,N ] − EPn [h(VN,A (βbA − βA )) | XA,N ]|

−1/2 1/2
−1/2
6 EPn [kVbN,A VN,A − Ik ∧ 1 · kVN,A (βbA − βA )k ∧ 1 | XA,N ] →Pn E[0 ∧ 1 · kN (0, I)k ∧ 1] = 0,

in order to conclude that
sup
h∈BL1

(Rd )

−1/2
EPn [h(VbN,A (βbA − βA )) | XA,N ] − E[h(N (0, I))] →P 0.

−1/2
−1/2
Moreover, since EPn [h(VbN,A (βbA − βA )) | XA,N ] = EPn [h(VbN,A (βbA − βA )) | DataA ], the first
−1/2
conclusion follows: PPn [VbN,A (βbA − βA ) ∈ R | DataA ] →Pn P(N (0, I) ∈ R), by the conventional
smoothing argument (where we approximate the indicator of a convex region by a smooth map
with finite Lipschitz coefficient). The second conclusion
−1/2
PPn [VbN,A (βbA − βA ) ∈ R | DataA ∈ AN ] → P(N (0, I) ∈ R)

follows from the first by
−1/2
PPn [VbN,A (βbA − βA ) ∈ R | DataA ∈ AN ] =
−1/2
= EPn [PPn [VbN,A (βbA − βA ) ∈ R | DataA ]1({DataA ∈ AN })/PPn {DataA ∈ AN }]

→ E[P(N (0, I) ∈ R) · 1],
using the definition of the weak convergence, implied by the convergence to the constants in probability.


GENERIC ML FOR FEATURES OF HETEROGENOUS TREATMENT EFFECTS

39

References
[1] Alberto Abadie. Semiparametric difference-in-differences estimators. The Review of Economic Studies, 72(1):1–19,
2005.
[2] Alberto Abadie, Matthew M Chingos, and Martin R West. Endogenous stratification in randomized experiments.
Technical report, National Bureau of Economic Research, 2017.
[3] James Albrecht, Anders Bjorklund, and Susan Vroman. Is there a glass ceiling in sweden? Journal of Labor Economics,
21(1):145–177, 2003.
[4] Manuela Angelucci, Dean Karlan, and Jonathan Zinman. Microcredit impacts: Evidence from a randomized microcredit program placement experiment by compartamos banco. American Economic Journal: Applied Economics,
7(1):151–182, 2015.
[5] Susan Athey and Guido Imbens. Recursive partitioning for heterogeneous causal effects. Proceedings of the National
Academy of Sciences, 113(27):7353–7360, 2016.
[6] Susan Athey and Guido W Imbens. The econometrics of randomized experiments. Handbook of Economic Field Experiments, 1:73–140, 2017.
[7] Orazio Attanasio, Britta Augsburg, Ralph De Haas, Emla Fitzsimons, and Heike Harmgart. The impacts of microfinance: Evidence from joint-liability lending in mongolia. American Economic Journal: Applied Economics, 7(1):90–122,
2015.
[8] Britta Augsburg, Ralph De Haas, Heike Harmgart, and Costas Meghir. Microfinance, poverty and education. 2012.
[9] Abhijit Banerjee, Emily Breza, Esther Duflo, and Cynthia Kinnan. Do credit constraints limit entrepreneurship?
heterogeneity in the returns to microfinance. Evanston, USA: Department of Economics Northwestern University, 2015.
[10] Abhijit Banerjee, Esther Duflo, Rachel Glennerster, and Cynthia Kinnan. The miracle of microfinance? evidence
from a randomized evaluation. American Economic Journal: Applied Economics, 7(1):22–53, 2015.
[11] Abhijit Vinayak Banerjee. Microcredit under the microscope: what have we learned in the past two decades, and
what do we need to know? Annu. Rev. Econ., 5(1):487–519, 2013.
[12] G. Barnard. Discussion of “Cross-validatory choice and assessment of statistical predictions” by Stone. Journal of
the Royal Statistical Society. Series B (Methodological), page 133?135, 1974.
[13] A. Belloni, V. Chernozhukov, and C. Hansen. Inference on treatment effects after selection amongst highdimensional controls. Review of Economic Studies, 81:608–650, 2014.
[14] Alexandre Belloni, Victor Chernozhukov, and Kengo Kato. Uniform post selection inference for lad regression
models. arXiv preprint arXiv:1304.0282, 2013.
[15] Alexandre Belloni, Victor Chernozhukov, and Lie Wang. Pivotal estimation of nonparametric functions via squareroot lasso. arXiv preprint arXiv:1105.1475, 2011.
[16] Yoav Benjamini and Yosef Hochberg. Controlling the false discovery rate: a practical and powerful approach to
multiple testing. Journal of the royal statistical society. Series B (Methodological), pages 289–300, 1995.
[17] P. J. Bickel, Y. Ritov, and A. B. Tsybakov. Simultaneous analysis of Lasso and Dantzig selector. Annals of Statistics,
37(4):1705–1732, 2009.
[18] Francine D. Blau and Lawrence M. Kahn. The gender wage gap: Extent, trends, and explanations. Journal of Economic
Literature, 55(3):789–865, September 2017.
[19] V. Chernozhukov, I. Fernández-Val, and A. Galichon. Improving point and interval estimators of monotone functions by rearrangement. Biometrika, 96(3):559–575, 2009.
[20] V. Chernozhukov, I. Fernandez-Val, and Y. Luo. The Sorted Effects Method: Discovering Heterogeneous Effects
Beyond Their Averages. ArXiv e-prints, December 2015.
[21] Victor Chernozhukov, Denis Chetverikov, Mert Demirer, Esther Duflo, Christian Hansen, Whitney Newey, and
James Robins. Double/debiased machine learning for treatment and structural parameters. The Econometrics Journal, 2017.

40

V. CHERNOZHUKOV, M. DEMIRER, E. DUFLO, I. FERNANDEZ-VAL

[22] Victor Chernozhukov, Denis Chetverikov, and Kengo Kato. Anti-concentration and honest, adaptive confidence
bands. The Annals of Statistics, 42(5):1787–1818, 2014.
[23] DR Cox. A note on data-splitting for the evaluation of significance levels. Biometrika, 62(2):441–444, 1975.
[24] Bruno Crépon, Florencia Devoto, Esther Duflo, and William Parienté. Estimating the impact of microcredit on those
who take it up: Evidence from a randomized experiment in morocco. American Economic Journal: Applied Economics,
7(1):123–150, 2015.
[25] Jonathan Davis and Sara B Heller. Rethinking the benefits of youth employment programs: The heterogeneous
effects of summer jobs. Technical report, National Bureau of Economic Research, 2017.
[26] Ruben Dezeure, Peter Bühlmann, and Cun-Hui Zhang. High-dimensional simultaneous inference with the bootstrap. arXiv preprint arXiv:1606.03940, 2016.
[27] Esther Duflo, Rachel Glennerster, and Michael Kremer. Using randomization in development economics research:
A toolkit. Handbook of development economics, 4:3895–3962, 2007.
[28] Jerome Friedman, Trevor Hastie, and Robert Tibshirani. The elements of statistical learning, volume 1. Springer series
in statistics New York, 2001.
[29] Christopher Genovese and Larry Wasserman. Adaptive confidence bands. The Annals of Statistics, pages 875–905,
2008.
[30] Evarist Giné and Richard Nickl. Confidence bands in density estimation. The Annals of Statistics, 38(2):1122–1170,
2010.
[31] Christian Hansen, Damian Kozbur, and Sanjog Misra. Targeted undersmoothing. arXiv preprint arXiv:1706.07328,
2017.
[32] John A Hartigan. Using subsample values as typical values. Journal of the American Statistical Association,
64(328):1303–1317, 1969.
[33] Keisuke Hirano, Guido W. Imbens, and Geert Ridder. Efficient estimation of average treatment effects using the
estimated propensity score. Econometrica, 71(4):1161–1189, 2003.
[34] Kosuke Imai and Marc Ratkovic. Estimating treatment effect heterogeneity in randomized program evaluation. The
Annals of Applied Statistics, 7(1):443–470, 2013.
[35] Guido W Imbens and Donald B Rubin. Causal inference in statistics, social, and biomedical sciences. Cambridge University Press, 2015.
[36] Dean Karlan and Jonathan Zinman. Expanding credit access: Using randomized supply decisions to estimate the
impacts. The Review of Financial Studies, 23(1):433–464, 2009.
[37] Dean Karlan and Jonathan Zinman. Microcredit in theory and practice: Using randomized credit scoring for impact
evaluation. Science, 332(6035):1278–1284, 2011.
[38] Leslie Kish and Martin Richard Frankel. Inference from complex samples. Journal of the Royal Statistical Society. Series
B (Methodological), pages 1–37, 1974.
[39] Max Kuhn. Caret package. Journal of Statistical Software, 28(5):1–26, 2008.
[40] Guillaume Lecué and Charles Mitchell. Oracle inequalities for cross-validation type procedures. Electronic Journal
of Statistics, 6:1803–1837, 2012.
[41] Mark G Low et al. On nonparametric confidence intervals. The Annals of Statistics, 25(6):2547–2554, 1997.
[42] Nicolai Meinshausen, Lukas Meier, and Peter Bühlmann. P-values for high-dimensional regression. Journal of the
American Statistical Association, 104(488):1671–1681, 2009.
[43] Frederick Mosteller and John Wilder Tukey. Data analysis and regression: a second course in statistics. AddisonWesley Series in Behavioral Science: Quantitative Methods, 1977.
[44] Natalia Rigol, Reshmaan Hussam, and Benjamin Roth. Targeting high ability entrepreneurs using community information: Mechanism design in the field, 2016.

GENERIC ML FOR FEATURES OF HETEROGENOUS TREATMENT EFFECTS

41

[45] Alessandro Rinaldo, Larry Wasserman, Max G’Sell, Jing Lei, and Ryan Tibshirani. Bootstrapping and sample splitting for high-dimensional, assumption-free inference. arXiv preprint arXiv:1611.05401, 2016.
[46] Donald B Rubin. Estimating causal effects of treatments in randomized and nonrandomized studies. Journal of
educational Psychology, 66(5):688, 1974.
[47] Charles J. Stone. Optimal global rates of convergence for nonparametric regression. Ann. Statist., 10(4):1040–1053,
1982.
[48] Alessandro Tarozzi, Jaikishan Desai, and Kristin Johnson. The impacts of microcredit: Evidence from ethiopia.
American Economic Journal: Applied Economics, 7(1):54–89, 2015.
[49] Stefan Wager and Susan Athey. Estimation and inference of heterogeneous treatment effects using random forests.
Journal of the American Statistical Association, (just-accepted), 2017.
[50] Larry Wasserman. Machine learning overview. In Becker-Friedman Institute, Conference on ML in Economics, 2016.
[51] Larry Wasserman and Kathryn Roeder. High dimensional variable selection. Annals of statistics, 37(5A):2178, 2009.
[52] Marten Wegkamp et al. Model selection in nonparametric regression. The Annals of Statistics, 31(1):252–273, 2003.
[53] Qingyuan Zhao, Dylan S Small, and Ashkan Ertefaie. Selective inference for effect modification via the lasso. arXiv
preprint arXiv:1705.08020, 2017.

42

V. CHERNOZHUKOV, M. DEMIRER, E. DUFLO, I. FERNANDEZ-VAL

Elastic Net

Random Forest

ATE

90% CB(ATE)

ATE

90% CB(ATE)

GATES

90% CB(GATES)

GATES

90% CB(GATES)

1

2

4000

Treatment Effect

Treatment Effect

4000

2000

2000

0

0

−2000

−2000
1

2

3

4

5

Group by Het Score

3

4

5

Group by Het Score

Figure 3. GATES of Microfinance Availability: Amount of Loans. Point estimates
and 90% adjusted confidence intervals uniform across groups based on 100 random
splits in half
Elastic Net
40000

ATE

90% CB(ATE)

GATES

90% CB(GATES)

Random Forest
40000

90% CB(ATE)

GATES

90% CB(GATES)

1

2

Treatment Effect

20000

Treatment Effect

20000

ATE

0

0

−20000

−20000
1

2

3
Group by Het Score

4

5

3

4

Group by Het Score

Figure 4. GATES of Microfinance Availability: Output. Point estimates and 90%
adjusted confidence intervals uniform across groups based on 100 random splits in
half

5

GENERIC ML FOR FEATURES OF HETEROGENOUS TREATMENT EFFECTS

Elastic Net
ATE

90% CB(ATE)

GATES

90% CB(GATES)

Random Forest

20000

10000

Treatment Effect

Treatment Effect

20000

43

0

−10000

ATE

90% CB(ATE)

GATES

90% CB(GATES)

1

2

10000

0

−10000

1

2

3

4

5

3

Group by Het Score

4

5

Group by Het Score

Figure 5. GATES of Microfinance Availability: Profit. Point estimates and 90% adjusted confidence intervals uniform across groups based on 100 random splits in
half
Elastic Net
ATE

90% CB(ATE)

GATES

90% CB(GATES)

500

Treatment Effect

Treatment Effect

500

Random Forest

0

−500

ATE

90% CB(ATE)

GATES

90% CB(GATES)

1

2

0

−500

1

2

3
Group by Het Score

4

5

3

4

Group by Het Score

Figure 6. GATES of Microfinance Availability: Consumption. Point estimates and
90% adjusted confidence intervals uniform across groups based on 100 random
splits in half

5

44

V. CHERNOZHUKOV, M. DEMIRER, E. DUFLO, I. FERNANDEZ-VAL
Elastic Net

Random Forest

ATE

90% CB(ATE)

ATE

90% CB(ATE)

GATES

90% CB(GATES)

GATES

90% CB(GATES)

−0.1

Treatment Effect

Treatment Effect

−0.1

−0.2

−0.3

−0.2

−0.3

−0.4

−0.4
1

2

3
Group by Het Score

4

5

1

2

3

4

Group by Het Score

Figure 7. GATES of gender wage gap for K = 5 groups defined by the quintiles of
S(Z). Point estimates and 90% adjusted confidence intervals uniform across groups
based on 1,000 random splits in half

5

