Double Machine Learning Density Estimation for Local Treatment Effects with Instruments
Year: 2022

Abstract:
Local treatment effects are a common quantity found throughout the empirical sciences that measure the treatment effect among those who comply with what they are assigned. Most of the literature is focused on estimating the average of such quantity, which is called the “local average treatment effect (LATE)” [31]). In this work, we study how to estimate the density of the local treatment effect, which is naturally more informative than its average. Specifically, we develop two families of methods for this task, namely, kernel-smoothing and model-based approaches. The kernel-smoothing-based approach estimates the density through some smooth kernel functions. The model-based approach estimates the density by projecting it onto a finite-dimensional density class. For both approaches, we derive the corresponding double/debiased machine learning-based estimators [13]. We further study the asymptotic convergence rates of the estimators and show that they are robust to the biases in nuisance function estimation. The use of the proposed methods is illustrated through both synthetic and a real dataset called 401(k).

heading: 1 Introduction
text: Controlled experimentation is one powerful tool used throughout the empirical sciences to infer the effect of a certain treatment on a given outcome. The idea is to randomize the treatment assignment so as to neutralize the effect of unobserved confounders. However, in some practical settings, it may be challenging to ascertain that individuals who are selected for treatment will follow their recommendations. Issues of non-compliance and unmeasured confounding are quite common and lead to the non-identification of treatment effects in many real-world cases [29, 50, 32, 56].
An approach known as instrumental variables (IVs) has been proposed to try to circumvent this issue [68]. The idea is to find a set of variables (possibly singleton) that are not the target of the analysis by itself but that will help to control for the unobserved confounding between the treatment and the outcome. In particular, IVs are special variables that (i) are correlated with the treatment, (ii) do not directly influence the outcome, and (iii) are not affected by certain unmeasured confounders. For concreteness, consider a study of the effect of 401(k) participation (X) on the distribution of net financial assets (Y ) [2]. This setting is represented in the causal graph in Fig. 1. Note that a dashed-bidirected arrow exists between X and Y , which in graphical language represents unobserved confounding affecting both X and Y . The variable Z in this model represents the eligibility of 401(k). We note that Z qualifies as an instrument in this case – (i) it does affect the participation of 401(k) (X) and (ii) has no direct influence on the net financial asset (Y ), (iii) is not affected by unmeasured confounders between X and Y . The variable W represents observed covariates (e.g., gender, age, ethnicity, income, family size).
We are interested in the particular setting where only individuals who were offered the treatment may have access to it [31]. For instance, in the case of 401(k) participation (X = 1), only eligible
35th Conference on Neural Information Processing Systems (NeurIPS 2021).
individuals (Z = 1) would be allowed to join the program. This assumption is known in the literature as monotonicity, which rules out the possibility that any units would respond contrary to the instrument. Under monotonicity, the causal effect in the subpopulation whose actual treatment X coincides with the assigned treatment Z (called compliers) is identifiable [31, 2]. The average treatment effect (ATE) for the compliers is called ‘Local ATE’ (LATE) (or Complier average causal effects, CACE) [31].
The most common quantification of these effects in IV settings found in practice is the average (e.g., LATE). The average is certainly an informative summary; however, it may fail to capture significant differences in the causal distributions of the outcome. For instance, consider Fig. 2 that shows the densities of outcomes Y under treatmentsX = 1 among compliers which are generated from samples drawn from four synthetic data generating processes represented by the IV graph in Fig. 1 (further discussed in Sec. 5). All of the four distributions have the same mean 0 and variance 2. However, the difference in the LTE distributions is self-evident.
Most of the prior work on quantifying distributions of treatment effects focuses on estimating cumulative distribution functions (CDFs) or quantiles, and little attention has been given to estimating densities (refer to Sec. 1.1 for further comparison). As a complement to CDFs, densities have various advantages, including a more interpretable visualization of the distribution and generative capability of producing samples. One challenge with estimating densities is that while CDFs are pathwise-differentiable and enjoy √ n-rate estimators (n is the size of data), densities are not (i.e., they are non-regular), and therefore possess no influence functions nor √ n-rate estimators without approximations [7, Ch. 3].
In this paper, our goal is to provide methods to estimate densities of local treatment effects in IV settings under the monotonicity assumption. We develop two families of methods for this task based on kernel-smoothing and model-based approximations. The former smooths the density by convolution with a kernel function; the latter projects the density onto a finite-dimensional density class based on a distributional distance measure. For both approaches, we construct double/debiased machine learning (DML) style density estimators [43, 54, 52, 70, 13]. We analyze the asymptotic convergence properties of the estimators, showing that they can converge fast (i.e., √ n-rate) even when nuisance estimates converge slowly (e.g., n−1/4 rate) (a property called ‘debiasedness’1). We illustrate the proposed methods on synthetic and real data.

heading: 1.1 Related work
text: Our work touches different areas, which we will discuss next.
Double/Debiased Machine Learning (DML) [13]-based causal effect estimators. The DML framework has been adapted for estimating the average causal effect under the setting where the back-door criterion [50, Sec. 3.3.1] (also known as ignorability [57]) holds (e.g., [12, 19]). Recently, DML-based causal effect estimators have been developed for any identifiable causal functionals in a given causal graph and equivalence class thereof [33, 34].
1Also known as ‘nonparametric doubly robust [37] or ‘rate doubly robust’ [59].
Local average & quantile treatment effect. The formal identification results for LATE under the monotonicity assumption in IV settings were developed by [31, 3]. Building on these results, semiparametric estimation for LATE has received remarkable attention [2, 60, 23, 62, 48], including robust LATE estimators that achieve debiasedness [47, 40, 38, 64]. As shown in Fig. 2, however, the average is sometimes insufficient to capture the effects of the treatment on the distributions of outcomes. To address this issue, the problem of estimating quantiles or CDFs has taken attention. A common approach to estimate quantiles or CDFs is based on the LATE estimation. Since the expectation of 1Y≤y(Y ), an indicator that outcome Y falls short of threshold y, reduces to the CDF (i.e., replacing Y in LATE with 1Y≤y(Y )), estimators for the LATE can be used to estimate quantiles or CDFs [1, 2, 15, 24, 16, 30, 45, 18, 69]. Non-regular target estimand. An estimand that possesses no influence functions nor √ n-rate estimators is called ‘non-regular’. Densities are an example of non-regular target estimands [7, Chap. 3]. One can approximate a non-regular target with a smooth one such that an influence function and√ n-rate estimators can be derived. Two broadly used approaches are kernel-smoothing-based (e.g., [52, 6, 42, 19, 35]) and model-based (e.g., [46, 52, 21, 41, 40, 39]).
Causal density estimation. There is limited literature on estimating the density of treatment effects. Most of the results assume that the ignorability/backdoor admissibility holds [55, 49]. [22] used the kernel-smoothing technique to estimate the density of a treatment effect, and [42] provided a kernelsmoothing-based density estimator that achieves doubly robustness and debiasedness building on top of the work in [53]. Recently, [39] investigated a model-based approach and developed estimators that achieve debiasedness properties. Under the IV setting, [10] provided a local polynomial regressionbased density estimator for local treatment effects; We are not aware of any work studying debiased density estimators. As mentioned, this paper investigates both kernel-smoothing and model-based approaches for estimating local treatment effects under IV settings and develops DML-style density estimators for both.

heading: 2 LTE Estimation – Problem setup
text: In our analysis, each variable is represented with a capital letter (X) and its realized value with a small letter (x). For a discrete (e.g., binary) random variable X , we use 1x(X) to represent the indicator function such that 1x(X) = 1 if X = x; 1x(X) = 0 otherwise. For a continuous variable X with a probability density p(x) of a distribution P and a function f(x), EP [f(X)] ≡ ∫ X f(x)p(x) d[x]
where X is the domain for X , and ∥f(X)∥ ≡ √ EP [(f(X))2]. f̂ is said to converge to f at rate rn if
∥f̂(x)− f(x)∥ = OP (1/rn). For a dataset D = {Vi}ni=1, we use ED [f(V )] ≡ (1/n) ∑n i=1 f(Vi) to denote the empirical mean of f(V ) with D. Structural Causal Models (SCMs). We use the language of SCMs as our basic semantic and inferential framework [50, 4]. An SCMM is a quadrupleM = ⟨U, V, P (U), F ⟩ where U is a set of exogenous (latent) variables following a joint distribution P (u), and V is a set of endogenous (observable) variables whose values are determined by functions F = {fVi}Vi∈V such that Vi ← fVi(pai, ui) where PAi ⊆ V and Ui ⊆ U . Each SCMM induces a distribution P (v) and a causal graph G = G(M) over V in which there exists a directed edge from every variable in PAi to Vi and dashed-bidirected arrows encode common latent variables (e.g., see Fig. 1). Within the structural semantics, performing an intervention and setting X = x is represented through the do-operator, do(X = x), which encodes the operation of replacing the original equations of X (i.e., fX(pax, ux)) by the constant x and induces a submodelMx and an interventional distribution P (v|do(x)). For any variable Y ∈ V , the potential response Yx(u) is defined as the solution of Y in the submodel Mx given U = u, which induces a counterfactual variable Yx.
Local Treatment Effect (LTE) with IV. We consider the IV setting represented by the causal graph G in Fig. 12, where Z is a binary instrument with domain {0, 1}, X is a binary treatment with domain {0, 1}, and Y is a (set of) continuous outcomes with bounded domain Y ⊂ Rd, and W is a set of covariates (continuous, discrete, or mixed). G satisfies the IV assumption that Z has no direct influence on outcome Y and is not affected by unmeasured confounders between X and Y .
2It is common in the literature to define IV assumptions in terms of conditional independences among counterfactual [51, 9, 8, 2, 60, 47, 64], whose connection with the causal graph in Fig. 1 is discussed in Assumption A.1
The causal density p(yx) is not identifiable from the observed density p(x, y, z, w) due to the unobserved confounders between X and Y . However, the effect is possibly recovered for certain subpopulation under additional assumptions. Formally, a unit in the population is an always-taker if XZ=1 = XZ=0 = 1, a never-taker if XZ=1 = XZ=0 = 0, a complier if XZ=1 = 1, XZ=0 = 0, and a defier if XZ=1 = 0, XZ=0 = 1 [3, 2]. We will make the following assumptions based on literature. Assumption 1 (Monotonicity). There are no defiers: XZ=1 ≥ XZ=0. Assumption 2 (Positivity). P (x|z, w) > 0, P (z|w) > 0 for any x, z, w.
Let C denote the event that a unit is a complier (i.e., a unit such that XZ=0 = 0 and XZ=1 = 1). For a given constant a and a variable X , let xa denote the event X = a. The LTE p(yx|C) is identifiable under monotonicity and is given by [31, 2]:
p(yx|C) = EP [ p(y|x, zx,W )P (x|zx,W )− p(y|x, z1−x,W )P (x|z1−x,W ) ] EP [P (x1|z1,W )− P (x1|z0,W )] , (1)
where the expectation is over W . In this paper, we aim to estimate the LTE density p(yx|C) in Eq. (1). We will make the following mild assumption on some densities, popularly employed in density estimation literature (e.g., [44, 25, 27, 61, 26, 42]). Assumption 3. For any x, z, w, y, densities p(y|w, z, x), p(y|z, x) and p(yx|C) are bounded, and p(yx|C) is twice differentiable.
DML method. Let ψ ≡ ψP ′ denote a functional of an arbitrary distribution P ′. We use P to denote the true distribution such that D ∼ P . Let ψ0 ≡ ψP denote the true parameter to be estimated. To estimate ψ0, DML-based estimators use a Neyman Orthogonal score φ(V ;ψ0, η) (where η ≡ ηP ′ is a set of nuisance parameters and η0 ≡ ηP denotes the true nuisances), a function such that EP [φ(V ;ψ0, η0)] = 0, (∂/∂η)|η=η0EP [φ(V ;ψ0, η)] = 0. Given φ, an DML estimator is constructed using the cross-fitting technique as follows: For randomly split halves of D denoted {D0,D1}, let η̂p for p ∈ {0, 1} denote the estimates for η fromDp. Let Tp denote a solution such that ED1−p [φ(V ;Tp, η̂p)] = oP (N−1/2). Then, T ≡ (T0 + T1)/2 is an DML estimator [13, Def. 3.1]. In addition to being consistent, the estimator T exhibits a robustness property called debiasedness: T converges to ψ0 in the root-N rate even when η̂ converges to η0 in slower N−1/4 rate [13, Thm. 3.1]. A Neyman Orthogonal Score can be derived by adding ψ to its influence function ϕ [14, Thm. 1]. An influence function of the functional ψP is defined as a solution satisfying EP [ϕ] = 0, EP [ ϕ2 ] <∞, and (∂/∂t)ψPt |t=0 = EP [ϕ(V ;ψ, η)St(V ; t = 0)] where Pt ≡ P (v)(1 + tg(v)) for t ∈ R and any bounded mean-zero functions g(·) over V , and St(v; t = 0) ≡ (∂/∂t) logPt(v)|t=0 [63, Chap. 25]. Due to space constraints, all the proofs are provided in Appendix B in suppl. material.

heading: 3 Kernel-smoothing-based approach
text: In this section, we develop a kernel-smoothing-based approach for estimating the LTE density. The kernel-smoothing technique approximates a non-pathwise-differentiable target estimand with a differentiable estimand by convoluting the density with a kernel function K(y). Properties of the kernel function includes symmetry about the origin (i.e., ∫ Y yK(y) d[y] = 0), non-negativity
(0 < K(y) <∞, ∀y ∈ Y) , and integrates to 1 (i.e., ∫ Y K(y) d[y] = 1) [66, Chap. 4.2].
We consider a product kernel Kh,y(y′) ≡ h−d ∏d j=1K((yj − y′j)/h) with given bandwidth h ∈ R and a fixed point y = {yj}dj=1 ∈ Rd. We assume that the kernel of interest has a bounded second moment and norm: i.e., κ2(K) ≡ ∫ Y y
2K(y) d[y] < ∞ and ∥K(y)∥ < ∞ following [27, 61]. Example of kernels include Gaussian kernel: K(u) = (1/ √ 2π) exp ( −u2/2 ) for u ∈ R,
Epanechnikov kernel: K(u) = (3/4)(1 − u2)1|u|≤1(u), Quadratic kernel: K(u) = (15/16)(1 − u2)21|u|≤1(u), Cosine kernel: k(u) = (π/4) cos (πu/2)1|u|≤1(u), etc.
For convenience, we denote the target estimand by ψ(y) ≡ p(yx|C). In the kernel-smoothing-based approach, we will aim to estimate a kernel-smoothed approximation for ψ(y) defined as follows:
ψh(y) ≡ ∫ Y ψ(y′)Kh,y(y ′) d[y′] = ψ[Kh,y(Y )], (2)
where ψ[f(Y )] is an expectation of a function f(Y ) w.r.t. ψ(y), which is specified as ψ[f(Y )] ≡ EP [ EP [f(Y )1x(X)|zx,W ]− EP [ f(Y )1x(X)|z1−x,W ]] EP [P (x1|z1,W )− P (x1|z0,W )] . (3)
The second equality in Eq. (2) is by Eq. (1). For a target estimand ψ[f(Y )], we will denote nuisances by πz(w) ≡ P (z|w), ξx(z, w) ≡ P (x|z, w), and θ(x, z, w)[f(Y )] ≡ EP [f(Y )1x(X)|z, w], shortly (π, ξ, θ).
We aim to construct a DML estimator for the estimand ψh. Toward this goal, we will first derive a Neyman orthogonal score for ψh. Since a Neyman orthogonal score can be constructed based on moment score functions (a function of parameters such that its expectation is 0 at the true parameters) [14, Thm. 1], we start by defining the moment score function. Let
ψX ≡ EP [ ξx1(z 1,W )− ξx1(z0,W ) ] , (4)
VX({π, ξ}) ≡ 1z1(Z)− 1z0(Z)
πZ(W ) {1x1(X)− ξx1(Z,W )}+
{ ξx1(z 1,W )− ξx1(z0,W ) } . (5)
Then, the following is a moment score function for ψh:
m(ψ′;ψh) ≡ 1
ψX (ψh − ψ′)VX , (6)
where ψh is given in Eq. (2) and ψ′ is an estimate of ψh.
Next, we derive an influence function for the moment score function m(ψ′;ψh). We first define the following function: for a bounded function f(Y ) <∞, let
ψY X [f(Y )] ≡ EP [ θ(x, zx,W )[f(Y )]− θ(x, z1−x,W )[f(Y )] ] , (7)
VY X({π, θ})[f(Y )] ≡ 1zx(Z)− 1z1−x(Z)
πZ(W ) {f(Y )1x(X)− θ(x, Z,W )[f(Y )]} + { θ(x, zx,W )[f(Y )]− θ(x, z1−x,W )[f(Y )] } , (8)
and
ϕ(η = {π, ξ, θ}, ψ)[f(Y )] ≡ 1 ψX (VY X({π, θ})[f(Y )]− ψ[f(Y )]VX({π, ξ})) , (9)
where VX is defined in Eq. (5). Then, the influence function for the expectation of the moment score function m(ψ′;ψh) in Eq. (6) is given as follows: Lemma 1 (Influence function for m(ψ′;ψh)). Let m(ψ′;ψh) be the score defined in Eq. (6). Then, the influence function for EP [m(ψ′;ψh)], denoted ϕm, is given by
ϕm(η = {π, ξ, θ}, ψ) ≡ ϕ(η, ψ)[Kh,y(Y )] (10) where ϕ is in Eq. (9).
For any score function (e.g., m in Eq. (6)), its addition to the influence function of the expected score (e.g., ϕm) is a Neyman orthogonal score3 ([14, Thm.1], [13, Sec. 2.2.5]). Specifically, Lemma 2 (Neyman orthogonal score for ψh). Let m(ψ′;ψh) be the score function in Eq. (6), and ϕm(η = {π, ξ, θ}, ψh) be the influence function for EP [m(ψ′;ψh)] given in Eq. (10). Then, a Neyman orthogonal score for ψh is given as φ(ψ′; η = {π, ξ, θ}) ≡ m(ψ′;ψh) + ϕm(η, ψ); Specifically,
φ(ψ′; η = {π, ξ, θ}) = 1 ψX (VY X({π, θ})[Kh,y(Y )]− ψ′VX({π, ξ})) . (11)
Given the Neyman orthogonal score φ(ψ′; η), an estimate ψ̂h satisfying ED [ φ(ψ̂h; η̂ = {π̂, ξ̂, θ̂}) ] = oP (n
−1/2) gives a DML estimator. Specifically, we propose the following kernel-smoothing based estimator for the LTE density, named ‘KLTE’ (kernel-based estimator for LTE):
3A Neyman orthogonal score is a function ϕ satisfying EP [ϕ(ψ, η0)] = 0 and ∂∂ηEP [ϕ(V ;ψ, η)]|η=η0 = 0, where η0 denotes the true nuisance [13, Def.2.2]. In words, a score function that is not sensitive to local errors in nuisance models.
Definition 1 (KLTE estimator for ψh). Let φ(ψ′; η = {π, ξ, θ}) be the Neyman orthogonal score for ψh given in Eq. (11). Let {D,D′} denote the randomly split halves of the samples, where |D| = |D′| = n. Let η̂ = {π̂, ξ̂, θ̂} denote the estimates for the nuisance η using D′. Then, the KLTE estimator for ψh(y) for all y ∈ Y , denoted ψ̂h(y), is given by
ψ̂h(y) ≡ ED [ VY X({π̂, θ̂})[Kh,y(Y )] ] /ED [ VX({π̂, ξ̂}) ] , (12)
where VX and VY X are given in Eqs. (5,8), respectively.
We will show that the KLTE is a DML estimator exhibiting debiasedness property. Detailed asymptotic properties are discussed next.

heading: 3.1 Asymptotic convergence
text: Now, we study the convergence rate of the estimator ψ̂h(y). For any fixed y ∈ Y , the error ψ̂h(y) − ψ(y) will be analyzed in two folds: we will first analyze the error between the estimator in Eq. (12) and the smoothed estimand in Eq. (2) (i.e., ψ̂h(y)− ψh(y)), and then analyze the error between the smoothed estimand and the true estimand (i.e., ψh(y)− ψ(y)).
The following result gives the error analysis for ψ̂h(y)− ψh(y): Lemma 3 (Convergence rate of ψ̂h to ψh). For any fixed y ∈ Y , suppose the estimators for nuisances are consistent; i.e., ∥ν − ν̂∥ = oP (1) for ν ∈ η = {π, ξ, θ} for all (w, z, x). Suppose h <∞, and nhd →∞ as n→∞. Then,
ψ̂h(y)− ψh(y) = OP ( 1/ √ nhd +Rk2 + 1/ √ n ) ,
where
Rk2 ≡ ∑ z ∥π̂z − πz∥ {∥∥∥θ̂z − θz∥∥∥+ ∥∥∥ξ̂z − ξz∥∥∥} , (13)
where πz ≡ πz(W ), ξz ≡ ξx(z,W ) and θz ≡ θ(x, z,W )[Kh,y(Y )].
The error analysis in Lemma. 3 implies the following:
Corollary 1 (Debiasedness property of ψ̂h toψh). If all nuisances {π̂, ξ̂, θ̂} for any given (w, z, x, y) converge at rate {nhd}−1/4, then the target estimator ψ̂h(y) achieves √ nhd-rate convergence to ψh.
We now analyze the gap between the smoothed estimand ψh and the true estimand ψ; i.e., ψh − ψ: Lemma 4 ([66, Thm. 6.28]). The following holds:
ψh(y)− ψ(y) = By ≡ 0.5h2κ2(K)(∂2/∂2y′)|y′=yψ(y′) +O(h2). (14)
Combining the results of Lemma. (3,4), we have the following result:
Theorem 1 (Convergence rate of ψ̂h to ψ). For any fixed y ∈ Y , suppose the estimators for nuisances are consistent; i.e., ∥ν − ν̂∥ = oP (1) for ν ∈ η = {π, ξ, θ} for all (w, z, x). Suppose h <∞, and nhd →∞ as n→∞. Then
ψ̂h(y)− ψ(y) = OP ( 1/ √ nhd +Rk2 + 1/ √ n ) +By, (15)
where By is defined in Eq. (14), and Rk2 is defined in Eq. (13).
Thm. 1 implies that ψ̂h(y) converges fast (see Corol. 1) to ψ(y)+By . A natural question is then how to choose the bandwidth h that minimizes the gap in Eq. (15). The following provides a guideline in choosing the bandwidth h: Lemma 5 (Data-adaptive bandwidth selection). The bandwidth h that minimizes the error in Eq. (15) is h = O(n−1/(d+4)). This choice of h satisfies the assumption in Lemma 3 (i.e., nhd →∞).
Recall that Corol. 1 states the debiasedness property of ψ̂h to ψh for any bandwidth h satisfying nhd →∞. With the choice of h as in Lemma 5, ψ̂h converges to ψ with the debiasedness property preserved.
Corollary 2 (Debiasedness property of ψ̂h to ψ). Let h = O(n−1/(d+4)). If nuisances {π̂, ξ̂, θ̂} converge at {nhd}−1/4 rate for any (w, z, x, y), then the target estimator ψ̂h(y) achieves √ nhd-rate convergence to ψ.
So far, we have analyzed the error ψ̂h(y) − ψ(y) pointwise for the fixed y ∈ Y . To analyze the difference between the two densities ψ̂h(y) and ψ(y) for all y ∈ Y , we consider the following divergence function of two densities: Definition 2 (f -Divergence Df [20]). Let f denote a convex function with f(1) = 0. Df (p, q) ≡∫ Y f(p(y), q(y))q(y) d[y], is a f -divergence function between two densities p, q.
f -divergence covers many well-known divergences. For example, Df reduces to KL divergence with f(p, q) = (p/q) log(p/q). We will assume that the function f(p, q) in Df is differentiable w.r.t. p and q.
We now analyze the distance between ψ̂h and ψ w.r.t. Df . The following result provides an upper bound for Df .
Lemma 6 (Upper bound of the divergence Df ). Suppose Df is a f -divergence such that f(p, q) = 0 if p = q. Then,
Df (ψ, ψ̂h) ≤ ∫ Y w(y) ( ψ̂h(y)− ψ(y) ) d[y],
where w(y) ≡ f ′2(ψ(y), ψ̃(y))ψ̂h(y), f ′2(p, q) ≡ (∂/∂q)f(p, q), and ψ̃h(y) ≡ tψ̂h(y)+ (1− t)ψ(y) for some fixed t ∈ [0, 1].
By invoking Thm. 1, we derive an upper bound for Df (ψ, ψ̂h) as follows:
Theorem 2 (Convergence rate of ψ̂h). Suppose the estimators for nuisances are consistent; i.e., ∥ν − ν̂∥ = oP (1) for ν ∈ η = {π, ξ, θ} for all (w, z, x, y). Suppose Df is a f -divergence such that f(p, q) = 0 if p = q. Suppose w(y) in Lemma 6 is finite. Then,
Df (ψ, ψ̂h) ≤ OP ( sup y∈Y { Rk2 +By } + 1/ √ nhd + 1/ √ n ) , (16)
where Rk2 is defined in Eq. (13) and By is defined in Eq. (14).
The following result asserts that the debiasedness property is exhibited w.r.t. Df :
Corollary 3 (Debiasedness property of ψ̂h w.r.t. Df ). Let h = O(n−1/(d+4)). Suppose Df satisfies f(p, q) = 0 if p = q. Suppose w(y) in Lemma 6 is finite. If nuisances {π̂, ξ̂, θ̂} converges at {nhd}−1/4 rate for any (w, z, x, y), then Df (ψ, ψ̂h) converges to 0 at √ nhd-rate.

heading: 4 Model-based approach
text: In this section, we develop a model-based approach for estimating the LTE density ψ(y) = p(yx|C). We will approximate ψ with a class of distributions or a density model G = {g(y;β) : β ∈ Rb} where g(y;β) ∈ G is differentiable w.r.t. β. Example density models include exponential family (e.g., Gaussian distribution), mixture of Gaussians, or more generally, mixture of exponential families. The choice of the density model may depend on domain knowledge. Alternatively, one may choose among a set of candidate density families using separate validation data or applying cross-validation. We adapt the model-based approach developed in [39] for estimating the causal density under the no unmeasured confounders assumption.
Given a density model G, the best approximation for ψ(y) is defined as g(y;β0) ∈ G that achieves the minimum f -divergence to ψ:
β0 ≡ arg min β∈Rb Df (ψ(y), g(y;β)), (17)
where Df is the f -divergence defined in Def. 2. Our goal is estimating β0.
Consider m(β;ψ) ≡ (∂/∂β)Df (ψ(y), g(y;β)). Definition of β0 given in Eq. (17) implies that m(β;ψ) = 0 at β = β0. We note that m(β;ψ) serves as a moment score function. The closed-form expression of the score is given by [39]:
m(β;ψ) ≡ ∫ Y g′(y;β) {f ′2(ψ(y), g(y;β))g(y;β) + f(ψ(y), g(y;β))} d[y], (18)
where g′(y;β) = (∂/∂β)g(y;β) and f ′2(p, q) ≡ (∂/∂q)f(p, q). To construct a DML estimator based on the score function m(β;ψ), we first derive an influence function for the score: Lemma 7 (Influence Function for m(β, ψ)). An influence function for m(β;ψ) in Eq. (18), denoted ϕm, is given by
ϕm(β; η = {π, ξ, θ}, ψ) ≡ ϕ(η, ψ)[Rf (Y ;β, ψ)], (19) where ϕ(η, ψ)[·] is defined in Eq. (9), and
Rf (Y ;β, ψ) ≡ g′(Y ;β) {f ′′21(ψ(Y ), g(Y ;β))g(Y ;β) + f ′1(ψ(Y ), g(Y ;β))} , where g′(y;β) ≡ (∂/∂β)g(y;β), f ′1(p, q) ≡ (∂/∂p)f(p, q) and f ′′21(p, q) ≡ (∂/∂p)f ′2(p, q).
We derive a Neyman orthogonal score based on the moment score m(β, ψ) and its influence function ϕm(β, η, ψ): Lemma 8 (Neyman orthogonal score for β). A Neyman orthogonal score for estimating β, denoted φ(β′; (η = {π, ξ, θ}, ψ)), is given by
φ(β′; (η = {π, ξ, θ}, ψ)) ≡ m(β′, ψ) + ϕm(β, η, ψ), (20) where ϕm(β, η, ψ) is defined in Eq. (19).
Given the orthogonal score φ(β′; (η, ψ)) in Eq. (20), we propose the following estimator for β, named ‘MLTE’ (model-based estimator for LTE): Definition 3 (MLTE estimator for β). Let φ(β′; η = {π, ξ, θ}, ψ) be the Neyman orthogonal score for β given in Eq. (20). Let {D,D′} denote the randomly split halves of the samples, where |D| = |D′| = n. Let η̂ = {π̂, ξ̂, θ̂} denote the estimators for the nuisance η using D′. Then, the MLTE estimator for β, denoted β̂, is given as a solution satisfying ED [ φ(β̂; η̂, ψ̂) ] = oP (n −1/2).
To illustrate, we exemplify Eq. (18) and Lemma (7, 8) for the case where Df is a KL-divergence and g(y;β = {µ, σ2}) is a normal distribution. First, m(β;ψ) = {mµ(µ;ψ),mσ(σ2;ψ, µ)}, where mµ(µ;ψ, σ) = (1/σ 2) (ψ[Y ]− µ) and mσ(σ2;ψ, µ) = (0.5/σ4) ( σ2 − ψ[(Y − µ)2] ) . We note
that µ̂m ≡ ψ̂[Y ] and σ̂2m ≡ ψ̂[(Y − µ̂)2] are estimators for β0 = {µ0, σ20} for the score m(β;ψ). Also, Rf (Y ;β, ψ) ≡ −(∂/∂β) log (g(Y ;β)) = {Rf (Y ;µ, ψ), Rf (Y ;σ2, ψ)}, where Rf (Y ;µ, ψ) ≡ (µ− Y )/σ2 and Rf (Y ;σ2, ψ) ≡ 0.5{σ2 − (Y − µ)2}/σ4. Then, the Neyman orthogonal score is given as φ(µ;σ2, η, ψ) = (1/σ2) {µ− ψ[Y ]− ϕ(η, ψ)[Y ]} and φ(σ2;µ, η, ψ) = (0.5/σ4) { σ2 − ψ[(Y − µ)2]− ϕ(η, ψ)[(Y − µ)2] } . Finally, solutions for φ(µ;σ2, η, ψ) and
φ(σ2, µ; η, ψ) are given by (µ̂, σ̂2), where, for ϕ[·] in Eq. (9), µ̂ = ψ̂[Y ] + ED [ ϕ(η̂, ψ̂)[Y ] ] and
σ̂2 = ψ[(Y − µ̂)2] + ED [ ϕ(η̂, ψ̂)[(Y − µ̂)2] ] .
The MLTE estimator in Def. 3 is consistent provided that nuisances estimates η̂ are consistent [14, Thm.4]. Such β̂ is known to achieve debiasedness [13], since β̂ is a DML estimator. Specifically,
Theorem 3 (Convergence rate of β̂). Let φ(β′; (η = {π, ξ, θ}, ψ) be given in Eq. (20). Let ϕm(β, η, ψ) be given in Eq. (19). Let β0, η0, ψ0 denote the true parameters. Let β̂ be the MLTE estimator for β defined in Def. 3. Suppose (1) Rf (y;β, ψ) is bounded and R′f (y;β, ψ) ≡ (∂/∂ψ)Rf (y;β, ψ) < ∞; (2) There exists a function H(y) < ∞ s.t. supβ,ψmax{Rf (y;β, ψ), R′f (y;β, ψ)} = O (H(y)); (3) {φ(β; (η, ψ))} is Donsker4 w.r.t. β for
4A function class where complexities are restricted. See Def. S.1 in the Appendix for the definition. Donsker class include Sobolev, Bounded monotone, Lipschitz class, etc.
the fixed η; (3) The estimators are consistent: β̂ − β0 = oP (1) and ∥ν − ν̂∥ = oP (1) for ν ∈ {πz(w), ξx(z, w), θ(x, z, w)[H(Y )]} for all (w, z, x, y); and (4) EP [φ(β; (η, ψ))] is differentiable w.r.t. β at β = β0 with non-singular matrix M(β0, (η, ψ)) ≡ (∂/∂β)|β=β0EP [φ(β; (η, ψ))] for all (η, ψ), where M(β0, (η̂, ψ̂)) P→M ≡M(β0, (η0, ψ0)). Then,
β̂ − β0 = −M−1ED [ϕm(β0; (ψ0, η0))] + oP (n−1/2) +OP (Rm2 ),
where Rm2 = ∑ z ( ∥π̂z − πz∥ {∥∥∥θ̂z − θz∥∥∥+ ∥∥∥ξ̂z − ξz∥∥∥}+ ∥∥∥ξ̂z − ξz∥∥∥2 + ∥∥∥θz − θ̂z∥∥∥2 + ∥∥∥ξ̂z − ξz∥∥∥∥∥∥θz − θ̂z∥∥∥) , where πz ≡ πz(W ), ξz ≡ ξx(z,W ), and θz ≡ θ(x, z,W )[H(Y )]. Corollary 4 (Debiasedness property for β̂). If nuisances {π̂, ξ̂, θ̂} converges at n−1/4 rate, then the target estimator β̂ converges to β0 at √ n-rate.
For the above example where Df is the KL divergence and g(y;β) is a normal distribution, H(Y ) = Y for Rf (y;µ, ψ), and H(Y ) = Y 2 for Rf (y;σ2, ψ).

heading: 5 Empirical applications
text: In this section, we apply the proposed methods to synthetic and real datasets. For the kernel-smoothing based approach, we compare KLTE with a baseline plug-in estimator (‘kernel-smoothing’), where estimates of nuisances η̂ = {π̂, ξ̂, θ̂} are plugged in the estimand Eq. (2). We use the Gaussian kernel. The bandwidth is set to h = 0.5n−1/5. In estimating the density, we choose 200 equi-spaced points {y(i)}200i=1 in Y and evaluate both estimators at Kh,y(i) for i = 1, · · · , 200. For the model-based approach, we compare MLTE (e.g., µ̂, σ̂2) with a moment-score-based estimator (called ‘moment’), defined as β̂m satisfying m(β̂m; ψ̂) = oP (n−1/2) (e.g., {µ̂m, σ̂2m}). We use KL divergence for Df and the normal distribution for g(y;β). For both approaches, nuisances are estimated through a gradient boosting model XGBoost [11], which is known to be flexible.

heading: 5.1 Synthetic dataset
text: We applied the proposed estimators to estimate the LTE p(yx|C) where the true densities are given as in the 4th plot in Fig. 2. As shown in the ground-truth in Fig. 3a, true densities p(yx0 |C), p(yx1 |C) are given as a mixture of four Gaussians. Estimated densities for Moment and MLTE are given in Fig. (3b, 3c). We note that model-based approaches fail to capture important characteristics (such as the number of modes) of the true density (‘ground-truth’ in Fig. 3a) because the assumed density class is misspecified. The ‘kernel-smoothing’ (Fig. 3d) captures only one of the modes having the highest densities, and this leads to misinterpretation of the true densities. KLTE (Fig. 3e) is able to capture the number, location, and scales of modes correctly.

heading: 5.2 Application to 401(k) data
text: We applied the proposed estimators (KLTE and MLTE) on 401(k) data, where the data generating processes corroborate with Fig. 1. Monotonicity assumption holds naturally, since ineligible units (Z = 0) cannot participate (X = 1) in 401(k). In our analysis, we used the dataset introduced
by [2] containing 9275 individuals, which has been studied in [2, 17, 5, 47, 58, 64], to cite a few. Model-based approaches (Moment in Fig. 4a and MLTE in Fig. 4b) and kernel-smoothing based approaches (kernel-smoothing in Fig. 4c and KLTE in Fig. 4d) are implemented to analyze the data.
The model-based (Fig. (4a,4b)) and kernel-smoothing based (Fig. (4c,4d)) estimates both capture important characteristics of the distribution, such as mode, location, and scale parameters. The results of proposed estimators (MLTE and KLTE in Fig. (4b,4d)) are consistent with findings from previous analyses [2, 17, 5, 58]: The effects of the 401(k) participation (i.e., X = 1) on net financial assets are positive over the whole range of asset distributions. To connect to CDF method, we provide in Fig. 4e the CDF estimate induced by KLTE density estimation (Fig. 4a). We note that the CDF in Fig. 4e captures the nonconstant impact trend of the 401(k) participation on the net financial assets, which has been also described in the previous analyses [2, 17, 5, 58].

heading: 6 Conclusion
text: In this paper, we develop kernel-smoothing-based and model-based approaches for estimating the LTE density in the presence of instruments. For each approach, we give Neyman orthogonal scores (Lemma (2,8)) and constructed corresponding DML estimators (KLTE in Def. 1 and MLTE in Def. 3), that exhibit debiasedness property (Corol. (3, 4)). We demonstrated our work through synthetic and real datasets. The performance of model-based estimators depends critically on the choice of the density class. Kernel-based estimators do not have to make assumptions about the true density class but will suffer from the curse of dimensionality. This work is limited to settings where the monotonicity assumption holds, i.e., there are no defiers. One could perform sensitivity analyses on the impact of potential defiers to the estimates as conducted in [65, 36].

heading: Acknowledgements
text: We thank the reviewers for their feedback helping to improve this manuscript. Elias Bareinboim and Yonghan Jung were partially supported by grants from NSF IIS-1750807 (CAREER). Jin Tian was partially supported by ONR grant N000141712140.

Review summaries:
This paper proposes two methods for estimating the density of the outcome variable for the compliers in the potential outcomes framework. The first is a kernel-smoothing based approach, while the second is a model-based approach, and in both cases, the double machine learning approach based on data splitting is employed. The nature of the contribution seems mostly theoretical, with a lot of effort into proving the convergence rates. 

This paper deals with the problem of the estimation of the local treatment effect (LTE) which is measure of the causal impact of a treatment. To this extent the authors develop kernel-smoothing and model-based methodologies to recover this quantity from the data.  The theoretical properties of the machinery are developed and the authors show that the estimator of the LTE is consistent. Two applications of the proposed methodology showcase the technique.

The authors propose kernel-smoothing-based and model-based approaches for estimating the LTE density in the presence of instruments. They derive Neyman orthogonal scores and construct the corresponding DML estimators (KLTE and MLTE), that exhibit debiasedness property. They experiment on synthetic and real datasets.

The most common problem in causal effect estimation is that of estimating average effects over populations. However, to get to individual level causal effects, one has to estimate the full *distribution* or densities in the continuous outcome case. The authors are estimating such a density for compliers in the IV setting. Two density estimation strategies are provided i) kernel-based, ii) model-based. Essentially the core contribution is that in the identifiable setting (under monotonicity), authors have extended density estimation methods to the IV setting for estimating density of causal effects for compliers. The key methodology used is double machine learning which has extensive precedent for designing robust estimators for causal inference problems. I think the contribution is valuable to the community and toward making causality realistic and using it for real-world applications.