Estimation and Inference about Conditional Average Treatment Effect and Other Structural Functions

Estimation and Inference about Conditional Average Treatment Effect and Other Causal Functions

\fnmsVira \snmSemenova\thanksrefm3label=e1]semenovavira@gmail.com [ \fnmsVictor \snmChernozhukov\thanksrefm1label=e1]vchern@mit.edu [ Harvard University\thanksmarkm1, MIT\thanksmarkm3

Abstract

Our framework can be viewed as inference on low-dimensional nonparametric functions in the presence of high-dimensional nuisance function (where dimensionality refers to the number of covariates). Specifically, we consider the setting where we have a signal $Y = Y (η_{0})$ that is an unbiased predictor of causal/structural objects like treatment effect, structural derivative, outcome given a treatment, and others, conditional on a set of very high dimensional controls/ confounders $Z$ . We are interested in simpler lower-dimensional nonparametric summaries of $Y$ , namely $g (x) = E [Y (η_{0}) | X = x],$ conditional on a low-dimensional subset of covariates $X$ . The random variable $Y = Y (η)$ , which we refer to as a signal, depends on a nuisance function $η_{0} (Z)$ of the high-dimensional controls, which is unknown.

In the first stage, we need to learn the function $η_{0} (Z)$ using any machine learning method that is able to approximate $η$ accurately under very high dimensionality of $Z$ . For example, under approximate sparsity with respect to a dictionary, $ℓ_{1}$ -penalized methods can be used; in others, tools such as deep neural networks can be used. To estimate $Y (η_{0})$ we would use $Y (ˆ η)$ , but to make the subsequent inference valid, we need to carefully construct $Y (\cdot)$ such that the signal is orthogonal to perturbations of $η$ , namely $\partial_{η} E [Y (η_{0}) | X] = 0$ . This property allows us to use arbitrary, high-quality machine learning tools to learn $η$ , because it eliminates the first order biases arising from biased/regularized estimation of $η$ necessarily occurring in such cases. As a result, the second-stage low-dimensional nonparametric inference enjoys the quasi-oracle properties, as if we knew $η_{0}$ .

In the second stage, we approximate the target function $g (x)$ by a linear form $p (x)^{'} β_{0}$ , where $p (x)$ is a vector of the approximating functions and $β_{0}$ is the Best Linear Predictor parameter, where the dimension of $p (x)$ is increasing slower than the sample size. We develop a complete set of results about estimation and approximately Gaussian inference on $x \mapsto p (x)^{'} β$ and $x \mapsto g (x)$ . If $p (x)$ is sufficiently rich and $g (x)$ admits a good approximation, then $g (x)$ gets automatically targeted by the inference; otherwise, the best linear approximation $p (x)^{'} β$ to $g (x)$ gets targeted. In some applications, $p (x)$ can be simply specified as a collection of group indicators, in which case $p (x)^{'} β$ describes group-average treatment/structural effects (GATEs), providing a very practical summary of the heterogeneous effects.

[

\kwd

{aug}

class=MSC] \kwd[Primary ]62G08 \kwd[; secondary ]62G05 62G20 62G35

High-dimensional statistics, heterogeneous treatment effect, conditional average treatment effect, orthogonal estimation, machine learning, double robustness

1 Introduction and Motivation.

This paper gives a method to estimate and conduct inference about a nonparametric function $g (x)$ that summarizes heterogeneous treatment/causal/structural effects conditionally on a small subset $X : X \subset Z$ of (potentially very) high-dimensional controls $Z$ . We assume that this function can be represented as a conditional expectation function

g (x)

= E [Y (η_{0}) | X = x],

(1.1)

where the random variable $Y = Y (η)$ , which we refer to as a signal, depends on a nuisance function $η_{0} (Z)$ of the controls. Examples of the nonparametric target function include the Group Average Treatment Effect, Continuous Treatment Effects, the Conditional Average Treatment Effect (CATE), regression function with Partially Missing Outcome, and the Conditional Average Partial Derivative (CAPD). Examples of the nuisance functions

η_{0} = η_{0} (Z)

include the propensity score (probability of treatment assignment), conditional density, and the regression function, among others. In summary,

dim (Z) is high, dim (X) is low.

Although there are multiple possible choices of signals $Y = Y (η)$ for the target function $g (x)$ , we focus on signals that possess the Neyman-type orthogonality property (Neyman (1959)). Formally, we require the pathwise derivative of the conditional expectation to be zero conditionally on $X$ :

\partial_{r} E [Y (η_{0} + r (η - η_{0})) | X = x] |_{r = 0} = 0, for all x and η .

(1.2)

If the signal $Y = Y (η)$ is orthogonal, its plug-in estimate $Y (ˆ η)$ is insensitive to the biased estimation of $ˆ η$ , resulting from application of modern adaptive learning methods in high dimensions, and delivers a high-quality estimator of the target function $g (x)$ under mild conditions.

We demonstrate the importance of the orthogonality property using an example of CATE. We define the CATE function as:

g (x) = E [Y^{1} - Y^{0} | X = x],

where $Y^{1} (Y^{0})$ are the potential outcomes corresponding to receipt(non-receipt) of a binary treatment $D \in {1, 0}$ and $X \in X \subset R^{r}$ is a covariate vector of interest. Consider an Inverse Probability Weighting (IPW, Horwitz and Thompson (1952)) type signal

Y (s) = \frac{D Y^{o}}{s (Z)} - \frac{(1 - D) Y^{o}}{1 - s (Z)},

where $Y^{o} = D Y^{1} + (1 - D) Y^{0}$ is the observed outcome, and the true value of the first-stage nuisance parameter $s (Z)$ is the propensity score $s_{0} (Z) = E [D = 1 | Z]$ . The function $s_{0} (Z)$ is estimated by modern regularized methods, and the bias of its estimation error $ˆ s (Z) - s_{0} (Z)$ converges slowly. The standard choice of signal $Y$ is not orthogonal to the first-stage bias:

\partial_{r} E [Y (s_{0} + r (s - s_{0})) | X] |_{r = 0}

= E [(- \frac{D Y^{o}}{s_{0}^{2} (Z)} - \frac{(1 - D) Y^{o}}{(1 - s_{0} (Z))^{2}}) (s (Z) - s_{0} (Z)) ∣ ∣ ∣ X] \neq 0.

Consequently, the bias of the estimation error, $ˆ s (Z) - s_{0} (Z)$ , in the propensity score translates into the bias of the estimated signal $Y_{ˆ s}$ . As a result, the estimate of the target function based on $Y_{ˆ s}$ is low-quality.

To overcome the transmission of the first-stage bias, we consider the signal of Robins and Rotnitzky (1995) type:

Y (η) = \frac{D (Y^{o} - μ (1, Z))}{s (Z)} - \frac{(1 - D) (Y^{o} - μ (0, Z))}{1 - s (Z)} + μ (1, Z) - μ (0, Z),

where $μ_{0} (D, Z) = E [Y^{o} | D, Z]$ is the regression function, and the true value of the first-stage nuisance parameter $η_{0} (Z) = {s_{0} (Z), μ_{0} (1, Z), μ_{0} (0, Z)}$ . Our proposed choice of signal $Y (η)$ is orthogonal with respect to the first-stage bias

	$\partial_{r} E [Y (η_{0} + r (η - η_{0})) \| X]$
	$= ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ \begin{matrix} E (- \frac{D (Y^{o} - μ_{0} (1, Z))}{s_{0}^{2} (Z)} - \frac{(1 - D) (Y^{o} - μ_{0} (0, Z))}{(1 - s_{0} (Z))^{2}}) (s (Z) - s_{0} (Z)) \| X E (1 - \frac{D}{s_{0} (Z)}) (μ (1, Z) - μ_{0} (1, Z)) \| X - E (1 - \frac{1 - D}{1 - s_{0} (Z)}) (μ (0, Z) - μ_{0} (0, Z)) \| X \end{matrix} ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦$
	$= 0.$

Consequently, the bias of the estimation error, $ˆ η (Z) - η_{0} (Z)$ , in the propensity score does not translate into the bias of the estimated signal $Y (ˆ η)$ . As a result, the estimate of the target function based on $Y (ˆ η)$ is high-quality.

In the second stage we consider a linear projection of an orthogonal signal $Y (η)$ on a vector of basis functions $p (X)$

β := arg min β \in R^{d} E (Y (η) - p (X)^{'} β)^{2} .

The choice of basis functions depends on the problem and the desired interpretation of CATE. For one example, perhaps the simplest choice is to take $X = ⊔_{k = 1}^{d} G_{k}$ be a partition of the support of $X$ into $d$ mutually exclusive groups. Setting

p_{k} (x) = 1_{{x \in G_{k}}}, k \in {1, 2, \dots, d}, x \in X

corresponds to $p (x)^{'} β_{0}$ that can be interpreted as Group Average Treatment Effect. For another example, let $p (x) \in R^{d}$ be a $d$ -dimensional dictionary of series basis functions (e.g., polynomials or splines). Then, $p (x)^{'} β$ corresponds to the best linear approximation to the target function $g (x)$ in the given dictionary. Under some smoothness conditions, as the dimension of the dictionary becomes large, $p (x)^{'} β$ will approximate $g (x)$ at the given point $x$ .

1.1 Overview of Main Results.

The first main contribution of this paper is to provide sufficient conditions for pointwise and uniform asymptotic Gaussian approximation of the target function. We approximate the target function $g (x)$ at a point $x$ by a linear form $p (x)^{'} β_{0}$ :

g (x) = p (x)^{'} β_{0} + r_{g} (x),

where $p (x)$ is a $d$ -vector of technical transformations of the covariates $x$ , $r_{g} (x)$ is the misspecification error due to the linear approximation¹, and $β_{0}$ is the Best Linear Predictor, defined by the normal equation:

E p (X) [g (X) - p (X)^{'} β_{0}] = E p (X) r_{g} (X) = 0.

(1.3)

The two-stage estimator $ˆ β$ of $β_{0}$ , which we refer to as Orthogonal Estimator, is constructed as follows. In the first stage we construct an estimate $ˆ η$ of the nuisance parameter $η_{0}$ , using any high-quality machine learning estimator. In the second stage we construct an estimate ${ˆ Y}_{i}$ of the signal $Y_{i}$ as ${ˆ Y}_{i} := Y_{i} (ˆ η)$ and run ordinary least squares of ${ˆ Y}_{i}$ on the technical regressors $p (X_{i})$ . We use different samples for the estimation of $η$ in the first stage and the estimation of $β_{0}$ in the second stage in a form of cross-fitting, described in the following definition.

Definition 1.1 (Cross-fitting).

For a random sample of size $N$ , denote a $K$ -fold random partition of the sample indices $[N] = {1, 2, . . ., N}$ by $(J_{k})_{k = 1}^{K}$ , where $K$ is the number of partitions and sample size of each fold is $n = N / K$ . Also for each $k \in [K] = {1, 2, . . ., K}$ define $J_{k}^{c} = {1, 2, . . ., N} ∖ J_{k}$ .
For each $k \in [K]$ , construct an estimator ${ˆ η}_{k} = ˆ η (V_{i \in J_{k}^{c}})$ ² of the nuisance parameter value $η_{0}$ using only the data from $J_{k}^{c}$ ; and for any observation $i \in J_{k}$ , define the estimate ${ˆ Y}_{i} := Y_{i} ({ˆ η}_{k})$ .

Definition 1.2 (Orthogonal Estimator).

Given $({ˆ Y}_{i})_{i = 1}^{N}$ , define Orthogonal Estimator as:

ˆ β

:= {(\frac{1}{N} N \sum i = 1 p (X_{i}) p (X_{i})^{'})}^{- 1} \frac{1}{N} N \sum i = 1 p (X_{i}) {ˆ Y}_{i} .

(1.4)

Under mild conditions on $η$ , Orthogonal Estimator delivers a high-quality estimate $p (x)^{'} ˆ β$ of the pseudo-target function $p (x)^{'} β_{0}$ with the following properties:

W.p. $\to 1$ , the mean squared error of $p (x)^{'} ˆ β$ is bounded by

$(E_{N} (p (X_{i})^{'} (ˆ β - β_{0}))^{2})^{1 / 2} = O_{P} (\sqrt{\frac{d}{N}}) + z_{d},$

where $z_{d}$ is the effect of the misspecification error $r_{g} (x)$ .
The estimator $p (x)^{'} ˆ β$ of the pseudo-target function $p (x)^{'} β_{0}$ is asymptotically linear:

$\sqrt{N} \frac{p (x)^{'} (ˆ β - β_{0})}{\sqrt{p (x)^{'} Ω p (x)}} = G_{N} (x) + o_{P} (1 / \sqrt{log N}),$

where the empirical process $G_{N} (x)$ is approximated by a Gaussian process with $N (0, 1)$ coordinates, uniformly over $x \in X$ , and the covariance matrix $Ω$ can be consistently estimated by a sample analog $ˆ Ω$ .
If the misspecification error $r_{g} (x)$ is small, the pseudo-target function $p (x)^{'} β_{0}$ can be replaced by the target function $g (x)$ :

$\sqrt{N} \frac{p (x)^{'} ˆ β - g (x)}{\sqrt{p (x)^{'} Ω p (x)}} = G_{N} (x) + o_{P} (1 / \sqrt{log N}) .$
The quantiles of the suprema of the approximating Gaussian process can be consistently approximated by bootstrap, using which valid uniform confidence bands for $x \mapsto p (x)^{'} β$ and $x \mapsto g (x)$ can be constructed.

The results of this paper accommodate the estimation of $ˆ η$ by high-dimensional/highly complex modern machine learning (ML) methods, such as random forests, neural networks, and $ℓ_{1}$ -shrinkage estimators, as well as estimation of $ˆ η$ by previously developed methods. The only requirement we impose on the estimation of $ˆ η$ is its mean square convergence to the true nuisance parameter $η_{0}$ at a high-quality rate $o_{P} (N^{- 1 / 4 - δ})$ for some $δ ⩾ 0$ . This requirement is satisfied under structured assumptions on $η_{0}$ , such as approximate sparsity of $η_{0}$ with respect to some dictionary, well-approximability of $η_{0}$ by trees or by sparse neural and deep neural nets.

Another important example, worth highlighting here, is that of Continuous Treatment Effects, studied in Kennedy et al. (2017), where

g (x) = E [Y^{x}],

where $X \in X \subset R$ is a continuous treatment and ${Y^{x}, x \in X}$ are the potential outcomes. We prove that the doubly robust signal of Kennedy et al. (2017)

Y (η)

:= \frac{Y^{o} - μ (X, Z)}{s (X | Z)} w (X) + \int μ (X, Z) d P (Z),

is conditionally orthogonal with respect to the nuisance parameter $η_{0} (X, Z)$ consisting of the conditional treatment density $s_{0} (X | Z)$ , regression function $μ_{0} (X, Z)$ , and marginal density $w_{0} (X) = \int s_{0} (X | Z) d P (Z)$ :

η_{0} (X, Z) := {s_{0} (X | Z), μ_{0} (X, Z), w_{0} (X)} .

The proposed estimator of the signal takes the form

Y_{i}^{†} (ˆ η) = \frac{Y_{i}^{o} - ˆ μ (X_{i}, Z_{i})}{ˆ s (X_{i} | Z_{i})} ˆ w (X_{i}) + \frac{1}{n - 1} \sum j \in J_{k}, j \neq i ˆ μ (X_{i}, Z_{j}), i \in J_{k},

where $ˆ η (X, Z)$ is estimated on $J_{k}^{c}$ and the sample average in the second summand is taken over data $(V_{j})_{j \in J_{k}}$ excluding observation $i$ . Having replaced $Y_{i} (ˆ η)$ by $Y_{i}^{†} (ˆ η)$ in Definition 1.2, we obtain an asymptotically linear estimator $p (x)^{'} {ˆ β}^{†}$ of the pseudo-target function:

\frac{\sqrt{N} p (x)^{'} ({ˆ β}^{†} - β_{0})}{} \sqrt{p (x)^{'} Ω^{†} p (x)} = G_{N} (x) + o_{P} (1 / \sqrt{log N}),

(1.5)

where empirical process $G_{N} (x)$ is approximated by a Gaussian process with the same covariance function uniformly over $x \in X$ .

To the best of our knowledge, the inferential approach via series approximation that we propose for this example and all other is novel and so are the inferential results.

1.2 Literature Review.

This paper builds on the three bodies of research within the semiparametric literature: orthogonal(debiased) machine learning, least squares series estimation, and treatment effects/missing data problems. Orthogonal machine learning (Chernozhukov et al. (2016a), Chernozhukov et al. (2016b)) concerns with the inference on a fixed-dimensional target parameter $β$ in presence of a high-dimensional nuisance function $η$ in a semiparametric moment problem. In case the moment condition is orthogonal to the perturbations of $η$ , the estimation of $η$ by ML methods has no first-order effect on the asymptotic distribution of the target parameter $β$ . In particular, plugging in an estimate of $η$ obtained on a separate sample, results in a $\sqrt{N}$ -consistent asymptotically normal estimate whose asymptotic variance is the same as if $η = η_{0}$ was known. This result allows one to use highly complex machine learning methods to estimate the nuisance function $η$ , such as $ℓ_{1}$ penalized methods in sparse models (Bühlmann,Peter and van der Geer,Sara (2011), Belloni et al. (2016)), $L_{2}$ boosting in sparse linear models (Luo,Ye and Spindler,Martin (2016)), and other methods for classes of neural nets, regression trees, and random forests. In many cases, the orthogonal moment conditions are also doubly robust (Robins and Rotnitzky (1995),Kennedy et al. (2017)) with respect to misspecification of one of the nuisance components.

The second building block of our paper is the literature on least squares series estimation (Newey (2007), Newey (2009), Belloni et al. (2015), Chen and Christensen (2015)), which establishes the pointwise and the uniform limit theory for least squares series estimation. We extend this theory by allowing the dependent variable of the series projection to depend upon an unknown nuisance parameter $η$ . We show that series properties continue to hold without any additional strong assumptions on the problem design.

Finally, we also contribute to the literature on missing data, estimation of the conditional average treatment effect and group average treatment effect (Robins and Rotnitzky (1995), Hahn (1998), Graham (2011), Graham et al. (2012), Hirano et al. (2003), Abrevaya and Lieli (2015), Athey and Imbens (2015), Kennedy et al. (2017), Grimmer et al. (2017)). It is also worth mentioning several references that appeared well after our paper was posted to ArXiv (1702.06240). They include Fan et al. (2019), Colangelo and Lee (2019), Jacob et al. (2019), Oprescu et al. (2019), Zimmert and Lechner (2019), and mainly analyze inference on (average, group, or local average) features of CATE(Fan et al. (2019)) and CTE (Colangelo and Lee (2019)). Our framework covers many more examples than features CATE (which was our original motivation), and uses series estimators instead of kernels for localization (e.g., in contrast to Fan et al. (2019), Colangelo and Lee (2019)).

2 Examples.

Examples below apply the proposed framework to study Continuous Treatment Effects, Conditional Average Treatment Effect, regression function with Partially Missing Outcome and Conditional Average Partial Derivative.

Here we introduce a new target function - average potential outcome for the case of continuous treatment $D$ - and derive an orthogonal signal $Y (η)$ for it.

Example 1 (Continuous Treatment Effects).

Let $X = D \in R$ be a continuous treatment variable, $Z$ be a vector of the controls, and ${Y^{d}, d \in D}$ stand for the potential outcomes corresponding to the subject’s response after receiving $d$ units of treatment. The observed data $V = (D, Z, Y^{o})$ consist of the treatment $D$ , controls $Z$ , and the treatment response $Y^{o} = Y^{D}$ . For a given value $x = d$ , the average potential outcome as defined in Gill and Robins (2001) and Kennedy et al. (2017) is

g (x) = E [Y^{x}] = E [Y^{d}] .

(2.1)

Let us define the following nuisance functions. Define the conditional density of $X$ given $Z$ as

s_{0} (x | z) = \frac{d P (X ⩽ t | Z = z)}{d t} {∣ ∣ ∣}_{t = x}

and the conditional expectation of the outcome $Y^{o}$ given $X, Z$ as $μ (x, z) = E [Y^{o} | X = x, Z = z]$ . Finally, let

w_{0} (x) = \int_{z \in Z} s_{0} (x | z) d P_{Z} (z)

be the marginal density of the vector $X$ .The orthogonal signal $Y (η)$ conditionally on $X = x$ is given by

Y (η)

:= \frac{Y^{o} - μ (X, Z)}{s (X | Z)} w (X) + \int μ (X, Z) d P (Z),

(2.2)

where the nuisance parameter $η_{0}$ is a vector-valued function:

η_{0} (X, Z) = {s_{0} (X | Z), w_{0} (X), μ_{0} (X, Z)} .

In particular, (2.2) coincides with the doubly robust signal of Kennedy et al. (2017). Corollary 5.1 shows that Equation (2.2) is orthogonal with respect to the nuisance parameter $η_{0} (D, Z)$ . Theorem 5.1 gives pointwise and uniform asymptotic theory for $g (x)$ . This theorem presents a novel asymptotic linearity representation of $\sqrt{N} α^{'} (ˆ β - β_{0})$ that accounts for the estimation error $\int μ (X, Z) d (P_{N} - P) (Z)$ of the integral $\int μ (X, Z) d P (Z)$ .

Example 2 (Conditional Average Treatment Effect).

Let $Y^{1}$ and $Y^{0}$ be the potential outcomes, corresponding to the response of a subject with and without receiving a treatment, respectively. Let $D \in {1, 0}$ indicate the subject’s presence in the treatment group. The object of interest is the Conditional Average Treatment Effect

g (x) := E [Y^{1} - Y^{0} | X = x] .

Since an individual cannot be treated and non-treated at the same time, only the actual outcome $Y^{o} = D Y^{1} + (1 - D) Y^{0}$ , but not the treatment effect $Y^{1} - Y^{0}$ , is observed.

A standard way to make progress in this problem is to assume unconfoundedness (Rosenbaum and Rubin (1983)). Suppose there exists an observable control vector $Z$ such that the treatment status $D$ is independent of the potential outcomes $Y^{1}, Y^{0}$ conditionally on $Z$ , and $X \subseteq Z$ .

Assumption 2.1 (Unconfoundedness).

The treatment status $D$ is independent of the potential outcomes $Y^{1}, Y^{0}$ conditionally on $Z$ : ${Y^{1}, Y^{0}} ⊥ D | Z$ .

Define the conditional probability of treatment receipt as $s_{0} (Z) = E [D = 1 | Z] .$ Consider an orthogonal signal $Y$ of Robins and Rotnitzky (1995) type:

Y (η)

:= μ (1, Z) - μ (0, Z) + \frac{D [Y^{o} - μ (1, Z)]}{s (Z)} - \frac{(1 - D) [Y^{o} - μ (0, Z)]}{1 - s (Z)},

(2.3)

where $μ_{0} (D, Z) = E [Y^{o} | D, Z]$ is the conditional expectation function of $Y^{o}$ given $D, Z$ . Corollary 5.2 shows that (2.3) is orthogonal to the estimation error of the nuisance parameter $η_{0} (Z) := (s_{0} (Z), μ_{0} (1, Z), μ_{0} (0, Z))$ .

Example 3 (Regression Function with Partially Missing Outcome).

Suppose a researcher is interested in the conditional expectation given a covariate vector $X$ of a variable $Y^{*}$

g (x) := E [Y^{*} | X = x]

that is partially missing. Let $D \in {1, 0}$ indicate whether the outcome $Y^{*}$ is observed, and $Y^{o} = D Y^{*}$ be the observed outcome. Since the researcher does not control the presence status $D$ , a standard way to make progress is to assume existence of an observable control vector $Z$ such that $Y^{*} ⊥ D | Z$ .

Assumption 2.2 (Missingnesss at Random).

The presence indicator $D$ is independent from the outcome $Y$ conditionally on $Z$ : $Y^{*} ⊥ D | Z$ .

Corollary 5.3 shows that the signal $Y$ , defined as

Y (η)

:= μ (Z) + \frac{D [Y^{o} - μ (Z)]}{s (Z)},

(2.4)

is orthogonal to the estimation error of the nuisance parameter $η (Z) := (s (Z), μ (1, Z), μ (0, Z))$ . Here, the function $μ_{0} (Z) = E [Y^{o} | Z] = E [Y | Z, D = 1]$ is the conditional expectation function of the observed outcome $Y^{o}$ given $Z$ and $D = 1$ .

Example 4 (Experiment with Partially Missing Outcome).

Let $Y^{1}$ and $Y^{0}$ be the potential outcomes, corresponding to the response of a subject with and without receiving a treatment, respectively. Suppose a researcher is interested in the Conditional Average Treatment Effect

g (x) := E [Y^{1} - Y^{0} | X = x] .

Conditionally on a vector of stratifying variables $Z_{D}$ , he randomly assigns the treatment status $T \in {1, 0}$ to measure the outcome $T Y^{1} + (1 - T) Y^{0}$ . In presence of the partially missing outcome, let $D \in {1, 0}$ indicate the presence of the outcome record $Y^{o} = D (T Y^{1} + (1 - T) Y^{0})$ in the data. Since the presence indicator $D$ may be co-determined with the covariates $X$ , estimating the treatment effect function $g (x)$ on the observed outcomes only without accounting for the missingnesss may lead to an inconsistent estimate of $g (x)$ .

Since the researcher does not control presence status $D$ , a standard way to make progress is to assume Missingnesss at Random, namely existence of an observable control vector $Z$ such that $Y ⊥ D | Z$ . Setting $Z$ to be a full vector of observables (in particular, $X \subseteq Z$ ) makes Missingness at Random (Assumption 2.2) the least restrictive. Define the conditional probability of presence

s_{0} (Z, T) = E [D = 1 | Y^{o}, Z, T] = E [D = 1 | Z, T]

and the treatment propensity score

h_{0} (Z) := E [T = 1 | Z] .

An orthogonal signal $Y$ for the CATE $g (x)$ can be obtained as follows:

Y (η)

= μ (1, Z) - μ (0, Z) + \frac{D T [Y^{o} - μ (1, Z)]}{s (Z, T) h (Z)} - \frac{D (1 - T) [Y^{o} - μ (0, Z)]}{s (Z, T) (1 - h (Z))},

(2.5)

where $μ_{0} (T, Z) = E [Y^{o} | T, Z]$ is the conditional expectation function of $Y^{o}$ given $T, Z$ .

Example 5 (Conditional Average Partial Derivative).

Let

μ (x, w) := E [Y^{o} | X = x, W = w]

be a conditional expectation function of an outcome $Y^{o}$ given a set of variables $X, W$ . Suppose a researcher is interested in the conditional average derivative of $μ (x, w)$ with respect to $w$ given $X = x$ , denoted by

g (x) = E [\partial_{w} μ (X, w) |_{w = W} | X = x] .

An orthogonal signal $Y (η)$ of Newey and Stoker (1993) type takes the form

Y (η) := - \partial_{w} log f (W | X) [Y^{o} - μ (X, W)] + \partial_{w} μ (X, W),

(2.6)

where $f (W | X = x)$ is the conditional density of $W$ conditionally on $X = x$ . The nuisance parameter $η = {μ (X, W), f (W | X)}$ consists of the conditional expectation function $μ (X, W)$ and the conditional density $f (W | X)$ . Corollary 5.4 shows that the $Y$ is orthogonal to the estimation error of the nuisance parameter $η (X, W) = {μ (X, W), f (W | X)}$ .

3 Empirical Application

To show the immediate usefulness of the method, we consider an important problem of inference about structural derivatives. We apply our methods to study the household demand for gasoline, a question studied in Hausman and Newey (1995), Schmalensee and Stoker (1999), Yatchew and No (2001) and Blundell et al. (2012). These papers estimated the demand function and the average price elasticity for various demographic groups. The dependence of the price elasticity on the household income was highlighted in Blundell et al. (2012), who have estimated the elasticity by low, middle, and high-income groups and found its relationship with income to be non-monotonic. To gain more insight into this question, we estimate the average price elasticity as a function of income and provide simultaneous confidence bands for it.

The data for our analysis are the same as in Yatchew and No (2001), coming from National Private Vehicle Use Survey, conducted by Statistics Canada between October 1994 and September 1996. The data set is based on fuel purchase diaries and contains detailed information about fuel prices, fuel consumption patterns, vehicles and demographic characteristics. We employ the same selection procedure as in Yatchew and No (2001) and Belloni et al. (2011), focusing on a sample of the households with non-zero licensed drivers, vehicles, and distance driven which leaves us with 5001 observations.

The object of interest is the average predicted percentage change in the demand due to a unit percentage change in the price, holding the observed demographic characteristics fixed, conditional on income. In context of Example 5, this corresponds to the conditional average derivative

	$g (x)$	$= E [\partial_{w} μ (X, Z, W) \| X = x],$
	$μ (w, x, z)$	$= E [Y^{o} \| X = x, Z = z, W = w],$

where $Y^{o}$ is the logarithm of gas consumption, $W$ is the logarithm of price per liter, $X$ is log income, and $Z$ are the observed subject characteristics such as household size and composition, distance driven, and the type of fuel usage. The orthogonal signal $Y$ for the target function $g (x)$ is given by

Y = - \partial_{w} log f (W | X, Z) (Y^{o} - μ (X, Z, W)) + \partial_{w} μ (X, Z, W),

(3.1)

where $f (w | x, z) = f (W = w | X = x, Z = z)$ is the conditional density of the price variable $W$ given income $X$ and subject characteristics $Z$ . The conditional density $f (w | x, z)$ and the conditional expectation functions $μ (w, x, z)$ comprise the set of the nuisance parameters to be estimated in the first stage.

The choice of the estimators in the first and the second stages is as follows. To estimate the conditional expectation function $μ (w, x, z)$ and its partial derivative $\partial_{w} μ (w, x, z)$ , we consider a linear model that includes price, price squared, income, income squared, their interactions with 28 time, geographical, and household composition dummies. All in all, we have 91 explanatory variables. We estimate $μ (w, x, z)$ using Lasso with the penalty level chosen as in Belloni et al. (2014), and estimate the derivative $\partial_{w} μ (w, x, z)$ using the estimated coefficients of $μ (w, x, z)$ . To estimate the conditional density $f (w | x, z)$ , we consider a model:

W = l (X, Z) + U, U ⊥ X, Z,

where $l (x, z) = E [W | X = x, Z = z]$ is the conditional expectation of price variable $W$ given income variable $X$ and covariates $Z$ , and $U$ is an independent continuously distributed shock with univariate density $ϕ (\cdot)$ . Under this assumption, the log density $\partial_{w} log f (W | X, Z)$ equals to

\partial_{w} log f (W = w | X = x, Z = z) = \frac{ϕ^{'} (w - l (x, z))}{ϕ (w - l (x, z))} .

We estimate $ϕ (u) : R \to R^{+}$ by an adaptive kernel density estimator of Portnoy and Koenker (1989) with Silverman choice of bandwidth. Finally, we plug in the estimates of $μ (w, x, z)$ , $\partial_{w} μ (w, x, z)$ , $f (w | z, x)$ into the Equation 3.1 to get an estimate of $ˆ Y$ and estimate $g (x)$ by least squares series regression of $ˆ Y$ on $X$ . We try both polynomial basis function and B-splines to construct technical regressors.

Figures 1 and 3 report the estimate of the target function (the black line), the pointwise (the dashed blue lines) and the uniform confidence (the solid blue lines) bands for the average price elasticity conditional on income, where the significance level $α = 0.05$ . The panels of Figure 1 correspond to different choices of the first-stage estimates of the nuisance functions $μ (w, x, z)$ and $f (w | x, z)$ and dictionaries of technical regressors. The panels of Figure 3 correspond to the subsamples of large and small households and to different choices of the dictionaries.

The summary of our empirical findings based on Figure 1 and 3 is as follows. We find the elasticity to be in the range $(- 1, 0)$ and significant for majority of income levels. The estimates based on $B$ -splines (Figures 0(c), 0(d)) are monotonically increasing in income, which is intuitive. The estimates based on polynomial functions are non-monotonic in income. For every algorithm on Figure 1 we cannot reject the null hypothesis of constant price elasticity for all income levels: for each estimation procedure, the uniform confidence bands contain the constant function. Figure 3 shows the average price elasticity conditional on income for small and large households.³. For majority of income levels, we find large households to be more price elastic than the small ones, but the difference is not significant at any income level.

To demonstrate the relevance of demographic data $Z$ in the first stage estimation, we have shown the average predicted effect of the price change on the gasoline consumption (in logs), without accounting for the covariates in the first stage. In particular, this effect equals to $E [\partial_{w} μ (X, W) | X = x]$ , where $μ (x, w) = E [Y | X = x, W = w]$ is the conditional expectation of gas consumption given income and price. This predictive effect consists two effects: the effect of price change on the consumption holding the demographic covariates fixed, which we refer to as average price elasticity, and the association of the price change with the change in the household characteristics that also affect the consumption themselves. Figure 2 shows this predictive effect, approximated by the polynomials of degree $k \in {1, 2}$ , conditional on income. By contrast to the results in Figure 1, the slope of the polynomial of degree $k = 1$ has a negative relationship between income and price elasticity, which present evidence that the demographics $Z$ confound the relationship between income and price elasticity.

(a) Step 2: $l (x, z)$ is estimated by Lasso. Step 3: polynomials of degree $3$ .

(a) Large Households, Polynomials of degree $3$ .

4 Main Theoretical Results.

Notation. For two sequences of random variables $a_{N}, b_{N}, N ⩾ 1 : a_{N} ≲_{P} b_{n}$ means $a_{N} = O_{P} (b_{N})$ . For two sequences of numbers $a_{N}, b_{N}, N ⩾ 1$ , $a_{N} ≲ b_{N}$ means $a_{N} = O (b_{N})$ . Let $a \land b = min {a, b}, a \lor b = max {a, b}$ . The $ℓ_{2}$ norm of a vector is denoted by $∥ \cdot ∥$ , the $ℓ_{1}$ norm is denoted by $∥ \cdot ∥_{1}$ , the $ℓ_{\infty}$ is denoted by $∥ \cdot ∥_{\infty}$ , and the $ℓ_{0}$ - norm denotes the number of non-zero components of a vector. Given a vector $δ \in R^{p}$ and a set of indices $T \subset {1, . . ., p}$ , we denote by $δ_{T}$ the vector in $R^{p}$ for which $δ_{T j} = δ_{j}, j \in T$ and $δ_{T j} = 0, j \notin T$ . Let $ξ_{d} := {sup}_{x \in X} ∥ p (x) ∥ = {sup}_{x \in X} (\sum_{j = 1}^{d} p_{j} (x)^{2})^{1 / 2}$ . For a matrix $Q$ , let $∥ Q ∥$ be the maximal eigenvalue of $Q$ . For a random variable $V$ , let $∥ V ∥_{P, q} := (\int | V |^{q} d P)^{1 / q}$ . Let $ζ_{d} := {sup}_{x \in X} ∥ p (x) ∥_{\infty}$ .

The random sample $(V_{i})_{i = 1}^{N}$ is a sequence of independent copies of a random element $V$ taking values in a measurable space $(V, A_{V})$ according to a probability law $P$ . We shall use empirical process notation. For a generic function $f$ and a generic sample $(X_{i})_{i = 1}^{N}$ , denote a sample average by

E_{N} f (X_{i}) := \frac{1}{N} N \sum i = 1 f (X_{i})

and a $\sqrt{N}$ -scaled, demeaned sample average by

G_{N} f (X_{i}) := \frac{1}{\sqrt{N}} N \sum i = 1 (f (X_{i}) - E f (X_{i})) .

All asymptotic statements below are with respect to $N \to \infty$ . Let

For an observation index $i \in J_{k}$ that belongs to a fold $J_{k}, k \in {1, 2, . ., K}$ , define $Y_{i} (ˆ η) = Y_{i} ({ˆ η}_{k}), i \in J_{k}$ , where ${ˆ η}_{k}$ is estimated on $J_{k}^{c}$ as in Definition 1.1. Define $ˆ η (Z_{i}) = {ˆ η}_{k} (Z_{i}), i \in J_{k}$ .

4.1 Asymptotic Theory for Least Squares Series

Assumption 4.1 (Identification).

Let $Q := E p (X) p (X)^{'} = Q_{d}$ denote population covariance matrix of technical regressors. Assume that $\exists 0 < C_{min} < C_{max} < \infty$ that do not depend on $d$ s.t. $C_{min} < min eig (Q) < max eig (Q) < C_{max} \forall d$ .

Assumption 4.2 (Growth Condition).

We assume that the $sup$ -norm of the technical regressors $ξ_{d} := {sup}_{x \in X} ∥ p (x) ∥ = {sup}_{x \in X} (\sum_{j = 1}^{d} p_{j} (x)^{2})^{1 / 2}$ grows sufficiently slow:

\sqrt{\frac{ξ_{d}^{2} log N}{N}} = o (1) .

Assumption 4.3 (Misspecification Error).

There exists a sequence of finite constants $l_{d}, r_{d}$ such that the norms of the misspecification error are controlled as follows:

∥ r_{g} ∥_{P, 2} := \sqrt{\int r_{g} (x)^{2} d P (x)} ≲ r_{d} and ∥ r_{g} ∥_{P, \infty} := sup x \in X | r_{g} (x) | ≲ l_{d} r_{d} .

Assumption 4.3 introduces the rate of decay of the misspecification error. Specifically, the sequence of constants $r_{d}$ bounds the mean squared misspecification error. In addition, the sequence $l_{d} r_{d}$ bounds the worst-case misspecification error uniformly over the domain of $X$ $X$ , where $l_{d}$ is the modulus of continuity of the worst-case error with respect to mean squared error.

Define the sampling error $U$ as follows:

U := Y - g (X) .

Assumption 4.4 (Sampling Error).

The second moment of the sampling error $U$ conditionally on $X$ is bounded from above by ${¯ σ}^{2}$ :

sup x \in X E [U^{2} | X = x] ≲_{P} {¯ σ}^{2} .

To describe the first-stage rate requirement, Assumption 4.5 introduces a sequence of nuisance realization sets $T_{N}$ for the nuisance parameter $η_{0}$ . As sample size $N$ increases, the sets $T_{N}$ shrink around the true value $η_{0}$ . The shrinkage speed is described in terms of the statistical rates $B_{N}$ and $Λ_{N}$ .

Assumption 4.5 (Small Bias Condition).

There exists a sequence $ϵ_{N} = o (1)$ , such that with probability at least $1 - ϵ_{N}$ , the first stage estimate $ˆ η$ , obtained by cross-fitting (Definition 1.1), belongs to a shrinking neighborhood of $η_{0}$ , denoted by $T_{N}$ . Uniformly over $T_{N}$ , the following mean square convergence holds:

	$B_{N} := \sqrt{N} sup η \in T_{N} ∥ E p (X) [Y (η) - Y (η_{0})]$	$∥ = o (1),$		(4.1)
	$Λ_{N} := sup η \in T_{N} (E ∥ p (X) [Y (η) - Y (η_{0})] ∥^{2})^{1 / 2}$	$= o (1) .$		(4.2)

In particular, $Λ_{N}$ can be bounded as

Λ_{N} ≲ ξ_{d} sup η \in T_{N} (E (Y (η) - Y (η_{0}))^{2})^{1 / 2} .

Assumption 4.5 is stated in a high-level form in order to accommodate various machine learning estimators. We demonstrate the plausibility of Assumption 4.5 for a high-dimensional sparse model for Example 3 (see Example 6).

Pointwise Limit Theory

Theorem 4.1 (Pointwise Limit Theory of Orthogonal Estimator).

Let Assumptions 4.1, 4.2, 4.3, 4.4, 4.5 hold. Then, the following statements hold:

The second norm of the estimation error is bounded as:

$∥ ˆ β - β_{0} ∥_{2} ≲_{P} \sqrt{\frac{d}{N}} + [\sqrt{\frac{d}{N}} l_{d} r_{d} \land ξ_{d} r_{d} / \sqrt{N}],$

which implies a bound on the mean squared error of the estimate $p (x)^{'} ˆ β$ of the pseudo-target function $p (x)^{'} β_{0}$ :

$(E_{N} (p (X_{i})^{'} (ˆ β - β_{0}))^{2})^{1 / 2} ≲_{P} \sqrt{\frac{d}{N}} + [\sqrt{\frac{d}{N}} l_{d} r_{d} \land ξ_{d} r_{d} / \sqrt{N}] .$
For any $α \in S^{d - 1} := {α \in R^{d} : ∥ α ∥ = 1}$ the estimator $ˆ β$ is approximately linear:

$\sqrt{N} α^{'} ($