Normalizing Flows for Interventional Density Estimation
Authors: Valentyn Melnychuk
Dennis Frauen
Stefan Feuerriegel
Venue: Submitted to ICLR 2023 (2023)
Keywords: causal inference
normalizing flows
treatment effect estimation
causal machine learning

Abstract: Existing machine learning methods for causal inference usually estimate quantities expressed via the mean of potential outcomes (e.g., average treatment effect). However, such quantities do not capture the full information about the distribution of potential outcomes. In this work, we estimate the density of potential outcomes after interventions from observational data. For this, we propose a novel, fully-parametric deep learning method called Interventional Normalizing Flows. Specifically, we combine two normalizing flows, namely (i) a teacher flow for estimating nuisance parameters and (ii) a student flow for a parametric estimation of the density of potential outcomes. We further develop a tractable optimization objective via a one-step bias correction for an efficient and doubly robust estimation of the student flow parameters. As a result our Interventional Normalizing Flows offer a properly normalized density estimator.  Across various experiments, we demonstrate that our Interventional Normalizing Flows are expressive and highly effective, and scale well with both sample size and high-dimensional confounding. To the best of our knowledge, our Interventional Normalizing Flows are the first fully-parametric, deep learning method for density estimation of potential outcomes.

Under review as a conference paper at ICLR 2023
NORMALIZING FLOWS FOR INTERVENTIONAL DENSITY ESTIMATION
Anonymous authors Paper under double-blind review
ABSTRACT
Existing machine learning methods for causal inference usually estimate quanti- ties expressed via the mean of potential outcomes (e.g., average treatment effect). However, such quantities do not capture the full information about the distribu- tion of potential outcomes. In this work, we estimate the density of potential outcomes after interventions from observational data. For this, we propose a novel, uully-parametric deep learning method called Interventional Normalizing Flows. Specifically, we combine two normalizing flows, namely (i) a teacher flow for estimating nuisance parameters and (ii) a student flow for parametric estimation of the density of potential outcomes. We further develop a tractable optimization objective based one-step bias correction for an efficient and doubly robust esti- mation of the student flow parameters. As a result our Interventional Normalizing Flows offer properly normalized density estimator. Across various experiments, we demonstrate that our Interventional Normalizing Flows are expressive and highly effective, and scale well with both sample size and high-dimensional con- founding. To the best of our knowledge, our Interventional Normalizing Flows are the first parametric, deep learning method for density estimation of potential outcomes.
INTRODUCTION
Causal inference increasingly makes use of machine learning methods estimate treatment effects from observational data (e.g.. van der Laan et al. 2011 Künzel et al. 2019 Curth & van der Schaar 2021: Kennedy| 2022). This relevant for various fields including medicine (e.g., |Bca et al. [002]] marketing (c.g., Yang al.) 2020), and policy-making (e.g.. Hünermund et al. (021). l). Here, causal
van
inference from observational data promises great value, especially when experiments for determining treatment effects costly even unethical.
The vast majority of the machine learning methods for causal inference estimate averaged quantities expressed by the (conditional) mean of potential outcomes. Examples of such quantities the average treatment effect (ATE) (e.g.. Shi et al. 2019 Hatt & Feuerriegel 2021), the individual treatment effect (ITE) (e.g., (halli al. 2017 Hassanpour & Greiner [019 Zhang al. 2020). and treatment-responee curves (e.g., Bica et al. 2220 Nie et al. (021). Importantly, these estimates only
describe averages without distributional properties.
However, making decisions based on averaged causal quantities can be misleading and, in some applications, even dangerous (Spiegelhalter [0177] van der Bles al. 2019). On the one hand, if potential outcomes have different variances or number of modes, relying on the average quantities provides incomplete information about potential outcomes, and may inadvertently lead local and not global - optima during decision-make making.. On the other hand, distributional knowledge is needed to account for uncertainty in potential outcomes, and thus informs how likely a certain outcome is. For example, in medicine, knowing the distribution of potential outcomes is highly important (Gische) Voelkle, 2021): it gives the probability that the potential outcome lies a desired range, and thus defines the probability of treatment success or failure. Motivated by this, aim to estimate the density of potential outcomes.
An example highlighting the need for estimating the density of potential outcomes is shown in Fig, Here, we simulated outcomes according to given structural causal model (SCM). The potential outcomes Y [[a] can be sampled by setting the treatment to specific value in the equation for (cf.

---

Under review as a conference paper at ICLR 2023
(4)
(1)
[d]
M1 M
F12 a
BYY vid
Minare (0.5 (3)
02-
711 -1
POYT -
POTT m1
A Bern
N(X:0.1)
0-3
((A +
64
(1 () 1.5):
Figure 1: Motivating example showing the densities of observational, interventional, and counter- factual distributions of outcome Y. These are simulated via the structural causal model on the right (here: N (x;; μ4, )²) are densities of the normal distribution; and 6 = 3 is a covariates shift, which regulates the probability of treatment assignment). Potential outcomes have different distributions but the same mean (0)) [{{[1]) = 4.77 and the same variance var(Y(o),
var x 4.06. Here, Y(a [[a is the potential outcome given treatment a. (a) Interventional distributions. (b) and (c) Ob- servational and counterfactual distributions for the same outcomes. As shown here, the observational, interventional, and counterfactual distributions can be vastly different.
Appendix B) At the same time, by flipping the treatment assignment in this equation, we obtain counterfactual outcomes Y [a] A a'. We observe that the potential outcomes have the same mean (i.e., (0) = [1(( and the same variance (i.e., var(Y[[)) = var(Y(1))) Hence, the ground-truth ATE equals zero. Nevertheless, the distributions of potential outcomes (i. e., P(Y(a)) are clearly different. Hence, in medical practice, acting upon the ATE without knowledge of the distributions of potential outcomes could have severe, negative effects. To show this, let us consider a "do nothing" treatment (a = 0) and some medical treatment (a 1). Further, let us consider an outcome be successful if some risk score Y is below the threshold of five. Then, the probability of treatment success (i. e. [1] < 5.0) R 0.63) is much larger than the probability of success after the "do nothing" treatment (i. PPY (0) < 5.0} x 0.51), highlighting the importance of treatment.
In this paper, we aim to estimate the density of potential outcomes after intervention a, i.e P(Y y). From this point on, we refer to this task as interventional density estimation (IDE). Estimating the density of interventions has several crucial advantages: it allows to identify multi-modalities in the distribution of potential outcomes; it allows to estimate quantiles of the distribution; and allows to compute the probability with which potential outcome lies in certain range. Importantly, traditional density estimation methods are not applicable for IDE due to the fundamental problem of causal inference: that is, the counterfactual outcomes are typically never observed, and, hence, the sample from ground-truth interventional distribution is also inaccessible.
In prior literature, Kennedy et al. (2021) introduced theory for efficient semi- cani-parametric IDE estimation, but without a flexible algorithmic instantiation in form of a method. Existing literature also offers some specific methods for IDE, which are either semi- or non-parametric. Examples are kernel density estimation (Kim et al., 2018) and kernel mean embeddings of distributions (Muandet et al. 2021). However, both methods neither scale well with the sample size nor with the dimensionality of covariates. Furthermore, both methods have an additional, crucial limitation: estimated densities could be unnormalized or even return negative values (which, by definition, is not possible). mully-parametric methods, on the other hand, have several practical advantages: they automatically provide properly normalized density estimators, they allow one to sample from the estimated density and typically scale well with large and high-dimensional datasets. However, to the best of our knowledge, there is no aally-parametrcc, deep learning method for IDE.
In this paper, we develop a novel, tilly--aramttric deep learning method: Interventional Normalizing Flows (INFs). Our INFs build upon normalizing flows (NFs) (Tabak & Vanden- Vanden-Eijnden 2010 Rezende & Mohamed 2013). but which carefully adapt for causal inference. This requires several non-trivial adaptations. Specifically, we combine two NFs: a (i) teacher flow for estimating nuisance parameters, and (ii) student flow for a parametric estimation of the density of potential outcomes. Here, we construct a novel, tractable optimization objective based on a one-step bias correction to allow for an efficient and doubly robust estimation. At the end, we develop two-step training procedure to train both the teacher and the student flows
Overall, our main contributions are following
'e distinguish the interventional distribution (i.c., P(Y[a)) and the counterfactual distribution (i.c., P(Y[a) A = a') d')), which are different in general. This can be seen by comparing plats (a) vs. (b) and (c) in Fig. M For further information, we refer to Appendix [
Cooee available at/httss: //tp////nnnnuusu/////// 4open. wcciencerr//e /tccieen//////nmooous Interrl ow R2F3
2

---

Under review as a conference paper at ICLR 2023
1. We introduce the first palameeraeetrcc, deep learning method for interventional density estimation, called Interventional Normalizing Flows (INFs). Our INFs provide a properly normalized density estimator.
2. We derive a tractable optimization problem with a one-step bias correction for efficient and doubly robust estimation. To solve, we propose a two-step training procedure with our INFs.
3. We demonstrate in various experiments that our INFs are highly expressive and effective. A major advantage owed to the parametric form of the student flow is that our INFs scale well to both large and high-dimensional datasets in comparison other non- and semi-parametric methods.
2 RELATED WORK
Recently, there has been great interest in using machine learning and, specifically, deep learning for estimating causal quantities. Examples are machine learning for estimating ATES (e.g.. Shi et al.
2019; Hatt& & Feuerriegel] 2 ITBB (e.g., Johansson al., 2016 Alaa & van der Schaar 2018
Wager& & | 2018)] urth & van der Schaar (021), and treatment teveent-reppee response curves (c.g.. BBica et al. 2020 Schwab al., 2020 Nie et al. 2021). In this regard, some papers proposed uncertainty-awe
methods, c. g., by using the variance of potential outcomes (Alaa & van der Schaar 2017 Jesson
et al. 2020). or the conditional outcome distribution (Jesson et al.] [2021] (022). However, the aforementioned works are all concerned with estimating averaged causal quantities expressed via the mean of potential outcomes. In contrast, there are only a few papers that estimate the density of outcomes after intervention.
2.1 INTERVENTIONAL DENSITY ESTIMATION
Kennedy et al. 2021 introduced theory for efficient semi semi-parametric estimation. The theory also lends to a hypothetical estimator as a solution to an integral equation, namely a bias-corrected moment condition. However, the theory comes without an algorithmic instantiation in form of method. We later adopt the theoretical framework and convert the bias-corrected moment condition into a tractable optimization objective, which we can then solve very effectively with deep learning.
T lists existing methods for IDE. Importantly, these are either conmmaaamttric or semi- parametric. Kim et al. 2018) developed a doubly robust kernel density estimation (KDE) via an efficient estimation of density functionals. MMandett et al.] (2021) proposed kernel mean embeddings of distributions (DKME), which provides a non-parametric plug-in estimator. However, both methods (Kim et al. 2018 Muandet et al., 202T) have limitations. First, they do not provide properly normalized density estimator. Hence, the estimated densities can be unnormalized or even negative, yet which, by definition, is not possible. Second, they do not offer direct sampling, which would allow one to sample from the estimated density without an additional algorithm. This may complicate computations of the test propobbability empirical Wasserstein distance during evaluation. Third, another limitation of both non-parametric and semi-paramen methods is that they typically scale not well. This is unlike parameera methods, which scale well both large and gigh-iimensional datasets. However, so far, there is no full-parametric, deep learning method for IDE.
The methods for IDE above (Kim et al. 2018 Muandet al. 2021 Kennedy et al. 2021 build upon general assumptions for causal identifiability. We later adopt the same assumptions for IDE (see deep learning method called INFs. Our method
Section 3) and we then develop a has three favorable properties: it yields a proper density estimator, it allows for direct sampling, and scales well.
Table 1: Overview of methods for interventional density estimation from observational data.
Method
Density model
Parametric
Estimator type
Proper dessity
Diseet sampling
Kim et al 2018
kernel density estimation (KDE)
semi parametric
A-IPTW
X
x
Mummar 3221
dieributioral kond mean embeddinnss (DKMF)
woe-paraeeric
plug-in
INFs (this poper)
noemalizing nows (NP4)
fully-parametic
A-IPTW
A-IPTW; nagieemed invense propersity treameent weighted
2.2 EFFICIENT ESTIMATION
In the context of treatment effect estimation, so-called augmented inverse propensity of treatment weighted (A-IPTW) estimators were developed for efficient, parameereet estimation of target estimands (parameters) (Robins 2000). Formally, A-IPTW estimation performs a first-order bias
correction of plug-in models (Bickel et al. 1993 Chernozhukov et al. 2018). A-IPTW estimation
3

---

Under review as a conference paper at ICLR 2023
also offers the property of being double robust, i. e., fast convergence rates even if one of the nuisance parameter estimators converges slowly Kennedy (2020). Although Kennedy et al. (2021) formulated an integral equation for semi- and fully- parametric efficient IDE estimation (see Eq. 19 therein), no
flexible algorithmic instantiations in form of a method have been implemented so far. Later, we reformulate this equation tractable optimization problem, and, thereby, turn our INFs into efficient and doubly robust IDE estimator.
2.3 NORMALIZING FLOWS
Normalizing flows were introduced for expressive variational approximations in variational autoen- coders (Tabak & Ganden-Eijnden 2010; Rezende & Mohamed| 2015). One practical benefit of NEs is that they yield universal density approximators (Dinh et al. 2014; 2017 Huang et al., 201 Durkan et al. 2019). Furthermore, NE's can be leveraged for conditional density estimation (e. g.. via so-called hypernetworks (Trippe & Turner 2018). Normalizing flows were used for causal inference, but in different setting from ours (see Appendix A We provide background on normalizing flows in Appendix B
Research gap:
Existing methods for IDE are either knowledge, our work is the first to propose a fully-parametricc deep learning method for IDE.
semi-parametro To the best of
SETUP: INTERVENTIONAL DENSITY ESTIMATION
Notation. Let P(Z) be a distribution of a random variable Z, and let
= P be its density or probability mass function. Let T. = P(A a X = x) denote the propensity score. Further, 1() the indicator function; P,, (x = n Ei= f(X;) is the sample average of random f(X): and I [[(X is the average evaluated on a minibatch of size b. For readability, we sometimes highlight random variables and the corresponding averaging operator in green color. Furthermore, P(Y X, A) is the conditional distribution the outcome Y.
Problem statement. In this work, aim at estimating the interventional density from observational data, namely P(Y y). To compare the goodness-of-fit of different estimators, we evaluate the distributional distance between the ground-truth interventional density and the estimated density. Such distributional distances include, e.g., the average log- probabbilityyy and the empirical Wasserstein distance.
We build upon the standard setting of potential outcomes framework (Rubin) 1974), where [[]] stands for the potential outcome after intervening on treatment by setting it to a. That is, we consider an observational sample with d -dimensional covariates X E C Rdx treatment A E (0, and a dy -immensionall continuous outcome Y E y C Rd drawn i.i.d. We consider dy = if not stated explicitly. We assume the treatment to be binary, but note that our INFs also work with categorical treatments. We denote D = (Xi, A;, Y ~ IP(X, A, Y), where is the sample size, and i the index of an observation. For example, in critical care, the patient covariates X are different risk factors (c.g., age, gender, weight, prior diseases), the treatment is whether ventilator is applied, and the outcome is the probability of patient survival. The covariates X are also called confounders if (al * TP(Y A= a).
Identifiability. To identify the interventional density, we make the following identifiability as- sumptions with respect to the data-generating mechanism of D: (1) Positivity: For some 6 >> 0, >
1. (2) Consistency: If A = a for some patient, then Y = [a). (3) Ex- [ X for all a. Note that these assumptions are standard in the literature (Kim et al. 2018 Kennedy al. 2021 Muandet et al. (021). Under assumptions (1)-(), the density [a] can be expressed terms of observational distribution with
P(1- > (X)
changeability: A
Y
of interventional distribution
back-door adjustment, i.c..
P(Y v) JEE
/ P(X =2
B X~P(X)
(P(Y
xx,A=A=a)) (I)
where PYY = y A) is the conditional density of the outcome. For more details on the potential outcomes framework and identifiability, we refer to Appendix B
Plug-in estimator. A straightforward approach for IDE (Robins & Rotnitzky 2001 is the follow-
ing: first, one estimates the conditional outcome distribution, P(Y X, A) (here, any method for
4

---

Under review as a conference paper at ICLR 2023
Studeing
CNF (CCCA
Pry s(x,A)
NF 2
Day
pestesalily
Yor
vl
Toster
bockin
Noli rgga attatin
Teacher flow
Student flow
Figure 2: Overview of Interventional Normalizing Flows. Our INFs combine two normalizing flows, which we call "teacher flow" and "student flow". 'The teacher flow estimates the nuisance parameters, Le., the propensity score it (X) and the conditional outcome distribution P(Y X, A). The student flow utilizes them estimate the projection parameters BA IPTW Both teacher and student flows are fitted via a two-step training procedure.
conditional density estimation could be used). Then, one takes sample average over covariates X:
1n [Y[aa y) = P. X
(2)
This estimator is unbiased but inefficient estimator of interventional density, which is known as semi-parametric plug-in estimator. Semi-parametric IDE, unlike, e. e.g., semi-parametric ATE estimation, is highly problematic. For large sample sizes, the semi-parametrie estimator requires averaging over the full sample for each evaluation point. Motivated by this, we aim to develop fully-parametric estimator.
THEORETICAL BACKGROUND FOR FULLY-PARAMETRIC IDE
In this section, we introduce theory for mully-parametric estimation of interventional density. First, we describe parametric plug-in estimator as a solution to the moment condition (Kennedy et al. 202T We call this estimator covariate-adjuster estimator. Second, develop one-step bias correction for efficient estimation.
We start by defining parametric model, {glv; 8a) Ba E Rd), where Sa IRdd are parameters of estimator, and g( Ba) is density, i. e., Syep I(1, Pa) dy = 1. For IDE, we approximate the interventional distribution [a]) with distribution from our parametric model. That is, we aim minimizing the distributional distance (specifically KL-divergence) between IP [al) and g(: 8a) via
3. arg min KL (((Y1))(((:..))) = arg min 9.
E (*)F(Y(a)) (- - gY ( (Ba))
(3)
Aa
where 8n are called projection parameters they project the true interventional density onto class (o(:B.)) 8a E Rd).
4.1 COVARIATE- ADJUSTED ESTIMATOR
Let the d-dimensional random variable Following Kennedy et al. (2021). the projection parameters can equivalently expressed as a solution
=
V (Y:; ;B,) denote the score function.
the moment condition m(33.) O, where
m(aa) = ("_)((()))
(a) x-F(X)
E (B(T(Y:B.) 11,11 a)). xx)
(4)
Here, the moment condition is the expected score function of the potential outcome. Throughout the paper, we assume that the moment condition has a unique solution, and, therefore, the minimization task in Eq. (3) and the root-finding task in Eq. (4) are equivalent.
In practice, have neither observations from the interventional distribution counterfactual outcomes. Therefore, we cannot use the ground-truth [a[(aa) but, instead, must use the plug-in estimator distribution from Eq. (2). Specifically, we can obtain a plug-in estimator of projection
parameters, i. e., EPI either by minimizing cross-entropy loss or by solving the moment condition, both of which are equivalent:
3P1 arg min 8..
E ya-P, (P(Y x.A=a))
(Y" vi" (Ba) yeaR, (P(1)X
E T(Y":3.) (3,) 1 0. (5)
5

---

Under review as a conference paper at ICLR 2023
Then, define parametric covariate-adjusten (CA) estimator PCA (Y[a] =) 9 ) By choosing a sufficiently expressive class of densities for both and the conditional density estimator P(Y X, A) (e. g., normalizing flows), CA can be shown to consistently estimate the interventional density (see Appendix B.5 in Kennedy al. (2021).
4.2 EFFICIENT ESTIMATION VIA ONE-STEP BIAS CORRECTION
In the following, aim develop efficient estimator of the projection parameter 3. from Eq. (3) or, equivalently, the moment condition rm(B))) ) at fixed B, from Eg. 44. For this, we make use of
semi-parametric efficiency theory (van der Laan & Robins 2003 Kennedy et 202T). We provide
background on efficiency theory in Appendix
Kennedy 2022) showed that the efficient influence function da (T, for the functional EE(((T
X,A a)) equals
6n (T; Ta (X) (T-
1(A=a)
E X, X-P(X)
(E(T) A,A = a)).
(6)
Here, red color to show the nuisance parameters of IP that are influencing the functional.
We emphasize that the nuisance parameters (i. c., the propensity score and conditional expecta- tions/probabilittes) can be either known estimated.
The efficient influence function in Eg. (6) allows to construct an efficient estimator of the moment l), we transform the plug-in estimator ( from Eq. (5) into efficient estimator with the help of one-step bias correction. In the bias-corrected
condition. Following (Kenneyy et al. 2021),
moment condition has the following form:
fit A-IPTW (3o) rn
where P(Y A)} are the estimated nuisance parameters of F. The estimated nuisance parameters are simultaneously used for plug-in estimation of the moment condition. We call the solution of the bias-corrected moment equation AA IPTW augmented inverse propensity of treafment weighted (A-IPTW) estimator the projection parameters. Unlike CA estimator, A-IPTW estimator achieves efficiency and possesses double robustness property.
We now transform the bias-corrected moment condition into the following tractable optimization task (see Appendix for details):
3 3APPTW
arg min Ra
E yamea (P (AAAa))
log
cross cross-entropy loss
+ L(A =()= (x(
(-logo(V;;
E Y~ x.A=
(logo(Y;; ;Ba))
(8)
ono-stop bins comection
Previously, Kennedy et al. 2021 proposed to directly solve bias clas-corrected moment condition, i. e., system of nonlinear equations, yet which is in general much harder to solve, even computationally. In contrast, we develop an optimization objective that can be directly incorporated into loss of deep learning density estimator.
INTERVENTIONAL NORMALIZING FLOWS
In the following, we describe our Interventional Normalizing Flows: a parly-parametric method for interventional density estimation via deep learning. First, we describe all the components of our architecture and, then, introduce an efficient estimation using one-step bias correction.
5.1 COMPONENTS
In our INFs, we combine two normalizing flows, which we refer to as (i) teacher flow and (ii) student flow (see Fig. 2). The rationale for this is based callmmearmmtticc IDE requires two models: (i) one for the estimation of nuisance parameters, and (ii) one for the subsequent optimization of the learning objective with respect to projection parameters.
derivations in Section 4 according which

---

Under review as a conference paper at ICLR 2023
Accordingly, both NFs in our INFs have thus different objectives: (i) the teacher flow estimates the nuisance parameters (i.e., the propensity score and the conditional outcome distribution); and (ii) the student flow uses the estimated nuisance parameters to estimate the projection parameters.
(i) Teacher flow. The teacher flow has three components: two fully-connected (FC) subnetworks and conditional normalizing flow parameterized by 0. The first FC subnetwork (FC1) takes the covariates X as input and, then, outputs a representation R E R" together with a propensity score (X). The second FC subnetwork (FC2) takes the representation It and the observed treatment A; as input and, then, outputs the parameters of flow, conditioned on X and A, i. i.. 0(X, A). Together, FC₁ and FC₂ form a so-called hypernetwork (Ha et al. 2017) for the conditional normalizing flow, which allows us to learn the conditional outcome distribution via back- pockkpropaggtion
Let L. be the loss of the teacher flow. Here, we combine a conditional negative log- likelkihioodd (ENLl) and binary cross-entropy loss for the propensity score (Cx), i.e., (IP, #a) Pr CNLL aCz} with LNLL = - log PY = Y A); C = ( A), where a > 0 is a hyper- parameter. In general, conditional normalizing flows are prone to overfitting when trained via conditional negative bog-likelihood. To address this, we later employ noise regularization (Rothfuss) al. 2019) in the conditional density estimation.
(ii) Student flow. The student flow uses the outputs of the teacher flow and then learns the inter- ventional distribution. We first describe the naive variant of the student flow without one-step bias correction (we introduce this later in Section 5.2p. Different from the conditional normalizing flow in the teacher flow, the student flow is a con-conditional normalizing flow, parameterized by Ba. Specifically, we consider two separate normalizing flows, that is, one for each potential outcome (i.e., a and a 1, respectively)
To fit the student flow, we must solve the moment condition from Eq. ( or, equivalently, minimize a cross-entropy loss (CCE). Here, we use a tractable approximation via numeric integration:
log "Y" = - [ | log (y: B.) (P(Y = y X, = a)} dy
Lcx(Bo) =
E yomp, ())))
h Eド= logg(y;:Bu) P. (É(Y Ui X, = a)}
dr =
(9)
Fx (log :
if dr 1.
where Ymin 81 ................................ < VK S Vmax is an equidistant grid of points on y with step size h, and (Yja) K is an i.i.d. sample drawn from P. (P A = a)). 13=1
Training. To train both components in our INFs, we make use of a two-step training procedure. Specifically, we first fit the nuisance parameters with the teacher flow. Then, we freeze the parameters of the teacher flow and fit the student flow. We additionally employ exponential moving average (EMA) of the student parameters with smoothing hyperparameter to stabilize the training for small minibatch sizes (Polyak& & Juditsky 1992). We show the full algorithm in Appendix and further implementation details in Appendix E
Inference time. One main advantage of our teacher-student model that the student flow has constant inference time (e.g.. during the evaluation phase). Hence, contrary to state-of-the-ae of-the-art baselines, the inference of our INFs do not depend on the dimensionality of covariates and the size of the training data. This is a major advantage over remi-parametric plug-in estimators. For a detailed runtime comparison, we refer to Appendix K As such, the student flow allows our method to scale well to large datasets as in medicine (Johnson al. 2016).
5.2 ONE-STEP BIAS CORRECTION
To provide an efficient estimation for the parameters of the student flow, we augment the cross-entropy loss (Eq. (9) with one-step bias correction. To evaluate the bias correction term, we need to compute an approximation of the conditional cross-entropy loss (Locr(X; Ba)) as in Eq. (9.. We thus compute
1=1 log ( 8.)) P(Y = ys X. X
df = 1,
Lece(X: v=
gog
Px (log x.; B.)):
dy >
where Ymui S 81 < yK S Fmax is an equidistant grid of points on y with step size h, and (Y, K J==1 an i.i.d. sample drawn from P(Y X, A a). Finally, we obtain the loss of the
'This is standard approach in neural conditional density estimation (Bihhop 1994 Kingma & Welling 2013 *One can use single normalizing flow with hypernewoor for categorical treatments.
7

---

Under review as a conference paper at ICLR 2023
Figure 3: Results for synthetic data based on the SCM from Figure train-test splits. Some runs for MDNs resulted with the log-prob out
Reported: mean over ten-fold
and, thus, are not shown.
student flow (C,), which is now suitable for A-IPTW estimation from Eq. (9. We thus yield
L.(B.) = (.. + Pn 1(A a) ( (- lYY Lcce(X:B.))]
(10)
Ra
EXPERIMENTS
6.1 OVERVIEW
To show the effectiveness of our INFs, we use established Seemi-ssynthetic datasets that have been previously used for treatment effect estimation (Shi et al. 2019 Curth & van der Schaar (021). The benefit of bemii-syythetii datasets that both factual and counterfactual outcomes are available (i.e., Y! and ). Therefore, we can obtain sample from the ground-truth interventional distribution, i.e Yaa: = (A, = a) Y + (A, + a) Y which we can then use for IDE benchmarking.
Evaluation metric. We use the average log- proobabiiiyy as our standard metric for comparing density estimators. is given by log-prob,
Σx= log P(Y[a) Y[a))), where higher values indicate better fit. The maximum value of the average log-probability is upper-bounded by the entropy, which, general, different for each potential outcome. Therefore, we separately report the results for each potential outcome. Of note, the tog-prbbability is equivalent to the empirical KL-divergence.
Baselines. We use state-of-the-art IDE baselines (see Sec. 2.T): (1) an extended TARNet (TARNet")
(Salit et al. 2017) estimating the mean of conditional homoscedastic normal distribution: (2) mix- ture density networks (MDNs) [Bishop| 1994 St (3) conditional normalizing flow (CNF) Trippe &
Turnrr 2018); (4) kernel density estimation (KDE) (Kim al. 2018): and (5) distributional kernel mean embeddings (DKME) (Muandet et al.) (021). TARNet", MDNs, and CNF are semi-pane eemi-parametric
plug-in estimators (see Eq. 2). Importantly, KDE and DKME do not guarantee a proper density estimation (unlike our INFs). We thus performed an additional re-normalization and negative values clipping, so that we can use the average log- proobabbbiyl as an evaluation metric. Details on the baselines are in Appendix F and hyperparameter tuning is reported in Appendix G
Ablation studies. We compare three variants of our INFs: (1) INFs (main): Our INFs introduced above using A-IPTW estimation. (2) INFs w/o stud flow: A simplified variant which uses only the conditional density estimation from the teacher flow as semi parammerreet plug-in estimator, and thus without student flow. This variant is identical to the CNF baseline. (3) INFs w/o bias corr: We use the adjjutee param-ra estimator, where the student flow only uses the cross-entropy loss from Eq. (9) but without one-step bias correction. The ablations have the same hyperparameters as our main method for better comparability.
6.2 RESULTS
Synthetic data. We generate synthetic data using the SCM (dx = 1) from Fig. [ Here, we vary the covariate shift b, which controls the overlap between the treated and non-treated population. Notably, low values b correspond to case, where both populations are similar or the same, while high values of b result in the violation of the positivity assumption. Further details the synthetic dataset are provided in Appendix H Fig. 3 shows the results. Our INFs achieve clear performance improvements over the baselines, especially for larger 6. Moreover, the ablation studies confirm that our proposed deep learning architecture with one-step bias correction is superior. In Appendix []] we additionally provide two-dimensional benchmark, where INFs prove their effectiveness.
IHDP dataset. The Infant Health and Development Program (IHDP) (Hill 2011) is a semi-synthetic dataset with two synthetic potential outcomes generated from real-world med-
'MDNs were previously used to estimate the conditional distribution of outcome for quantifying the ignorance regions of ITE estimation (Jssson al. 2021 2022) However, this different from our IDE task.
8

---

Under review as a conference paper at ICLR 2023
ical covariates (n 747, dx fold train/test splits (90%/10%) and Results are in Table 2 TARNet* is known to entail ground-truth condi- tional distribution model and should thus not be interpreted as a baseline but as an upper performance bound. Our INFs reach an equally good per- formance and, importantly, outper- form all the other baselines for both potential outcomes. The ablation study again confirms that our main INFs are superior over the other variants without the student flow we repeat the evaluation using the empirical Wasserstein
25, see details in Appendix H Here, we used ten- perform hyperparameter tuning based on the first split.
Table 2: Results for IHDP Mataset.Reportee mean + sd.
peg.
loggrah
keg prift..
leggrmh
TARNG (4 grousd outh fo (HDP)
11919 a 0.011
0.928
0.88%
-0.35 a: 8.81
R.634 : R.875
MDNS
1022 = 0.02
-8842
-L,879
2: 0.068
-0 0044 00077
CNF
Bow)
0,943 0.002
-0.900
-0.579
0.061
0:091
XDE
Kinstaa [onE]]
0L042 2: 0.010
0.068
0.060
0,200 +: 0.044
0.708 0.096
DKME Pararne 2021
0.540 0.010
00002
0.665 0.015
0066
INFS wh has cor
-1.52 4: 2-0.013
00.07 4: 0.008
-0670 +: 0006
Higher beter (hear hold scood bad underleed))
and without bias correction. In Appendix J distance with similar findings.
Table 3: Results for ACIC 2016 and ACIC 2018. Reported: % of runs with the best performance.
ACIC 2016 & 2018 datasets. ACIC 2016 & 2018 provide collection of semi-synthetic datasets with various data gettaggnerating mecha- nisms (Dorie al. 2019 (Shimoni al. (018) (see details in Appendix H. We select 15 random datasets from ACIC 2016 and 24 random datasets (4 of each of 6 sizes) from ACIC 2018. We per- form five random train/test splits (80%/20%) for each dataset, tune hyperparameters on the first split and evaluate the average in- and out-sample log- prepoabbbillty on every split. Table providee the performance comparison. Again, our INF's have a clear performance improvement over both baselines and other model variants. Compared to MDNs as the second-best method, our INFs scale much better in terms of runtime, especially for large sample sizes (see Appendix [))
ACIC 2016 (77 dazsett)
ACIC 2018 (24 datsetts
St heep
14 &bbo
TARNE MDN;
4.58%
K.3155
10.425
31.42%
21.67%
19.17%
flos]
KDB Kinestal F011
1.041%
1.1758
10.2256
9.58%
DRME Patrarte ar 202
0.395
0.78%
8.75%
10.33%
INFs 200
6.74%
8.53%
5,0%%
2,003
INF' (malt)
41.04%
34.00%
30.42%
28.33%
Highe beller (boat bol))
HC-MNIST dataset. Hidden confounding MNIST dataset is semi-synthetic dataset (Jes- son et al. 2021), constructed on top of the canon- ical image dataset of handwritten digits (MNIST)
Table 4: Results for HC-MNIST. Reported: mean
- sd. over ten random train-test splits.
log probus
grob
(LeCun 1998). To satisfy the exchangeability assumption, we add hidden confounder to the set of all covariates, i. i.e., 28x28 images (dxx = 784 + 1). For dataset details, see Appendix H For our experiments, we use only the train subset of the original MNIST (n 42, 000). Here we used ten random train/test splits (80% and tune hyperparameters on the first split. Table 4 shows the results of the experiments. Note that the non- and cami-paraeetric baselines suffer from scalability issues and were thus excluded. Further, our INFs outperform the variant without a bias correction, i. e.e., the only other available baseline.
INFs with Nos cor INF: (mnza) 1.34 0.00
-1.3 a 0.18
-143 u1s
-1.40 0.17
-1.00 0.37
1.34 = 001
1.33 4 0,00
1.33 0.01
Higher botee (beal bold)
fablee[]shows
Scalability. Experiments with ACIC 2018 and HC-MNIST datasets showed high effectiveness of our INFs for datasets with large sample sizes (n >> 25, 000) and with high-dimensional covariates (dx > 100). We provide runtime comparison in Appendix K For HC-MNIST, non- and semi- parametric methods even become completely impractical due to memory and time constraints. Importantly, this is a major advantage of our fully- pallmmparamttric IDE estimator, INFs, over semi- paramentric plug-in etimators and other baselines.
Case study. We performed case study using data from California's 's tobacco control program to estimate its effect on tobacco sales. Previously, the evidence was primarily based on point estimates without information on the interventional density (Abadie al. 2010 Our INFs suggest that the program led to a large reduction in tobacco sales..
Discussion. Interestingly, both components are important for the final performance (see our ablation studies). First, the teacher flow with the help of noise regularization performs consistent estimation of the nuisance parameters. Second, the student flow uses estimated nuisance parameters to solve the optimization objective. The student flow is not redundant but crucial for computational performance. While simple NFs have a similar estimation performance in terms of goodness-or offffe,,--o-fit,, only our INFs have constant inference time (e.g., during the evaluation phase regardless of the data size). This is a major advantage of parametric treatment effect estimators over semi- parameerritiiicc plug-in estimators.
9

---

Under review as a conference paper at ICLR 2023
REFERENCES
Alberto Abadie, Alexis Diamond, and Jens Hainmueller. Synthetic control methods for comparative case studies: Estimating the effect of California's 's tobacco control program. Journal of the American statistical Association, 105(490):493-500 2010.
Ahmed M. Alaa and Mihacla van der Schaar. Bayesian inference of individualized treatment effects using multi-task Gaussian processes, In Advances in Neural Information Processing Systems, 2017.
Ahmed M. Alaa and Mihaela van der Schaar. Bayesian nonparametric causal inference: Information rates and learning algorithms. IEEE Journal of Selected Topics in Signal Processing, 12(5): 1031-1046, 2018.
Sourabh Balgi, Jose M. Pena, and Adel Daoud. Personalized public policy analysis in social sciences using causal-graphical normalizing flows. In Association for the Advancement of Artificial Intelligence, 2022.
Elias Bareinboim, Juan Correa, Duligur lheling. and Thomas Icard. On Pearl's hierarchy and the foundations of causal inference. In Probabilistic and Causal Inference: The Works of Judea Pearl, pp. 507-556. Association for Computing Machinery. 2022.
Alexis Bellot and Mihaela van der Schaar. Policy analysis using synthetic controls in continuous-time. In International Conference on Machine Learning, 2021.
loana Bica, James Jordon, and Mihaela van der Schaar. Estimating the effects of continuous-valued interventions using generative adversarial networks. In Advances in Neural Information Processing Systems, 2020.
Ioana Bica, Ahmed M. Alaa, Craig Lambert, and Mihaela der Schaar. From real-world patient data to individualized treatment effects using machine learning: current and future methods to address underlying challenges. Clinical Pharmacology Therapeutics, 109:87-100, 2021.
Peter J. Bickel, Chris A.J. Klaassen, Ya' acov Ritov, and Jon A. Wellner. Efficient and adaptive estimation for semiparametric models, volume 4. Springer New York, NY, 1993.
Christopher M. Bishop. Mixture density networks. 1994.
Philippe Brouillard, Sébastien Lachapelle, Alexandre Lacoste, Simon Lacoste-Julien, and Alexan- dre Drouin. Differentiable causal discovery from interventional data. Advances in Neural Information Processing Systems, 2020.
Victor Chernozhukov, Denis Chetverikov, Mert Demirer, Esther Duflo, Christian Hansen, Whit- ney Newey, and James Robins. Double/debiased machine learning for treatment and structural parameters. The Econometrics Journal, 21(1):C1-C68 2018.
Alicia Curth and Mihaela van der Schaar. Nonparametric estimation of heterogeneous treatment effects: From theory to learning algorithms. In International Conference Artificial Intelligence and Statistics, 2021.
I Laurent Dinh, David Krueger, and Yoshua Bengio. NICE: Non-linear independent components
estimation. arXiv preprint arXiv:1410.8516, 2014.
Laurent Dinh, Jascha Sohl-Dickstein, and Samy Bengio. Density estimation using real NVP. In International Conference on Learning Representations, 2017.
Hadi Mohaghegh Dolatabadi, Sarah Erfani, and Christopher Leckie. Invertible generative modeling using linear rational splines. In International Conference on Artificial Intelligence and Statistics, 2020.
Vincent Dorie, Jennifer Hill, Uri Shalit, Mare Scott, and Dan Cervone. Automated versus do-it- yourself methods for causal inference: Lessons learned from a data analysis competition. Statistical Science, 34(1):43 68, 2019.
10

---

Under review as a conference paper at ICLR 2023
Conor Durkan, Artur Bekasov, lain Murray, and George Papamakarios. Neural spline flows. In Advances in Neural Information Processing Systems, 2019.
Damien Garreau, Wittawat Jitkrittum, and Motonobu Kanagawa. Large sample analysis of the median heuristic. arXiv preprint arXiv:1707.07266 2017.
Christian Gische and Manuel C. Voelkle. Beyond the mean: A flexible framework for studying causal effects using linear models. Psychometrika, 87:868-901, 2021.
Steffen Grünewälder, Guy Lever, Arthur Gretton, Luca Baldassarre, Sam Patterson, and Massimiliano Pontil. Conditional mean embeddings as regressors. In International Conference on Machine Learning, 2012.
David Ha, Andrew M. Dai, and Quoe V. Lc. HyperNetworks. In International Conference on Learning Representations, 2017.
Negar Hassanpour and Russell Greiner. Learning disentangled representations for counterfactual regression. In International Conference on Learning Representations, 2019.
Tobias Hatt and Stefan Feuerriegel. Estimating average treatment effects via orthogonal regularization. In CIKM, 2021.
Jennifer L. Hill. Bayesian nonparametric modeling for causal inference. Journal of Computational and Graphical Statistics, 20(1):217-240, 2011.
Chin-Wei Huang, David Krueger, Alexandre Lacoste, and Aaron Courville. Neural autoregressive flows. In International Conference on Machine Learning, 2018.
Paul Hünermund, Jermain Christopher Kaminski, and Carla Schmitt. Causal machine learning and business decision making. In Academy of Management Proceedings, 2021.
Maximilian llse, Patrick Forré, Max Welling, and Joris M. Mooij. Combining interventional and observational data using causal reductions. arXiv preprint arXiv:2103.04786, 2021.
Andrew Jesson, Sören Mindermann. Uri Shalit, and Yarin Gal. Identifying causal-effect inference failure with uncertainty-aware models. In Advances in Neural Information Processing Systems, 2020.
Andrew Jesson, Sören Mindermann, Yarin Gal, and Uri Shalit. Quantifying ignorance in individual- level causal-effect estimates under hidden confounding. In International Conference on Machine Learning, 2021.
Andrew Jesson, Alyson Douglas, Peter Manshausen, Nicolai Meinshausen, Philip Stier, Yarin Gal, and Uri Shalit. Scalable sensitivity and uncertainty analysis for causal-effect estimates of continuous-valued interventions. arXiv preprint arXiv:2204. 2022.
Fredrik Johansson, Uri Shalit, and David Sontag. Learning representations for counterfactual inference. In International Conference on Machine Learning, 2016.
Alistair W.W. Johnson, Tom J. Pollard, Lu Shen, H. Lehman Li-Wei, Mengling Feng, Mohammad Ghassemi, Benjamin Moody, Peter Szolovits, Leo Anthony Celi, and Roger G. Mark. MIMIC-III A freely accessible critical care database. Scientific Data, 3rt):ere 160035. 2016.
Edward H. Kennedy. Optimal doubly robust estimation of heterogeneous causal effects. arXiv preprint arXiv:2004.1444 2020.
Edward H. Kennedy. Semiparametric doubly robust targeted double machine learning: A review.
arXiv preprint arXiv:2203.0644 06469, 2022.
Edward H. Kennedy, Sivaraman Balakrishnan, and Larry Wasserman. Semiparametric counterfactual density estimation. arXiv preprint arXiv:2 20..102..220440 12034, 2021.
Ilyes Khemakhem, Ricardo Monti, Robert Leech, and Aapo Hyvarinen. Causal autoregressive flows. In International Conference on Artificial Intelligence and Statistics, 2021.

---

Under review as a conference paper at ICLR 2023
Kwangho Kim, Jisu Kim, and Edward H. Kennedy. Causal effects based on distributional distances. arXiv preprint arXiv: 18X.0238 2018.
Diederik P. Kingma and Jimmy Ba. Adam: method for stochastic optimization. In International Conference on Learning Representations, 2015.
Diederik P. Kingma and Max Welling, Auto-encoding variational Bayes. arXiv preprint arXiv:1312.6114 14, 2013.
Sören R Künzel, Jasject S. Sekhon, Peter J. Bickel, and Bin Yu. Metalearners for estimating heterogeneous treatment effects using machine learning. Proceedings of the National Academy of Sciences, 116(10):415-44665, 2019.
Yann LeCun. The MNIST database of handwritten digits. http:///ann.eccc..oom////////
1998.
M.F. MacDorman and J.O. Atkinson. Infant mortality statistics from the linked birth/infant death data set-1995 period data. Monthly Vital Statistics Report, 46(6 Suppl 2):1-22, 1998.
Krikamol Muandet, Motonobu Kanagawa, Sorawit Saengkyongam, and Sanparith Marukatat. Coun- terfactual mean embeddings. Journal of Machine Learning Research, 22:1-71, 2021.
Jens Müller, Robert Schmier, Lynton Ardizzone, Carsten Rother, and Ullrich Köthe. Learning robust models using the principle of independent causal mechanisms. In DAGM German Conference on Pattern Recognition, 2021.
Lizhen Nie, Mao Ye, Qiang Liu, and Dan Nicolae. VCNet and functional targeted regularization for learning causal effects of continuous treatments. In International Conference on Learning Representations, 2021.
Kenneth R. Niswander. The collaborative perinatal study of the National Institute of Neurological Diseases and Stroke. The Woman and Their Pregnancies, 1972.
Judea Pearl. Causality. Cambridge university press, 2009.
Boris T. Polyak and Anatoli B. Juditsky. Acceleration of stochastic approximation by averaging. SIAM Journal on Control and Optimization, 30(4):838-855, 1992.
Danilo Rezende and Shakir Mohamed. Variational inference with normalizing flows. In International
Conference on Machine Learning, 2015.
James M. Robins. Robust estimation in sequentially ignorable missing data and causal inference models. In Proceedings of the American Statistical Association, 2000.
James M. Robins and Andrea Rotnitzky. Comment on the Bickel and Kwon article, "Inference for semiparametrie models: Some questions and an answer". Statistica Sinica, 1(44::200-93, 2001.
Jonas Rothfuss, Fabio Ferreira, Simon Bochm, Simon Walther, Maxim Urich, Tamim Asfour, and Andreas Krause. Noise regularization for conditional density estimation. arXiv preprint arXiv: 09982 2019.
Donald B. Rubin. Estimating causal effects of treatments in randomized and nonrandomized studies. Journal of Educational Psychology, 66(5):688, 1974.
Patrick Schwab, Lorenz Linhardt, Stefan Bauer, Joachim M. Buhmann, and Walter Karlen. Learning counterfactual representations for estimating individual dose-response curves. In AAAI Conference on Artificial Intelligence, 2020.
Uri Shalit, Fredrik D. Johansson, and David Sontag. Estimating individual treatment effect: general- ization bounds and algorithms. In International Conference on Machine Learning, 2017.
Claudia Shi, David Blei, and Victor Veitch. Adapting neural networks for the estimation of treatment effects. In Advances in Neural Information Processing Systems, 2019.
Yishai Shimoni, Chen Yanover, Ehud Karavani, and Yaara Goldschmnidt. Benchmarking framework for performance-evallation of causal inference analysis. arXiv preprint / 002. 2018.
12

---

Under review as a conference paper at ICLR 2023
David Spiegelhalter. Risk and uncertainty communication. Annual Review of Statistics and Its Application, 4(1):31-60, 2017.
Esteban G. Tabak and Eric Vanden-Eijnden. Density estimation by dual ascent of the log-likelihood. Communications Mathematical Sciences, 8(1):217-233, 2010.
D. Michael Titterington, Smith Afm, Adrian F.M. Smith, U.E. Makov, et Statistical analysis of finite mixture distributions, volume 198. John Wiley & Sons Incorporated, 1985.
Brian L. Trippe and Richard E. Tumnr. Conditional density estimation with Bayesian normalising
flows. arXiv preprint arXiv:1802.04908 04908, 2018.
Anne Marthe van der Bles, Sander van der Linden, Alexandra L.J. Freeman, James Mitchell, Ana B. Galvao, Lisa Zaval, and David Spiegelhalter. Communicating uncertainty about facts, numbers and science. Royal Society Open Science, 6(5):181870, 2019.
Mark J. van der Laan and James M. Robins. Unified methods for censored longitudinal data and causality, volume Springer New York, NY, 2003.
Mark J. van der Laan and Daniel Rubin. Targeted maximum likelihood learning. The International Journal of Biostatistics, 2rt)cetticle 11, 2006.
Mark J. der Laan, Sherri Rose, et al. Targeted learning: Causal inference for observational and experimental data, volume 10. Springer, 2011.
Matthew J. Vowels, Sina Akbari, Necati Cihan Camgoz, and Richard Bowden. A free lunch with influence functions? Improving neural network estimates with concepts from semiparametric statistics. arXiv preprint arXiv:2202.09096. 2022.
Stefan Wager and Susan Athey. Estimation and inference of heterogencous treatment effects using random forests. Journal of the American Statistical Association, 113(523):1225 1228-124222 1242, 2018.
Rongguang Wang, Pratik Chaudhari, and Christos Davatzikos. Harmonization with flow-based causal inference. In International Conference Medical Image Computing and Computer-Assisted Intervention, 2021.
Antoine Wehenkel and Gilles Louppe. Graphical normalizing flows. In International Conference on Artificial Intelligence and Statistics, 2021.
Jeremy Yang, Dean Eckles, Paramveer Dhillon, and Sinan Aral. Targeting for long-term outcomes. arXiv preprint ar axXiv:2010.15835 2020.
Yao Zhang, Alexis Bellot, and Mihaela der Schaar. Learning overlapping representations for the estimation of individualized treatment effects. In International Conference on Artificial Intelligence and Statistics, 2020.
13

---

Under review as a conference paper at ICLR 2023
A RELATED WORK: NORMALIZING FLOWS FOR CAUSAL INFERENCE
NFs have been used the wider area of causal inference, yet in vastly different tasks than ours. Examples include, e. g., robust prediction by employing causal mechanisms (Müller et al. 2021 )
combining interventional and observational datasets (Ilse et al. 2021); and causal discovery (Brouil- lard al. 2020). Further, several works aim to model Bayesian networks structural causal models (SCMs) with known or unknown causal diagrams. For example, NFs were used as probabilistic model for Bayesian networks aimed causal discovery, as well as downstream interventional and counterfactual inference (Khemakhem et al. 2021 Wang et al., 2021 Wehenkel & Louppe (021). Balgi al. (2022) build upon a temporal SCM with exogenous noise, where NFs are used for inter- ventional and counterfactual queries. Importantly, all the aforementioned methods assume continuous variables in SCMs and independence of exogenous noise Hence, these methods are not applicable in our case, which considers eemi-Markovian SCMs and which is thus different inference task In sum, NFs have not yet been adapted to IDE, which is our novelty.
TThis commonly known as causal Markov condition.
Thhi is stated in our identifiability assumptions: there is limitation the exogenous noise independence hetween outcome and covariates. Hence, setting more general.
14

---

Under review as a conference paper at ICLR 2023
B BACKGROUND MATERIALS
B.I NORMALIZING FLOWS
Normalizing flows (NFs) TTabak & Vanden-Eljnden 2010 Rezende & Mohamed, 2013) are flexible probabilistic models with tractable density. A normalizing flow describes the change of the density of a continuous random variable after applying a sequence of invertible transformations. Given
random variable Z with some known density PZZ +), e.g., normal or uniform, we define transformed variable
( 2 P(Z),
(1))
where t(.) X' denotes an invertible forward transformation with inverse (-) - Z. Importantly, the transformation is defined between spaces of same dimensionality dz = dx. To find a distribution of X, we can apply the multivariate change of variables formula
dz
dl-1
P(X =) P(Z =) det
= P(Z = (x)) det
dx (z)/
(12)
dx
where det ax (x) is the Jacobian determinant the inverse transformation t-) Then, using the inverse function theorem, we obtain
dt
dt
(13)
dx
dz
so that the Jacobian of the inverse transformation can be substituted with the inverse Jacobian of forward transformation. Using the properties of the determinant, Eq. 12 can be simplified to
dt
P(X = 2) =) = P( = ==(==)-+(=))
det dz (t- (t)
(14)
The name normalizing comes from the fact that any regular continuous distribution X can be transformed to a normal Z with a specific t- (-).
We can construct arbitrarily complex densities by applying a composition of K transformations t1, .2. tx: .x:
(x ( (xx(x- -... lxo..... o...... (Zo).
(15)
where Zo is called base distribution. One calls this chain of transformations flow. Finally, the density of can be recursively found as
dtk (=K-1)
det dtp
P(ZK =k) P(ZZ-1 = == det dzk-1
P(Zo =0) II
(=k-1)
dZ
(16)
where 2g, 21...., 2K are found via Eq. L5 Consequently, we now can directly evaluate the log-likehood of observation X; Zki and, with a proper parametrization of transformations, back-propagate trough it. Examples of simple transformations include affine, planar, and radial (Rezende & Mohamed, 2015).
B.2 CAUSAL MODEL AND INDENTIFICATION
In this section, we provide a brief background on the underlying causal model in this paper, using both the potential outcomes and the structural causal model framework. These frameworks equivalent in the sense that they both allow for identification of the interventional density and yield the same statistical estimand.
Potential outcomes framework. The observed variables our model are covariates X E C Rdx a treatment E 1), and dd-dennnsio continuous outcome Y E y C R the main paper, we used the potential outcomes framework (Rubin 1974) to define the causal estimates. In particular, defined Y [a] the potential outcome after intervening on treatment by setting By imposing Assumptions (1)(3) in Section this allows to define the interventional density (our causal estimand) via
P(Y =)
PY y X A= P(X =) JEX
dr
E xxP(X)
(P(Y=) =y X, d)). (17)
15

---

Under review as a conference paper at ICLR 2023
SCM framework. Equivalently to the potential outcomes framework,
also define the
interventional density within the structural causal model (SCM) framework (Pearl 2009 Bareinboim et al. 2022). More precisely, we can define 'smmiMaarovvinn SCM by introducing independent exogenous latent variables UA ~ P(UA) and UXy ~ IP(Uxy): ); and the functional assignments fx (Uxxx,), A == fx(X, UA), and Y Ty(X,Uxx) Here, X, A and Y observed endogenous variables, satisfying Assumptions (1)-(), We show a corresponding causal graph in Fig.
X
A
Figure 4: Causal graph corresponding to the potential outcome framework assumptions.
Interventions vs. counterfactuals. We follow Pearl's hierarchy on causal inference (Bareinboim al. 2022) and distinguish the interventional and the counterfactual distribution. In SCM language, we can use Pearl's do-notation do(A) = to denote an intervention on the treatment A. This corresponds setting A in graph G, where all arrows from parent nodes of A A are removed.
We can then define the potential outcome [[a] via its interventional density P(Y y) P(Y = y do (A a))
(18)
and obtain the identification result from Eq.
In contrast, counterfactual queries aim answer individualized questions "what would have happened we had used a different treatment given already treated or untreated population". We can then define the counterfactual density
P(([a') = -
(19)
for some different treatment a' * a. This is the distributional equivalent to the average treatment effect of the treated (ATT). However, most of the treatment effect estimation literature focuses on interventional causal estimands (such as the ATE), Our paper is therefore in line with previous work. We acknowledge that other papers oftentimes call the interventional distribution counterfactual distributions for simplicity.
Comparison to other identification strategies. For the identification of the interventional density, we mainly rely on the three main assumptions positivity. consistency, and exchangcability (or, cquivalently, on the back-door adjustment from Eq. (I7)). This is common setup treatment effect estimation (van der Laan Rubin 2006 Shalit et al. 2017] Wager & Athey 2018). More complex adjustment rules (e.g., front-door adjustment, adjustment for napkin graph) have the following limitations: (1) they require more unusual, complex assumptions which are often violated in practice; and (2) they require a complex efficient estimation theory Vowels al. (002). Nevertheless, this could be an interesting direction for future research.
B.3 EFFICIENCY THEORY AND INFLUENCE FUNCTIONS
In this section, we give a brief background on semi parametric efficiency theory and influence functions. Our background builds upon Kennedy et al. (022), and we thus refer to it for mathematical details and further explanations.
Let us consider parimerraettric statistical model (IP E P), where P is family of probability measures. We are interested in estimating functional s: P - TR. sb sufficiently smooth, admits the so-called von Mises- distributional Taylor expansion
(1P) (PP) ( ((t. (t) R2 (PP
/
(20)
16

---

Under review as a conference paper at ICLR 2023
where (I, P) is second-order remainder term and ott, P) is the so-called efficient influence function of n, satisfying J ((t, P))P(t) = and ott d(P(t) < oc.
The efficient influence function -(:--) plays an important role in the theory of efficient semi camim-arttttii estimation. Under certain assumptions, can be shown that, for any sequence of estimators &... holds true that
inf lim sup n Eq (v - ((2))) 2 var (T(,,))),
(21)
850 n-20 Tx
where TV denotes total variation. Hence, e characterizes the best possible variance an estimator can achieve (in a local min-max sense).
Let now P be estimator of IP and (P) the so-called plug-in in estimator of (P)). The Mises expansion from Eq. (20) implies that (IP) yields first-order plug-in in bias because
eb (P) (P) - f ott, PP) dI(()) + PP)
(22)
due to that [(l,, P) (I) = 0. simple way to correct for the plug-in bias to estimate the bias term from the right-hand side of Eq. (22) and add it to the plug-in estimator via
GA-IPTW = (f) + (o(,
(23)
Under certain assumptions it, can be shown that the bias-corrected estimator @A-IPPTW is asymptotically normal with mean zero and variance var (T P)). Hence, by Eq. (2I), @A-IPTWW is (asymptotically) efficient in the sense that it consistent with the best possible variance.
Application to interventional density estimation: We now return to the specific statistical model in our paper, i.c., we aim interventional density estimation. In other words, the estimand (PP) we are interested in is the function
P( =-) = x~P(x) (P(Y = - (X,A
(24)
Given an initial estimator P(y = | X, = a) and the marginal empirical probability measure
Pu (-), the plug-in estimator becomes
P. PP
(25)
As described above, this estimator suffers from plug-in bias and is not efficient. However, one-step bias correction for our setting is not as simple due to the fact that the interventional density is
functional target estimand and, hence, infinite dimensional. As remedy, |Knnedy al. (022 proposes an elegant solution by introducing the finite dimensional projection parameter
3. = arg min KL (P(Y(a)) (((:B))).
(26)
which is equivalent to solving the moment condition
m(3,) = x-P(x) (B((();8)) X, |x, a)) -0
(27)
where T = P.) - v. log gY;; 3.). The advantage of this approach is that the moment m(3.) is a finite dimensional quantity, which means efficiency theory can be applied. The plug-in
estimator for the isment
m m
E Yo (P
T(Y"; (.))
(28)
Kennedy et al.
[2021) also derived the efficient influence function for the moment:
(T;
1 A = a)
(T- B(T (
-5 a)) + E(T x, --
X~P(X) (E(T) X, A-
a)).
n (X)
(29)
Hence, a bias-corrected estimator for the projection parameter can be obtained by solving
APTW (3.) Pl (B.) + PP (P.(('TT() (T.): _ 0.
(30)
Estimating the projection parameter via Eq. 30 requires solving a (potentially digh-dimension)
system of non-linear equations, which is often infeasible in practice. Hence, remedy, we propose
in this paper to reformulate Eq. (30) as an optimization problem which can be incorporated directly
into loss of neural network (see Appendix C
17

---

Under review as a conference paper at ICLR 2023
C BIAS-CORRECTED MOMENT CONDITION AS AN OPTIMIZATION TASK
We aim to transform the bias-corrected moment condition into an optimization objective:
A-PT (Ba) = m (8.) + ( P)} - 0.
(31)
We first note that the plug-in estimator of moment condition rn (B.) can be rewritten as
T(Y"; = / T(y: Pa) Pp, (p = x, A a)} dy
E vit (Ba) (
(32)
-Pa T(y: &.y
x, o) du)
= = ( - X | })). (33)
where the last equality follows from the definition of the conditional expectation. Notably, we see that the moment condition could be equivalently solved with either the conditional distribution, P(Y X, A = a), or with the functional regression, IE (TY;; X,A= a).
Let us unroll the bias correction term of Eq. 7
PT]( Pn
1 1(A a)
(T B x. x= E(T X, A A-A=00. a) - PP,((T X, X,A=
(X)
(34)
((A=aa) a)
Na (X) (T- X, ( + E(T X, A- a) } Pr (B X, a))}
=
(35)
where nuisance parameters are marked with red color. Here, the last term is in fact the plug-in in estimator of the moment condition, i.e., (B.). Therefore, we can simplify one-step bias corrected moment condition via
1 1(A (A)=) (T-B(T) B(T x.A= a)) B((T) X, aa at) (X)
A-PPT (Ba) P
(36)
T(Y": Ba) + PP (1AA= = a)
(T(:;
E x~
:(Y:Ba)
You7, A=))
(X)
(37)
where we use the conditional density estimator not an estimator for the functional regression. This allows to transform the A-IPTW moment condition into optimization objective (Eq. ( by taking antiderivative with respect to 3a.
18

---

Under review as a conference paper at ICLR 2023
D TWO-STEP TRAINING PROCEDURE
Our INFs are trained with a two-step procedure. The procedure is shown in Algorithm [ we use noise regularization as the main regularization technique for the teacher flow, and exponential moving average (EMA) for the student flow to stabilize training. A-IPTW estimation is also known become unstable finite sample setting (Shi al. 2019), so that inverse values of propensity score become too large. Thus, we manually discard observations with too small propensity score (fra < 0.05) from bias correction.
Recall that
Algorithm 1 Training procedure of INFs
Input: number of iterations Niter,ss Tlier,s minibatch sizes b,, ba; learning rates T) Taa intensities of the noise regularization o o; EMA smoothing 7: grid size K.
Init: parameters of the teacher flow: FC F
D. Fitting the teacher flow
for i 0 to Pitter do
= A, + minibatch of size by
#,, FC (X)
&, N(0, N(0 og); R R+ + Y
Y Ey
Noise regularization
0(X, =FC," (A,
P(Y A) normalizing flow with parameters 0(X, A)
ENLL f- 4 Po
Lx t (C)
L((P p {NNLL + aL
rC((+1) ..cc+t))
optimization step wrt. Le R.) with learning rate T
end for
Output: nuisance parameters: P((Y) A), Aa (X)
Init: parameters of the student flows: ( (((0) a.E
(000
Fitting the student flows
for i Vliter, do B = (X, Y + minibatch of size b.
for a E {0. do
( + hh === log (z0i (00) p (P = = X = a))
LOCE(X X; x;b(1)) - Σ j=1 log gu;; (1)
JP(Y = Vj X, 1 = a.
bias correction (3(")) p
(A=o&& (X)20.05)
log gY:: (1(0)) LOce (X;
(x(X)
(,(3)()) (B() +bias correction (3(5))
a(++1) + optimization step wrt. C, (,(())) with learning rate 7)
g(++l)
a.EMA < 70(11) a.EMA (1 7)
b EMA update
end for end for Output:
3a
A-IPTW
(
EME
19

---

Under review as a conference paper at ICLR 2023
E INFS IMPLEMENTATION DETAILS
Implementation. We implemented our INFs using Py PyTorch and Pyro. For both teacher and student flow, we employ neural spline flows (Durkan al., 2019) with standard normal (NV (0; 1)) as a base distribution. Neural spline flows construct an invertible transformation of the base distribution with the help of monotonic rational-quis -uaadratie splines. They are characterized by two main hyperparameters: number of knots Ramoss and a span of the transformation interval [-B; Bj. Namoos controls the smoothness estimated density and B affects the support of the transformation. our experiments, we heuristically set B = Vmax Vmin + 5.
For the teacher flow, we use fully- fully-connected subnetworks each with one hidden layer (with h 10 hidden units), and the dimensionality of representation is set to r 10.
Training. During training (see full algorithm in Appendix Dp. we adopt noise regularization (Roth-
fuss et al. 2019) and add an independent Gaussian noise &.r N(0, o2), ~ N(0, oz)) to the representation and output of the teacher flow, i. e., R = + f Ý Y + Ey. For faster learning, we approximate a full-sample average (. with minibatch average p (.}
for all the losses, where b is the minibatch size. We use stochastic gradient descent (SGD) for fitting the parameters of the teacher flow, and Adam optimizer (Kingma Ba 2015) for the student flow with learning rates Ti and Te. respectively. We fix the weighting hyperparameters of the loss to 1 and the EMA smoothing hyperparameter to = 0.995. The grid size for approximating the cross-entropy loss is set to K 100. Both Ymui and Vnax set to the empirical minimum and maximum of the train sub-sample.
Note that we would need sample splitting for training both flows to guarantee the asymptotic properties, e.e efficiency and double robustness (see Kennedy et al. 2021 Remark 5). Nevertheless, used all data for the both components and trained INFs with auxiliary regularization because sample splitting can affect the performance in settings with limited data. This is consistent with previous work on deep learning for efficient treatment effect estimation (Curth & der Schaar 2021.
Hyperparameter tuning. We perform extensive hyperparameter tuning only for the teacher flow. Hy- perparameters for tuning include, e. g.. number of knots of neural spline flows It the minibatch size b, the learning rate 7/+, and the intensities of the noise regularization o² o On the other hand, we discovered, that the student flow works well with the same plain set of hyperparameters in almost all the experiments. Those include the minibatch size b, = 64 and the learning rate = 0.005. The number of knots Rinass, is chosen hand for each dataset. Further details on hyperparameter tuning are provided in Appendix G
20

---

Under review as a conference paper at ICLR 2023
F BASELINES
In the following, we describe the baseline methods in detail. These are two naive semi-parametric plug-in estimators: mixture density networks (MDNs) (Bishop) 1994) and conditional normalizing flow (CNF) (Trippe & Turner, 2018). Further, we two sfate-of-the art IDE baselines: kernel density estimation (KDE) (Kim et al. 2018) and distributional kernel mean embeddings (DKME) (Muandet et al. 2021.
F.1
NAIVE SEMI- PLUG-IN ESTIMATORS
camimparaettric plug-in in estimators estimate the conditional outcome distribution and perform aver- aging over covariates during evaluation, as introduced in Eq. (2)
TARNet", MDNs, and CNF make use of hypernetworks (Ha et al. 2017), which take covariates X and treatment A as an input and output parameters, e.e., ((X,A) of the estimated conditional distribution (Y X, A). Hypernetwork architectures are considered to be state-of-the thate-t for neural conditional density estimation and can be found in, c. g., Gaussian mixtures (Bishop] 1194). variational autoencoders (Kingma & Welling 2013). and normalizing flows (Trippe & Turner 2018). For comparability. we use the same network structure of the teacher flow in our INFs as the hypernetwork for the conditional distribution parameters. This gives two fully-connected subnetworks stacked on each other, i. e. FC1 and FC2, introduced in Section 5.1 To regularize both conditional distribution estimators, we use noise regularization (Rothfuss et al. 2019).
The treatment-agnostic representation network (TARNet) (Shalit et al. 2017) was proposed to estimate nuisance parameters for ITE, i. c., conditional means of outcomes. To obtain density estimates as outputs, report results from extended variant which we refer to as TARNet". Specifically, we extended the original TARNet by modeling conditional outcome distribution as homoscedastic normal distribution. For this, add one unconditional parameter of standard deviation, a, so that the conditional density equals to
TARNet".
P-yY = X, = (XX,
(38)
where Ny; pt, )²) is a density of the normal distribution, and (xX,AA) A) is conditional mean of outcome. Notably, we do not use the two separate outcome heads (as in original TARNet) but only one, i. e., FC2. This is crucial to ensure a fair comparison with other plug-in estimators. We estimate the standard deviation using maximum-likclioood Note also that TARNet restricted to normal conditional outcome distributions and thus is not universal density estimator. In contrast to our INFs, TARNet* unable to capture heavy-tailed, multi-modal, and skewed distributions.
MDNs. Mixture density networks (Bishop 1994) are built on top of mixture of normal distributions, and can approximate any density arbitrarily well (Titterington et 1985). i.e.
no
?-y = X A)
Σ wy wx N(y (x,
(39)
where ne is a number of mixture components, wy 2 0, Σド j= w = are mixture weights, and N N Hj, of) is density of the normal distribution. In the case of MDNs, the hypernetwork outputs logits of mixture weights and parameters the normal distribution (i.e., mean and logarithm of the standard deviation), e., 0 = (logits(w;), Mjj log }. Here, the number of mixture components ne controls the smoothness the estimator and represents the main hyperparameter for tuning.
CNE. We implement conditional normalizing flow (Trippe & Turner 2018) with the help of neural spline flows (Durkan et al. 2019). Neural spline flows construct an invertible function parameterized by 0, e., f(::6): spline transforms the density of a base distribution on the interval [- B: B). Outside of the interval, f(.) equals to the identity function. This allows to perform flexible parametric density estimation with the help of the change of variables formula, i.e.,
R. which is monotonic rational-quadess
spline with Tiknoss knots. This
P(Y = X, A) N(5 0(X, A)):0,) dy
df (f- (g:
(40)
21

---

Under review as a conference paper at ICLR 2023
where (=)0) is the inverse transformation, and the density of standard normal distribution (y; 0, 1) serves as a base distribution. As already discussed in Appendix E B affects the support of transformation, and the number of knots Tiknots controls the smoothness of the estimator and represents the main hyperparameter for tuning.
F.2 KERNEL DENSITY ESTIMATION (KDE)
Kernel density estimation (KDE) is a cemi-parametric method for IDE (Kim et al. upon the idea of density functional, namely T, ha), to transform random variable Y into proper density via
2018 It builds
1 K
|| (IY-
yll2
1
|| vll2 2h a
TY;Y; ha)
exp
(41)
ha
ha
ha
where ((x) = 22n (xpp-=2²// /2) is a radial basis function (RBF) with treatment-specine smoothing
parameter ha called bandwidth, and |1-I12 is the L2-norm.
Robins & Rotnitzky (2001) proposed a semi-parametric plug-in estimator of interventional density
(Y y) P. P ((T, x, | x
(42)
where B XX) hY: A a) is a functional regression of X and A on TY;Y; ha). Kim et al. (0088) further extended this estimator to an efficient, A-IPTW-style semi-parametric
estimator
p.A-IPTW (Y|a]= =) Pmn
(1(A=0) {
(T,(Y;ha) - (x)x())) + Denu(X)3. (X)
(43)
where R (X) is an estimator of the propensity
The main challenge here is building a functional regression (X). Unfortunately. the work by Kim al. 2018) does not provide effective, practical solutions. Even more so, Eq.. (43) does not guarantee that the estimated density is proper, i. e., integrates to and positive, especially in small sample regime when the propensity score has extremely low values.
To estimate the nuisance parameters, namely, the propensity score and the functional regression, we the same network structure as for the teacher flow of our INFs (see Section 5.1). In this way, estimate the propensity score and perform functional regression with two joined, fully-connected subnetworks (i.e.. FC₁ and FC2). The first subnetwork. FC1, outputs a representation and estimates the propensity score. The second subnetwork, FC2, then takes the representation R and the treatment A, and performs an outcome regression: Ý ÉE( X, A). The functional expression, i. e.e Eq. (41 predicted via (X) (f(() X, A): h.). Although, this is biased estimator of Paay (X), ensures a proper normalization, i.e., Jy p dy 1
To fit FC, and FC2, we use the sum of mean-squared error (CMSE) and binary cross-entropy (Cx) losses via
LRIDE na) P., LASE aLr) with LMSE 1Ý YY); L. BCE(#A
(44)
where or is hyperparameter. In experiments, set and fit the nuisance parameters (i.e., f, and B(Y X, A)) using the Adam optimizer with fuerr 10000 iterations. Both learning rate n
and minibatch size b are subject to hyperparameter tuning.
We employ median heuristic (Garreau al., 2017) for choosing the bandwidth has i.e.
Median (I 11 A A=a). isi<jsn.
(45)
where ||:1|2 is the Ly-norm, and where Y, Y, are observations from the train subset. conditioned on = a. To address the numeric instability of the A-IPTW estimator, we discard observations with small propensity scores (#. (X) 0.05) from averaging in Eq. (43). similarly our INFs.
22

---

Under review as a conference paper at ICLR 2023
F.3 DKME
Distributional kernel mean embeddings (DKME) (Muandet et al. 2021 is a non-parametric plug-in estimator of interventional densities. This methods builds kernel mean embedding (KME), namely, |Y|X, for the conditional distribution P(Y X, A = a) via
E PY|X,A=a (V) Y~M(Y|X.A=a) (V,
(46)
where (a(-.) -) is a measurable positive definite kernel associated with reproducing kernel Hilbert space H. so that AAi provides a mapping from the space of conditional distributions to the space of functions H. (::-) properly normalized, then ( is in fact a conditional density estimator.
To estimate the KME of the conditional outcome distribution (conditional mean embedding). we use the i..dd. sample D = (X,, A,, Y) and split it into control and treated subsamples: D {X), yo ( y Then, Ay|X,A=a can estimated via
Ra
Av (1) = 2
w (x),
(47)
i=1
(wf(X)
na (X)) = (K" + + (okel)-k((X)) (
(48)
ks kxX) ((x R
(49)
where I E R is an identity matrix, E > is regularization hyperparameter, K E R X
is kernel matrix with elements K :) = ,X and k(.,-) is a second kemel representing
conditional dependencies between X and Y (Grünewälder et al. 2012).
Muandet et al. 2021) further developed KME for interventional distribution, i.e., [Y[a] and its empirical estimate, [Y[4
(Y([(1)- = E x-P(X)
(Y|X,A=a(1))
(50)
na
(y)a = P, TA 2 (y,
(51)
( = (K" + n==l) K 1m E R"
(52)
where K E Rn.xx is kernel matrix with elements KG k(;; X,), and 1m (1/n,...,1/o)
For our experiments, we choose both kernels, e., outcome kernel, (a(:)),, and conditional kernel, k(,)). -), to be RBF kernels with bandwidth parameters ha, and hg. respectively. Therefore, ( (v) represents a valid interventional density estimator. Nevertheless, due to small sample sizes, some B could be negative and the estimated density ends up having negative values.
We set the bandwidth of the outcome kernel, ha,i according to the median heuristic from Eq. (45). The bandwidth of the conditional kernel hą and the regularization hyperparameter are subjects to the hyperparameter tuning. Motivated by the interpretation of conditional mean embedding as kernel ridge regression (Grünewälder et al. 2012), use out-sample MSE of the ridge regression with parameters h2 and & as tuning criterion.
23

---

Under review as a conference paper at ICLR 2023
G HYPERPARAMETER TUNING
We performed hyperparameters tuning for all the baselines based on five-fold cross-validation using the train subset. For each baseline, we performed grid search with respect different tuning criteria, evaluated on the validation subsets. parameters, such as tuning criteria, number of training iterations, and optimizers. We aimed for fair comparison and thus kept the number of parameters, network structures, and grid size similar across models. For the sake of reproducibility, we make the chosen hyperparameters for all the experiments public (see YAML. files in our GitHut
shows grids for hyperparameter tuning and other
Importantly, facilitate the convergence of baseline methods, we additionally perform standard normalization of both factual and counterfactual outcomes for all the datasets.
Table 5: Hyperparameter tuning for baselines.
Model
Sub-model
Hyperparameter
Range / Value
Intensity noise regularization (a Intensity of noise regularization (ar)
0.01, 0.05, 0.0, 0.01, 0.05, 0.1
Learning rate (1))
0.001, 0.005
Minibatch size (b)
32, 64
TARNet
Tuning strategy Tuning criterion
random grid search with 50 runs
CxL1
Number of train iterations (7aer)
5000
Optimizer
SGD (momentum 0.9)
Number of mixture components (nc)
5, 10, 20
Intensity of noise regularization (o
0.0. 0.01, 0.05, 0.1
Intensity of noise regularization (o
0.0. 0.01, 0.05, 0.1
Learning rate (1)
0.001, 0.005
MDNs
Minihatch size (h) Tuning strategy Tuning criterion
32, 64
random grid search with runs
Lx11
Number of train iterations (Narr
5000
Optimizer
SGD (momentum = 0.9)
Learning rate
0.001, 0.005,
Minibatch size (6)
32, 64, grid search
KDE
Tuning strategy
Tuning criterion
LMSE + oxC,
Number of train iterations (Naer)
10000
Optimizer
Adam (betas= 0.999))
Kemel smoothness (d), 22h))
0.0001, 0.001, 0.01, 1, 20
Regularization parameter (e)
0.0001, 0.001, 0.01, 0.1, 1,
DKMR
Tuning strategy Tuning criterion
grid search MSE of ridge regression
Number of knots (Mlooos.)
5. 10, 20
Intensity of noise regularization (o
0.0. 0.01, 0.05, 0.1
Intensity of noise regularization (a.)
0.01, 0.05,
Learning rate (n)) Minibatch size (b.) Tuning strategy Tuning criterion Number of train iterations (Rars.) Optimizer
0.001, 0.005
Teacher flow
32, 64
random grid search with 50 runs
Lx11
INFs
5000
SGD (momentum 0.9)
Number of knots (Minsss.)
dataset specitic"
Learning rate (1)))
0.005
Minibatch size (b.) Tuning strategy Number of train iterations (Nass.a)
64
Student flow
w/o tuning
4000
Optimizer
Adam (betas=(0.9, 0.999)
Anons 5 (synthetic data), 10 (IHDP, HC-MNIST datasets), = Fikmots, (ACIC 2016 & 2018 datasets)
https:://ancnymouss 4open wtt.eeccee////oomooo/uu.//// interrl
24

---

Under review as a conference paper at ICLR 2023
H DATASET DETAILS
H.1 SYNTHETIC DATASET
We sample n 1000 observations from the SCM (Fig. D and use a ten-fold split for train/test samples (90% We separately perform hyperparameter tuning based on the first split for each baseline and each level b. We then report an average out-sample og-likelihood over ten folds.
H.2 IHDP DATASET
The IHDP dataset (Hill 2 uses a real-world dataset with 25 covariates (6 continuous, 19 binary) and one binary treatment, capturing aspects related to children and their mothers. Both treated and untreated, synthetic outcomes of IHDP are sampled from different conditional normal distributions. These distributions are homoscedastic (o² 1) but have substantially different conditional means. We used the setting "B" in (Hill 2011) with following SCM:
X (oald-orldd-);
A Real- Real-World(X))
(53)
Y (A w) (1 A) (exp((X + W))));1),
where 3, W, w are constant parameters of the simulation. For the further details, we refer to (Hill 2011
H.3 ACIC 2016 & 2018 DATASETS
Covariates of ACIC 2016 are taken from large study of developmental disorders (Niswander 1972.
and covariates of ACIC 2018 are derived from the linked birth and infant death data (MacDorman & Atkinson 1998). ACIC 2016 and ACIC 2018 differ in the number of true confounders, the varying level of overlap, and the form of conditional outcome distributions. ACIC 2016 has 77 different data- generating mechanisms with 100 equal- sized samples for each mechanism (n 4802, dx 82) ACIC 2018 provides 63 distinct mata--ennraatin mechanisms with around 40 non-equal-sized samples for each mechanism (n ranges from 1, 000 50, 000, dx 177). Notably, ACIC 2018 has constant ITE for most of the datasets, but heterogeneous propensity scores.
H.4 HC-MNIST
Jesson et al. (022)) introduced a complex high- high-dimensional, semi-synthetic dataset based on the MNIST image dataset LeCun (1998), namely HC-MNIST. This dataset maps high-dimensional images onto one-dimensional manifold, where potential outcomes depend a complex way on the average intensity of light and the label of an image. The treatment also uses this one-dimensional summary, Φ, together with an additional (hidden) synthetic confounder, U. This is described by the following SCM:
U (0.5(.
X MNIST-image())
(clip
UN, dc
1.4, 1.4
Min,
Max. Min. 1-4(-114)
a(0; I+) F. sigmmid(0.75 756
(54)
((: [*)
r"
Sggmoid0.750+tul
A Bern a((TT) A(o:T*)
Y N ((2A 1)Φ + (2A 1) (2(A 1)6) 2(2u 1)(1 + 0.55);;),
where is label of the digit from the sampled image X; K,N, is the average intensity of the sampled image; He and o. are the mean and standard deviation of the average intensities of the images with the label c and Min, factor influences the treatment assignment to a larger extent, i.c., the additional confounder or the one-dimensional summary. We set I cxp(1). For further details, we refer to (Jesson et al. 2021)
-2 + c, Max. - + 10 (c + 1). The parameter I" defines what
After one-hot-cneoding of categorical covariates.
25

---

Under review as a conference paper at ICLR 2023
For the experiments with HC-MNIST, we larger network size for our INFs (compared to other benchmarking experiments) to allow for more flexibility. We set the number of hidden units in fully-comnected subnetworks to h 30, and the dimensionality of representation increase the number of training iterations to Tfter, 15, 000 and fiter,s = 5000.
30. We also
Fig. shows both ground-truth interventional (P(Y(a))) and observational (P(YY A = a)) distri- butions together with our INFs A-IPTW estimator (IPA-IPTT (Y [a])).. Remarkably, the interventional distributions in HC-MNIST are multi- modal and differ lot from observational distributions.
207
0.6
7Y y 4-
FLY -
0.5
34-1 (V 10/00110 -
(y
0.4
0.2
0.1
00
Figure 5: Empirical ground-truth interventional and conditional distributions of the HC-MNIST synthetic outcome. We also plot our INFs density estimator, i.e., JAA-PTW (Yaa)).
26

---

Under review as a conference paper at ICLR 2023
1 SYNTHETIC TWO-DIMENSIONAL DATA
In the following, we benchmark our INFs for estimating an interventional density of the multidimen- sional outcome, dy 2.
Noisy moons synthetic data, We used a standard tvo-dimensional toy data generator, namely moons data 10 It draws samples from two interleaving half-cireles with different noise levels The noise level controls the level of the overlap between two half-circles (a higher corresponds to better overlap, and, thus, a satisfaction of the positivity assumption). We drawing of two-dimensional covariates, i. e.e dx = 2, and use an inclusion to the top or bottom semi-circle as a treatment. Finally, we generate the synthetic outcome by rotating the covariates by a random treatment-specifee angle, i. e., og and 01:
1000 observations
X, A Make-Moons(:-0
N(; V0.1",
(55)
0g- ] 5, 0g - S, Y == R(ai A ++00 ao (1 - A)) + =12,
[cos
sin a cos
where 12 (1, and R(a)
is a-angle rotation matrix. We set 0.75.
sin a
For the benchmarking with the noisy moons data, we increased the number of the training iterations for all the plug-in in methods (Naer 10000) and for our INFs, (Fier, 10000, Miter 5000). To model two-dimensional (conditional) density, we employed an auto-regressive extension of neural spline flows (Dolatabadi et al. 2020). We decreased the number of sampled points for approximating the cross-entropy, K = 70, to speed up the training, and set the number of knots for the student flow
7 5.
Results. Table 6 shows the results. Here, our INFs (main) scores second best in terms of in-sample performance, but, more importantly, best in out-sample performance. MDNs, although scoring the best with in in-same average log- pogorbbbbillty. do not generalize well. Finally, we again confirmed, that our INFs are superior over their ablations and other existing methods, e.g., KDE and DKME.
Table 6: Results for synthetic experiments using the noisy moons synthetic data. The performance is benchmarked using the empirical in in-sample / samlse average log probability for the two potential outcomes (i.e., a = and 1). Reported: mean standard deviation over ten-fold train-test splits.
log-prob.
log-prob...
log-prob,
log-probou
TARNet
2.907
0.121
-3.005 0.263
-2.781 0.092
-2955 : 0.222
MDNs
- 2.698
0.050
-2887 0.173
-2,83
1: 0.051
-2827 1 0.165
CNF INFs w/o stud flow]
- 2.767 0.087
-2935 0.239
-2807 + 0.162
-2.00 + 0.183
KDE Kim et al. 2018
2993 = 0.015
22916 1 0.052
-2898 1: 0.013
-2,01 1 0.049
Muandel at
2.822 0.016
22875 0.056
-2847
0.012
-2849 1 0.067
INF's w/o bias cot
- 2.787 = 0.057
-2.994 0.130
-27885
4: 0.048
-2.788 + 0.135
INFs (main)
2,764 0.030
2.766 - 0.102
-2.780
: 0.022
-2.85 + 0.134
Higher betler (best in bold, second best underlined)
att.r//aaeeee.eeeeeeae.ttattat dataseta.
make_moons.. html
27

---

Under review as a conference paper at ICLR 2023
J ADDITIONAL RESULTS
J.1 IHDP DATASET
Here, we provide additional results for the IHDP dataset with an alternative evaluation metric, that is, the empirical Wasserstein distance.
Evaluation metric. For dneedimensional outcomes, the Wasserstein distance between two distribu- tions can be simply expressed via quantile functions
WP(P1,PP2) = -(11
-F (2" 1/P
(56)
where FT' (4) and ' are quantile functions of 11 and P2. respectively. The Wasserstein distance is not upper-bounded and equals zero if and only if both distributions are the same. Here, we compute the empirieal Wasserstein distance, i. e., W based on empirical quantile functions. This requires two samples: one from the ground-truth interventional distribution and another from the estimated density. Therefore, methods which do not provide direct sampling (e. g., KDE and DKMH) cannot be used for evaluation.
Table 7 shows the results. Note that TARNet", e.e the plug-in with the ground-truth conditional density estimator for this specific dataset due to the fact how the data was constructed. Hence, we do not interpret TARNet* as baseline but rather interpret it as a bound for the best performance. We see that all baselines (MDNs, CNF, INFs w/o bias correction) are inferior by a large margin. In contrast, our INFs achieve a performance similar the bound. In particular, our INFs perform overall best: our INFs are superior over the two other na've plug-in estimators and the variant of INFs without bias correction. In sum, the results corroborate our findings from the main paper and add to the effectiveness of our INFs.
Table 7: Additional results for semi- mynihsynthetic experiments using the IHDP dataset. The performance benchmarked using the empirical in-sample out-sample Wasserstein distance (i.e., W and W W for the two potential outcomes (i.e., a 0 and a 1). Reported: mean + standard deviation over ten-fold train-test splits.
a
w
W11
Wi2
W11
TARNet" { ground-buth for IHDP)
0.048 + 0.014
0.131 + 0.040
0.046 :: 0.024
0.126 + 0.065
0.066 0.053
0.071
CNF [= INFs wlo stud flow]
0.118 = 0.048
0.192 + 0.069
0,111 = 0.087
0,146 = 0.082
wio
0.075
0.137
INFs (main)
0.040 - 0.009
0.132 + 0.051
0,100 + 0.037
0.117 = 0.055
Higher (best hold, second underlined)
J.2 ACIC 2016 & 2018 DATASETS
In the following, present detailed results of the experiments with ACIC 2016 and ACIC 2018 datasets. Fig. rrpports the median performance for the individual datasets in ACIC 2016, and Fig. 7 for ACIC 2018. In the latter, datasets are grouped by sample size. We also show the performance gain of our INF' (when INF' score better than the baselines). The percentage of the datasets with the positive performance gain for our INF' roughly correspond to the percentage reported in the Table For ACIC 2016, our INFs are the best method for 33 + 21 = 54 out of 77 = 154 potential outcomes of individual datasets (35%) with respect to out-of-sample average log- probbbbbbillty. For comparison, the second-best baseline (MDNs) are only best for 23 + 25 48 out of 77 potential outcomes of individual datasets (31% %), and thus inferior. For ACIC 2018, our INFs are the best method for + 9 out of 2 * 24 potential outcomes of individual datasets (35%). For comparison, the second-best baseline (also MDNs) are only best for 5 4 out of * 24 potential outcomes of individual datasets (19%), and thus again inferior. This thus provides consistent performance that our INFs are highly effective.
28

---

Under review as a conference paper at ICLR 2023
TAZINN
-06
-12
1400
1.00
-06
C,000
-1.
1.01
-16
-1000 SOEL
-12
111000
-10
C.C.CC
-12
1010
Figure 6: Detailed results for ACIC 2016. For each dataset, we perform five random train-test splits, tune the baselines on the first split, and evaluate the average in-sample out-sample log-probayility for each of the two potential outcomes separately. Shown: median over five runs and improvement of our INFs (main), when they score better than other baselines.
29

---

Under review as a conference paper at ICLR 2023
TARNA
XLDNN
CNF w/oaed Dan)
KTDB
DXME
DVFs
(mett)
16332
0
+0031
001
20001
10012
10 0.0
1022
3-
0.121
0035
190
2011
10 0003
EOLEEOONN 1800 www. aukottttoo
NU EA www oethhuttion
5500
7000
10000
15000
Dutect
Figure 7: Detailed results for ACIC 2018. sorted with respect sample sizes. For each dataset, perform five random train-test splits, tune the baselines on the first split, and evaluate the average iamsep sample prebbabbbyyyty for each of the two potential outcomes separately. Shown: median over five runs and improvement of our INFs (main), when they score better than other baselines.
30

---

Under review as a conference paper at ICLR 2023
K RUNTIME COMPARISON
INFs are marly-parametric model and, therefore, provide a decent speed up at the inference time. This is particularly important for scalability, that is, for datasets with large sample size and high- dimensional covariates. In Table S we report the total runtime of the baselines and the different variants of INFs. We see that the runtime of both full INFs and INFs (CA) stay relatively constant, but for the other baselines, it grows polynomially. This demonstrates the benefits of our INFs for scalability.
Table 8: Total runtime (in minutes) of the experiments using ACIC 2018 datasets with different sample sizes. Reported: mean + standard deviation four datasets and five runs for each size (lower better). Experiments are carried out on Intel(R) Xeon(R) Silver 4316 CPU @ 2.30GHz..
Sample size
1000
2500
5000
10000
25000
50000
TARNet"
0.30 0.01
0.33 + 0.01
0.43 = 0.01
0.64 + 0.03
2.24 + 0.12
8.16 + 0.22
MDNs
0.29 0.05
0.34
=
1.02
5.06
0.71
16.24 1.83
CNF [ INFs w/o stud flow]
0.58 0.07
0.82 0.11
1.14 = 0.17
2.33:...2
8.95 + 1.71
29.38 + 5.55
KDE Kim 2018
+
DKME Muandel el al. 2021
0.03 0.01
0.07 1 0.04
0.13 1 0.05
0.29 1: 0.09
2.28 1 0.35
11,51 1: 0.80
INPs (main)
2.47 0.07
2.47 1 0.09
2.48 +: 0.08
2.52 4: 0.07
2.59 1: 0.09
2,70 + 0.12
31

---

Under review as a conference paper at ICLR 2023
L CASE STUDY: CALIFORNIA'S TOBACCO CONTROL PROGRAM
Overview, To show a real-world application of our INFs, we provide additional results using a case study where we evaluate the effect of California's 's tobacco control program (Abadie et al.) 2010). This refers to the effect of Proposition 99, a large-scale tobacco control program introduced in California after 1988. Proposition 99 increased California's ' cigarette tax by 25 cents per pack, and earmarked the tax revenues to health and anti-smoking education. The main conclusion of[Abe et al. (2010) is
that the effects of the tobacco control program are much larger than previously reported. The dataset has also found widespread use in causal inference ever since (e.g.. Bellot & van der Schaar 2021 In the original paper (Abadie et al. 2010), the results were based on a synthetic control method but without providing density estimates.
Dataset. After an initial pre-processing, the dataset consists of the 39 states, including California. For each state, observe several covariates (e. g., beer consumption per capita, GDP per capita, retail price, and percent of people aged 15-24) and the outcome, i. e., cigarette sales per capita. These recorded annually for each year from 1970 to 2001. Further details on the datasets are in (Abadie ett al. 2010
To apply our INFs, we make several gross assumptions. First, as there is only one treated state, it impossible to satisfy the positivity assumption. Therefore, we consider a tuple (state, year) as an independent unit of measurement, thus obtaining n 1209 observations with 12 treated observations (i. c., those of the state of California after 1989). We also add year as a covariate, which gives dx = 4++1 1. We acknowledge that, even after the previous pre-processing, we still cannot formally guarantee the independence between units of measurement, as the observations of one state over time are not independent. Second, we assume the consistency holds, and there is no spillover effect between neighboring states, so that the potential outcome of one state is independent of the others.
Results. We plot the empirical conditional and the estimated interventional distributions in Fig. g The results go in line with the conclusion in (Abadie et al., 2010). Our main finding is that the introduction of the Proposition 99 (a 1) to all the states from 1970 would substantially reduce tobacco sales. In particular, the mass of the interventional density is shifted to the left which accounts for the reduction of the consumption.
As robustness check, we analyze the role of the smoothness hyperparameter. Our conclusion remains consistent if one specifics different smoothness hyperparameter for the student flow, i. e. 7kkooss 5 and 10. The specification of this hyperparameter is based on the prior knowledge of researcher and cannot be chosen via observational data. However, we find consistent evidence of a positive effect.
5
Rinete,
0.030
0.080
P-Y v =
2(Y 4 0)
P.Y v 1-
0.025
2(Y -pl--
pa-1oow (Y100 v)
34-10rw (1'10
$4-1PIW (Y'
34-1P1W (1'11) M)
0.015
0.015
0.010
0.010
0.005
0.005
0.000
0.000
30
100
1.50
250
300
W. per-capita cigurette sales from 1970 2000 (in packs)
per-capita cigurette sales from 1970 to 2000 (in packs)
Figure& 8: Empirical ground-truth conditional and estimated interventional distributions of cigarette sales per capita from 1970 to 2001. Treatment a = corresponds to the introduction the Proposition 99, that is, comprehensive tobacco tax along with educational programs. We plot our INF' density estimator, IPA (Y(a)) with different smoothness hyperparameter values Agnots of the student flow.
32


Reviews:
[{"review_id": "pkHvXvaWAE", "reviewer": "Reviewer_A7uk", "summary": "", "questions": "", "limitations": "", "rating": 5, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "The authors propose a fully-parametric, deep-learning method for interventional density estimation, called Interventional Normalizing Flows (INFs). INFs provide a properly normalized density estimator. The authors further develop a two-step training procedure with a one-step bias correction for efficient and doubly robust estimation.", "summary_of_the_review": "The method proposed in this work is novel. The authors estimate the density of potential outcomes after interventions from observational data and propose a fully-parametric deep learning method INF. The authors also provide extensive numerical studies to support the performance of INF. ", "strength_and_weaknesses": "To improve the readability, the authors may want to: \n1. Add more introduction to the structural causal model.\n2. In the example provided in Introduction Section,  explain interventional distributions, observational distributions, and counterfactual distributions. Explicitly define the notation Y [a].\n\n\nQuestion: Why does Table 3 compare % of runs with best performances, while other tables compare the log-probability? And why are the columns for Table 3 named log-prob?\n", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "4: All of the claims and statements are well-supported and correct.", "technical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work.", "empirical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work."}, "extra_text": {}}, {"review_id": "W6bLZk-6_-", "reviewer": "Reviewer_uUEZ", "summary": "", "questions": "", "limitations": "", "rating": 8, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "The paper introduces an optimization framework with a model consisting of a teacher and student normalizing flows for estimating the interventional distribution for the case of backdoor adjustment. The developed optimization uses one-step bias correction and utilizes influence functions for doubly robust estimation of the interventional density. Results on multiple synthetic and semi-synthetic datasets are reported showing that INFs outperform nonparametric and ablated methods.", "summary_of_the_review": "The paper has good quality, including further information (summarized under comments) can make it stronger.", "strength_and_weaknesses": "**Strengths**\n* The presented methods fill a gap in the literature on estimating interventional densities using parametric models. Therefore the problem considered is timely and significant.\n* The motivation is clear as the example in Fig. 1 is representative of why the problem considered is important.\n* I found the schematic in Fig. 2 very helpful, it provides a nice summary of the components of the model.\n* The presentation and the logical flow of the paper are clear, I specifically appreciated the use of colors in the text for determining what corresponds to nuisance parameters and what does not.\n* The background contains a nice summary of the existing literature for estimating interventional effects for the case of backdoor adjustment. The literature has focused on average treatment effect (ATE) and the non-ATE literature is mostly non-parametric or non-algorithmic.\n\n**Comments**\n* The problem of estimating $\\mathbb{E}_{X \\sim P(X)} \\big[ P(Y|X, A) \\big]$ given samples from $X, Y, A$ sounds like a well studied estimation problem in statistics regardless of the context of interventional distributions (here, the interventional or counterfactual data is not used and we are interested in the above estimation problem). To me, it is surprising that Normalizing Flows are not used for this purpose before so I cannot confirm the novelty of the methodology. Given that (conditional) normalizing flows, the theory of semi-parametric estimation of the interventional distribution, doubly robust estimation tools existed before the main contribution seems to lie in the addition of the student network, and framing and solving the optimization problem.\n\n* Given that the work is applicable only to the backdoor adjustment, even for the cases in which we can identify the effect of intervention we cannot use this method. However, the extension to the more general case sounds feasible. This would largely broaden the impact and scope of the paper. Can the authors discuss the possibility of this extension?\n\n* What happens if $Y$ is not one-dimensional? It seems that the method relies on creating a grid on $\\mathcal{Y}$. Is this applicable and does this scale to higher dimensional $Y$ variables?\n\n* Can the authors explain why in the case of the synthetic dataset in Fig. 3 KDE is doing very well (even better than the proposed method) for $a=0$ but fails for $a=1$? Where does this asymmetry come from?\n\n* Can the authors include the SCMs used for other examples in the supplementary? \n\n* I think the contributions of this paper by themselves are enough for the acceptance of the paper if these techniques are not used in the context of causal inference and interventional density estimation before. I am happy to change my score to acceptance if the authors include the above information.\n\n* Can the authors include the famous tobacco control example in the experiments?\n\n* To make this a standalone paper, I encourage the authors to include a short summary of the influence functions and doubly robust estimation. More specifically, to me, it wasn't quite clear how we go from Eq. 5 to 7, 8.\n\n\n**Post rebuttal**\n\nI thank the authors for clarifying their contributions and distinguishing between the original estimation theory and the semi-parametric giving rise to the optimization framework and the suggested student-teacher architecture. Now it's clear why the classical density estimation techniques aren't as powerful as semi-parametric methods and why finite-dimensional estimands are not applicable. I think the new additions (adding more background on semi-parametric theory and SCMs, including the SCMs associated with the experiments, adding results on the tobacco control experiment) make the paper more standalone and powerful. Hence I changed my score to acceptance.", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "4: All of the claims and statements are well-supported and correct.", "technical_novelty_and_significance": "2: The contributions are only marginally significant or novel.", "empirical_novelty_and_significance": "4: The contributions are significant, and do not exist in prior works."}, "extra_text": {}}, {"review_id": "DzdO-TimGf", "reviewer": "Reviewer_quXx", "summary": "", "questions": "", "limitations": "", "rating": 6, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "The paper considers an interesting problem in interventional causal inference, which is estimating the distribution of a potential outcome corresponding to an intervention. It proposes a fully parametric deep learning model for this task called Interventional Normalizing Flow that uses 2 flows. \n\nThe teacher flow computes the propensity score $\\pi\\_{a}(X)$ from $X$ and parameters of a conditional normalizing flow network (CNF) from $(X, A)$. This conditional flow network maps Gaussian noise to $P(Y|X, A)$.\n\nThe student flow is a set of unconditional normalizing flows for each intervention $A=a$. Each unconditional flow maps Gaussian noise to $g(Y, \\beta\\_a)$ where $\\beta\\_a$ is the set of parameters corresponds to $A=a$. Each unconditional flow w.r.t. $A=a$ is learned by mimicking $P(Y|do(A=a))$ obtained from the teacher flow. This problem is transformed into an equivalent problem of solving moment conditions. This paper also introduces an (optional) one-step bias correction for doubly robust estimation of $\\beta\\_a$ by leveraging the propensity score prediction from the teacher flow.\n", "summary_of_the_review": "Overall, this is a good paper which proposes a novel method to address an interesting problem. The only problem is that the model is more complex than it should be. This can limit the practical applications of the work. I give this a paper a weak accept.", "strength_and_weaknesses": "### Strengths\n1) The problem considered in the paper is interesting and the method is novel.\n2) The paper is well-written with a good discussion of related works especially those related to interventional density estimation and normalizing flows for causal inference (Appendix A). I also appreciate the authors' attempt to make things clear for the readers by using colors in their formulas.\n3) I think the experiment is quite extensive with various baselines for density estimation. The authors also compare different variants of their method.\n\n### Weaknesses\n1) The main weakness is the complexity of this method as in my opinion, it can be much simpler. For example, the student flow seems redundant as it just approximates $P(Y|do(A=a))$ though I know the authors want to avoid taking an average over $X$ during testing. In the ablation studies, the authors also consider a variant of the proposed method without the student flow, and when looking at the results in Table 2, I saw that this variant is not too bad compared to the full model.\n2) Please correct me if I am wrong but I think we can use simple methods like Balancing Linear Regression [1], TARNet [2], DragonNet [3], … for interventional density estimation (IDE) problem in the paper. These methods are mainly designed for ITE estimation but we can easily modify their outputs to model $P(Y|X, A)$. Once we can model P(Y|X, A), we can simply estimate $P(Y|do(A)) = \\mathbb{E}_{p(X)}[P(Y|X, a)]$ via backdoor adjustment. I would like to hear the authors’ opinions about this.\n3) The authors should compare their method with existing ITE estimation methods to support their argument in the paper saying that “estimating the potential outcome density is better than just estimating the average potential outcome“.\n4) In the TeacherFlow, I would like to understand why did the authors use $X, A$ to compute the weight $\\theta$ of the conditional normalizing flow? Why don’t just model $P(Y|X, A)$ directly via a neural network with $X, A$ are inputs to this network? Since the distribution of outcomes $Y$ is usually simple. Do we really need a normalizing flow with hyper-network for its parameters?\n5) The authors should provide more insights and explanations about why each component is useful for the method to work. Currently, they only use one metric, which is the log-likelihood, for comparison, which I think is not enough.\n\n[1] Learning Representations for Counterfactual Inference, Johansson et al., ICML-2016\n\n[2] Estimating individual treatment effect: generalization bounds and algorithms, Shalit et al., ICML-2017\n\n[3] Adapting Neural Networks for the Estimation of Treatment Effects, Shi et al., NIPS-2019\n", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "3: Some of the paper’s claims have minor issues. A few statements are not well-supported, or require small changes to be made correct.", "technical_novelty_and_significance": "4: The contributions are significant, and do not exist in prior works.", "empirical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work."}, "extra_text": {}}]

Author rebuttal:

Decision: Reject