Introduction

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 treatments X = 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. 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., p n-rate) even when nuisance estimates converge slowly (e.g., nfoot_0/4 rate) (a property called 'debiasedness' 1 ). We illustrate the proposed methods on synthetic and real data.

Related work

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].

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 Y y (Y ), an indicator that outcome Y falls short of threshold y, reduces to the CDF (i.e., replacing Y in LATE with Y 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 p 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 p 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.

LTE Estimation -Problem setup

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 x (X) to represent the indicator function such that x (X) = 1 if X = x; x (X) = 0 otherwise. For a continuous variable X with a probability density p(x) of a distribution P and a function f

where X is the domain for X, and kf (X)k ⌘ p

Structural Causal Models (SCMs). We use the language of SCMs as our basic semantic and inferential framework [50,4]. An SCM M is a quadruple M = hU, V, P (U ), F i 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 = {f Vi } Vi2V such that V i f Vi (pa i , u i ) where P A i ✓ V and U i ✓ U . Each SCM M 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 P A i to V i 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., f X (pa x , u x )) by the constant x and induces a submodel M x and an interventional distribution P (v|do(x)). For any variable Y 2 V , the potential response Y x (u) is defined as the solution of Y in the submodel M x given U = u, which induces a counterfactual variable Y x . Local Treatment Effect (LTE) with IV. We consider the IV setting represented by the causal graph G in Fig. 1 2 , 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 ⇢ R d , 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 .

The causal density p(y x ) 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 [3,2]. We will make the following assumptions based on literature. Assumption 1 (Monotonicity). There are no defiers: X Z=1 X Z=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 X Z=0 = 0 and X Z=1 = 1). For a given constant a and a variable X, let x a denote the event X = a. The LTE p(y x |C) is identifiable under monotonicity and is given by [31,2]:

where the expectation is over W . In this paper, we aim to estimate the LTE density p(y x |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(y x |C) are bounded, and p(y x |C) is twice differentiable.

DML method. Let ⌘ P 0 denote a functional of an arbitrary distribution P 0 . 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 0 is a set of nuisance parameters and ⌘ 0 ⌘ ⌘ P denotes the true nuisances), a function such that E P ['(V ; 0 , ⌘ 0 )] = 0, (@/@⌘)| ⌘=⌘0 E P ['(V ; 0 , ⌘)] = 0. Given ', an DML estimator is constructed using the cross-fitting technique as follows: For randomly split halves of

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 b ⌘ 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 E P [ ] = 0, E P ⇥ 2 ⇤ < 1, and (@/@t) Pt | t=0 = E P [ (V ; , ⌘)S t (V ; t = 0)] where P t ⌘ P (v)(1 + tg(v)) for t 2 R and any bounded mean-zero functions g(•) over V , and S t (v; t = 0) ⌘ (@/@t) log P t (v)| t=0 [63,Chap. 25].

Due to space constraints, all the proofs are provided in Appendix B in suppl. material.

Kernel-smoothing-based approach

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). For convenience, we denote the target estimand by (y) ⌘ p(y x |C). In the kernel-smoothing-based approach, we will aim to estimate a kernel-smoothed approximation for (y) defined as follows:

where [f (Y )] is an expectation of a function f (Y ) w.r.t. (y), which is specified as

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

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

Then, the following is a moment score function for h :

where h is given in Eq. ( 2) and 0 is an estimate of h .

Next, we derive an influence function for the moment score function m( 0 ; h ). We first define the following function: for a bounded function

and

where V X is defined in Eq. ( 5). Then, the influence function for the expectation of the moment score function m( 0 ; h ) in Eq. ( 6) is given as follows:

Lemma 1 (Influence function for m( 0 ; h )). Let m( 0 ; h ) be the score defined in Eq. ( 6). Then, the influence function for E P [m( 0 ; h )], denoted m , is given by

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 scorefoot_2 ([14, Thm.1], [13, Sec. 2.2.5]). Specifically, Lemma 2 (Neyman orthogonal score for h ). Let m( 0 ; h ) be the score function in Eq. ( 6), and m (⌘ = {⇡, ⇠, ✓}, h ) be the influence function for E P [m( 0 ; h )] given in Eq. (10). Then, a Neyman orthogonal score for h is given as '( 0 ; ⌘ = {⇡, ⇠, ✓}) ⌘ m( 0 ; h ) + m (⌘, ); Specifically,

Given the Neyman orthogonal score '( 0 ; ⌘), an estimate ˆ h satisfying

) gives a DML estimator. Specifically, we propose the following kernel-smoothing based estimator for the LTE density, named 'KLTE' (kernel-based estimator for LTE): Definition 1 (KLTE estimator for h ). Let '( 0 ; ⌘ = {⇡, ⇠, ✓}) be the Neyman orthogonal score for h given in Eq. (11). Let {D, D 0 } denote the randomly split halves of the samples, where |D| = |D 0 | = n. Let ⌘ = {⇡, ⇠, ✓} denote the estimates for the nuisance ⌘ using D 0 . Then, the KLTE estimator for h (y) for all y 2 Y, denoted ˆ h (y), is given by

where V X and V Y 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.

Asymptotic convergence

Now, we study the convergence rate of the estimator ˆ h (y). For any fixed y 2 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 2 Y, suppose the estimators for nuisances are consistent; i.e., k⌫ ⌫k = o P (1) for ⌫ 2 ⌘ = {⇡, ⇠, ✓} for all (w, z, x). Suppose h < 1, and nh d ! 1 as n ! 1. Then,

where

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 {nh d } 1/4 , then the target estimator ˆ h (y) achieves p nh d -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:

Combining the results of Lemma. (3,4), we have the following result: Theorem 1 (Convergence rate of ˆ h to ). For any fixed y 2 Y, suppose the estimators for nuisances are consistent; i.e., k⌫ ⌫k = o P (1) for ⌫ 2 ⌘ = {⇡, ⇠, ✓} for all (w, z, x). Suppose h < 1, and

where B y is defined in Eq. ( 14), and R k 2 is defined in Eq. ( 13).

Thm. 1 implies that ˆ h (y) converges fast (see Corol. 1) to (y) + B y . 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., nh d ! 1). So far, we have analyzed the error ˆ h (y) (y) pointwise for the fixed y 2 Y. To analyze the difference between the two densities ˆ h (y) and (y) for all y 2 Y, we consider the following divergence function of two densities: Definition 2 (f -Divergence D f [20]). Let f denote a convex function with f (1) = 0. D f (p, q) ⌘ R 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, D f reduces to KL divergence with f (p, q) = (p/q) log(p/q). We will assume that the function f (p, q) in D f is differentiable w.r.t. p and q.

Recall that

We now analyze the distance between ˆ h and w.r.t. D f . The following result provides an upper bound for D f . Lemma 6 (Upper bound of the divergence

where w(y) ⌘ f 0 2 ( (y), ˜ (y)) ˆ h (y), f 0 2 (p, q) ⌘ (@/@q)f (p, q), and ˜ h (y) ⌘ t ˆ h (y) + (1 t) (y) for some fixed t 2 [0, 1].

By invoking Thm. 1, we derive an upper bound for D f ( , b h ) as follows:

Theorem 2 (Convergence rate of ˆ h ). Suppose the estimators for nuisances are consistent; i.e., k⌫ ⌫k = o P (1) for ⌫ 2 ⌘ = {⇡, ⇠, ✓} for all (w, z, x, y). Suppose D f is a f -divergence such that f (p, q) = 0 if p = q. Suppose w(y) in Lemma 6 is finite. Then,

where R k 2 is defined in Eq. ( 13) and B y is defined in Eq. ( 14).

The following result asserts that the debiasedness property is exhibited w.r.t. D f :

Suppose w(y) in Lemma 6 is finite. If nuisances {⇡, ⇠, ✓} converges at {nh d } 1/4 rate for any (w, z, x, y), then D f ( , b h ) converges to 0 at p nh d -rate.

Model-based approach

In this section, we develop a model-based approach for estimating the LTE density (y) = p(y x |C).

We will approximate with a class of distributions or a density model G = {g(y; ) : 2 R b } where g(y; ) 2 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 ) 2 G that achieves the minimum f -divergence to :

where D f is the f -divergence defined in Def. 2. Our goal is estimating 0 .

Consider m( ; ) ⌘ (@/@ )D f ( (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]:

where g 0 (y; ) = (@/@ )g(y; ) and f 0 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

where (⌘, )[•] is defined in Eq. ( 9), and

))} , where g 0 (y; ) ⌘ (@/@ )g(y; ), f 0 1 (p, q) ⌘ (@/@p)f (p, q) and f 00 21 (p, q) ⌘ (@/@p)f 0 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 '( 0 ; (⌘ = {⇡, ⇠, ✓}, )), is given by

where m ( , ⌘, ) is defined in Eq. (19).

Given the orthogonal score '( 0 ; (⌘, )) in Eq. ( 20 To illustrate, we exemplify Eq. ( 18) and Lemma (7,8) for the case where D f is a KL-divergence and g(y; = {µ, 2 }) is a normal distribution. First, m( ; ) = {m µ (µ; ), m ( 2 ; , µ)}, where

] are estimators for 0 = {µ 0 , 2 0 } for the score m( ; ). Also, R f (Y ; , ) ⌘ (@/@ ) log (g(Y ;

) 4 . Then, the Neyman orthogonal score is given as

. Finally, solutions for '(µ; 2 , ⌘, ) and '( 2 , µ; ⌘, ) are given by (μ, ˆ 2 ), where, for [•] in Eq. ( 9),

i and

i .

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 '( 0 ; (⌘ = {⇡, ⇠, ✓}, ) 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) R f (y; , ) is bounded and R 0 f (y; , ) ⌘ (@/@ )R f (y; , ) < 1; for all (⌘, ), where M ( 0 , (⌘, ˆ ))

Empirical applications

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) } 200 i=1 in Y and evaluate both estimators at K h,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 ; ˆ ) = o P (n 1/2 ) (e.g., {μ m , ˆ 2 m }). We use KL divergence for D f 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.

Synthetic dataset

We applied the proposed estimators to estimate the LTE p(y x |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(y x 0 |C), p(y x 1 |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. (e) KLTE Figure 3: LTE estimation with a synthetic dataset. The ground-truth density is in (a). Red and Green for x 0 and x 1 , respectively.

Application to 401(k) data

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].

Conclusion

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].