# Advancing Counterfactual Inference through Nonlinear Quantile Regression

Shaoan Xie  Biwei Huang  Bin Gu  Tongliang Liu  Kun Zhang

###### Abstract

The capacity to address counterfactual "what if" inquiries is crucial for
understanding and making use of causal influences. Traditional counterfactual
inference, under Pearls’ counterfactual framework, typically depends on having
access to or estimating a structural causal model. Yet, in practice, this
causal model is often unknown and might be challenging to identify. Hence,
this paper aims to perform reliable counterfactual inference based solely on
observational data and the (learned) qualitative causal structure, without
necessitating a predefined causal model or even direct estimations of
conditional distributions. To this end, we establish a novel connection
between counterfactual inference and quantile regression and show that
counterfactual inference can be reframed as an extended quantile regression
problem. Building on this insight, we propose a practical framework for
efficient and effective counterfactual inference implemented with neural
networks under a bi-level optimization scheme. The proposed approach enhances
the capacity to generalize estimated counterfactual outcomes to unseen data,
thereby providing an upper bound on the generalization error. Furthermore,
empirical evidence demonstrates its superior statistical efficiency in
comparison to existing methods. Empirical results conducted on multiple
datasets offer compelling support for our theoretical assertions.

Machine Learning, ICML

##  1 Introduction

Understanding and making use of cause-and-effect relationships play a central
role in scientific research, policy analysis, and everyday decision-making.
Pearl’s causal ladder (Pearl, 2000) delineates the hierarchy of prediction,
intervention, and counterfactuals, reflecting their increasing complexity and
difficulty. Counterfactual inference, the most challenging level, allows us to
explore what would have happened if certain actions or conditions had been set
to a different value, providing valuable insights into the underlying causal
relationships between variables.

Figure 1: Illustration of our proposed quantile-based counterfactual
estimation ($Z$ is omitted for illustration purpose). For a sample of interest
($X=x,Z=z,Y=y$), we estimate the quantile $\tau=P(Y\leq y|X=x,Z=z)=0.70$ with
factual observations. Then the counterfactual outcome $Y_{X=x^{\prime}}$ is
equal to the value $y^{\prime}$ which satisfy $P(Y\leq
y^{\prime}|X=x^{\prime},Z=z)=\tau$.

Conventional approaches to counterfactual inference often rely on having
access to or estimating a structural causal model (SCM). For instance, within
Pearls’ framework (Pearl, 2000), if the SCM is given, the typical three-step
procedure involves first estimating the noise value of the sample of interest,
then modifying the model, and finally computing the counterfactual value using
the estimated noise within the adjusted SCM. Unfortunately, the structure
causal model is often unavailable. To tackle this issue, various approaches
have been proposed for estimating the SCM using observational data (Pawlowski
et al., 2020; De Brouwer, 2022; Sanchez & Tsaftaris, 2022; Song et al., 2020;
Chao et al., 2023; Lu et al., 2020a; Nasr-Esfahany et al., 2023; Khemakhem et
al., 2021; Xia et al., 2022; De Lara et al., 2021; Melnychuk et al., 2022).
However, general causal models without specific functional constraints may
lack identifiability (Zhang et al., 2015). In other words, there is no
guarantee that the learned model is the true SCM, making the counterfactual
predictions unreliable. Moreover, current methods often rely on specific
distributional assumptions regarding the noise terms, thereby restricting the
model class.

This paper aims to tackle the challenge of reliable counterfactual inference
without relying on a predefined structural causal model or even direct
estimations of conditional distributions. To this end, we establish a novel
connection between counterfactual inference and quantile regression.
Specifically, suppose there are three random variables $X,Y,Z$ and a latent
noise variable $E_{Y}$, where $X,Z,E_{Y}$ cause the outcome $Y$, and suppose
that the observed evidence is $\langle X=x,Y=y,Z=z\rangle$. We find that the
counterfactual outcome of $Y$ if $X$ had been intervened to $x^{\prime}$ given
the observed evidence, denoted by $Y_{X=x^{\prime}}|X=x,Y=y,Z=z$, is equal to
the $\tau$-th quantile of the conditional distribution
$P(Y|X=x^{\prime},Z=z)$, where $Y=y$ is the $\tau$-th quantile of
$P(Y|X=x,Z=z)$, i.e., $\tau=P(Y\leq y|X=x,Z=z)$, under mild conditions.
Consequently, we reframe counterfactual inference as an extended quantile
regression problem, which yields improved statistical efficiency compared to
existing methods. An illustration is provided in Fig. 1. Furthermore, our
approach enhances the generalization ability of estimated counterfactual
outcomes for unseen data and provides an upper bound on the generalization
error. Our contributions can be summarized as follows.

  * •

We introduce a novel framework that formulates counterfactual inference as an
extended quantile regression problem. This framework is implemented using
neural networks within a bi-level optimization scheme, offering superior
statistical efficiency.

  * •

We assess the generalization capacity of our proposed approach to previously
unseen data and establish an upper bound on the generalization error.

  * •

We conduct extensive experiments to validate our theories and showcase the
effectiveness of our proposed method in diverse scenarios.

##  2 Problem Formulation and Related Work

In this section, we introduce key concepts relevant to our study, including
Pearl’s three-step procedure for counterfactual inference, the technique of
quantile regression, and recent works in counterfactual inference under
Pearl’s procedure. Below, we first give a formal definition of counterfactual
outcomes.

######  Definition 1 (Counterfactual outcomes (Pearl, 2000)).

Suppose $X$, $Y$, and $Z$ are random variables, where $X$ causes $Y$, and $Z$
is a set of common causes to $X$ and $Y$. Given observations $\langle
X=x,Y=y,Z=z\rangle$, the counterfactual outcome of $Y$ is defined as the value
of $Y$ if $X$ had been set to a different value $x^{\prime}$ and is
mathematically represented as $Y_{X=x^{\prime}}|Y=y,X=x,Z=z$.

#### Pearl’s Three-Step Procedure for Counterfactual Inference

In the context of a SCM, Pearl (2000); Pearl et al. (2016) introduced a three-
step procedure to address such counterfactual reasoning.

Suppose the SCMs $Y=f_{Y}(X,Z,E_{Y})$, $X=f_{X}(Z,E_{X})$, and $Z=E_{Z}$ are
given, denoted by $M$, and that we have evidence $\langle X=x,Y=y,Z=z\rangle$.
The following steps outline the process of counterfactually inferring $Y$ if
we had set $X=x^{\prime}$ (Pearl, 2000; Pearl et al., 2016):

  * •

Step 1 (abduction): Use the evidence $\langle X=x,Y=y,Z=z\rangle$ to determine
the value of the noise $E_{Y}=e$.

  * •

Step 2 (action): Modify the model, $M$, by removing the structural equations
for the variables in $X$ and replacing them with the functions $X=x^{\prime}$,
thereby obtaining the modified model, $M_{x^{\prime}}$.

  * •

Step 3 (prediction): Use the modified model, $M_{x^{\prime}}$, and the
estimated noise $e$ to compute the counterfactual of $Y$ as
$Y_{X=x^{\prime}}=f_{Y}(x^{\prime},z,e)$.

Note that Step 1 performs deterministic counterfactual reasoning, focusing on
counterfactuals of a single unit of the population, where the value of $E_{Y}$
is determined.

#### Deep Counterfactual Inference

Various deep-learning approaches have been proposed for estimating the SCM and
noise values using observational data and accordingly perform the three-step
procedure for counterfactual inference. (Khemakhem et al., 2021) proposes to
use autoregressive flow to perform causal discovery by comparing the
likelihood and infer the noise for counterfactual inference by inverting the
flow. (Javaloy et al., 2023) presents a mechanism to embed additional causal
knowledge when learning the causal autoregressive flow. CFQP (De Brouwer,
2022) considers a setting when the background variables are categorical and
employs the Expectation-Maximization framework to predict the cluster of the
sample and perform counterfactual inference with the regression model trained
on the specific cluster. CTRL (Lu et al., 2020a) and BGM (Nasr-Esfahany et
al., 2023; Nasr-Esfahany & Kiciman, 2023) show that the counterfactual outcome
is identifiable when the SCM is monotonic w.r.t. the noise term. In
particular, BGM uses conditional spline flow to mimic the generation process
and performs counterfactual inference by reversing the flow. DeepSCM
(Pawlowski et al., 2020) proposes to use variational inference and normalizing
flow to infer the noise variable. DiffSCM (Sanchez & Tsaftaris, 2022) proposes
to match the observation distribution with a conditional diffusion model and
infer the noise by reversing the diffusion process, but it only allows one
cause. DCM(Chao et al., 2023) generalizes DiffSCM to support multiple causes
of a single variable. G-NCM (Xia et al., 2022) extends the neural causal model
(Xia et al., 2021) to estimate the counterfactual distribution. (Melnychuk et
al., 2022) proposes a counterfactual domain confusion loss to address
confounding bias and uses a transformer to perform counterfactual inference
for long-range time-series data. (Ribeiro et al., 2023) proposes a
hierarchical latent mediator model for counterfactual image generation task.
Individual treatment effect (ITE) (Johansson et al., 2016; Yoon et al., 2018;
Bica et al., 2020; Yao et al., 2018; Li & Yao, 2022; Lu et al., 2020b; Zhou et
al., 2021b, a) is also deeply connected to counterfactual inference while the
former focuses on the differences between expected outcomes over the
population before and after intervention. We present more ITE works in the
supplement.

#### Quantile Regression

Traditional regression estimation focuses on estimating the conditional mean
of $Y$ given $X$, typically represented by the function $f(X)$. On the other
hand, quantile regression (Koenker & Hallock, 2001) is concerned with
estimating conditional quantiles, specifically the $\tau$-th quantile
$\mu_{\tau}$, which is the minimum value $\mu$ such that $P(Y\leq\mu|X)=\tau$,
where $\tau$ is a predefined value. Some quantile regression settings
encounter the quantile crossing problem (Takeuchi et al., 2009), and
Tagasovska & Lopez-Paz (2019) proposes a loss function to learn all the
conditional quantiles of a given target variable with neural networks to
address this issue. Fortunately, as our framework only requires learning a
single quantile for the observation which contains important information about
the noise term, this issue can be avoided.

##  3 Quantile-Regression-based Counterfactual Inference: Theoretical Insights

Conventional counterfactual inference approaches typically rely on estimating
both structural causal models and noise values. However, the simultaneous
estimation of these elements can be challenging. In this paper, we establish a
novel connection between counterfactual inference and quantile regression,
inspired from the reinforcement learning literature Lu et al. (2020a). This
connection allows us to bypass the need for estimating both structural causal
models and noise values. Specifically, the counterfactual outcome
$Y_{X=x^{\prime}}|X=x,Y=y,Z=z$ is equal to the $\tau$-th quantile of the
conditional distribution $P(Y|X=x^{\prime},Z=z)$, with $Y=y$ representing the
$\tau$-th quantile of $P(Y|X=x,Z=z)$. Exploiting this connection, the
counterfactual outcome can be directly estimated through quantile regression
alone. This fundamental relationship is formalized in the theorem below.

######  Theorem 1.

Suppose a random variable $Y$ satisfies the following structural causal model:

$\vspace{-1mm}Y=f(X,Z,E)\vspace{-.5mm}$

where $X$ and $Z$ cause $Y$, with $Z$ being a cause to $X$. $E$ is the noise
term, indicating some unmeasured factors that influence $Y$, with
$E\mathchoice{\mathrel{\hbox to0.0pt{$\displaystyle\perp$\hss}\mkern
2.0mu{\displaystyle\perp}}}{\mathrel{\hbox
to0.0pt{$\textstyle\perp$\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptstyle\perp$\hss}\mkern
2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptscriptstyle\perp$\hss}\mkern
2.0mu{\scriptscriptstyle\perp}}}X,Z$. We assume that $f$ (which is unknown) is
smooth and strictly monotonic in $g(E)$ for fixed values of $X,Z$. Suppose we
have observed $\langle X=x,Y=y,Z=z\rangle$. Then for the counterfactual
inquiry, what would be the outcome of $Y$ if $X$ had been set to $x^{\prime}$,
given the observed evidence $(X=x,Y=y,Z=z)$, the counterfactual outcome
$Y_{X=x^{\prime}}|X=x,Y=y,Z=z$ is equal to the $\tau$-th quantile of the
conditional distribution $P(Y|X=x^{\prime},Z=z)$, where $Y=y$ represents the
$\tau$-th quantile of $P(Y|X=x,Z=z)$.

The theorem establishes a connection between counterfactual outcomes and
quantiles, providing the groundwork for identifying counterfactual outcomes
through quantile regression from purely factual observations under mild
assumptions. In particular, it relies on the assumption that the function $f$
is strictly monotonic w.r.t. $g(E)$, where $g(\cdot)$ is any arbitrary
function of the latent factor $E$. Note that, equivalently, $Y=f(X,Z,E)$ can
be represented as $f_{1}(f_{2}(X,Z),g(E))$. Furthermore, this theorem
additionally brings the following two important theoretical and practical
advantages compared to existing approaches for counterfactual inference.

Counterfactual outcome is identifable even when $f$ is not identifiable.
Numerous methods have attempted to identify the true generating function $f$
and infer the noise value by inverting the estimated function. However, the
structural causal model may not always be identifiable (Zhang et al., 2015).
In contrast, even if the structural causal model lacks identifiability,
quantile regression can still be employed to identify the counterfactual
outcome under our monotonicity assumption. Additionally, note that the theorem
remains valid regardless of whether the variables involved are continuous or
discrete.

No need to recover the true noise $E$. An empirical approach for
counterfactual inference based on the above theorem will be introduced in
Section 4, which avoids the estimation of the value of the noise term.
Moreover, introducing monotonicity on $g(E)$, rather than directly on $E$, not
only relaxes the constraints on structural causal models but also on the
distribution of the unobserved term. In contrast, recent counterfactual
methods such as BGM (Nasr-Esfahany & Kiciman, 2023) and DCM (Chao et al.,
2023), although considering an invertible function of the noise term, often
assume a specific distribution for the noise term (e.g., Gaussian) and rely on
the estimation of such noise values. This assumption makes them inflexible
when encountering different forms of $g$.

For instance, consider the following two scenarios: in the first scenario, let
$Y=\frac{1}{4}(Z+\sin(2\pi X)g(E)+2g(E))$, $E$ satisfies a Gaussian
distribution, and $g(E)=E$, while in the second scenario,
$g(E)=\min(\max(E,-0.5),1)$. For a particular sample $\langle
X=0.5,Z=0.5\rangle$ (the unknown noise $E=0.5$), we have $Y=1.5$ for both
scenarios. Now we estimate the counterfactual outcome when $X$ had been set to
different values, with results illustrated in Figure 2. We can see our
approach, which relies on Theorem 1, remains unaffected by the form of $g(E)$,
since in both scenarios, it relies on $P(Y\leq 1.5|X=0.5,Z=0.5)$. However, BGM
and DCM exhibit vastly different performances in the two scenarios because
they assume Gaussian noise terms, while it is impossible to recover a noise
term that is both Gaussian and satisfies invertibility when
$g(E)=\min(\max(E,-0.5),1)$.

(a) $g(E)=E$.

(b) $g(E)=min(max(E,-0.5),1)$

Figure 2: Comparisons of the counterfactual predictions when traversing $X$
under two different forms of $g(E)$. Our method, based on Theorem 1,
demonstrates resilience to various forms of $g(E)$.

Theorem 1 relies on the monotonicity condition on $g(E)$ and covers a wide
range of cases. Below, we present a compilation of commonly encountered
special cases where this condition remains valid.

  * •

Linear causal models: $Y=aX+bZ+g(E)$.

  * •

Nonlinear causal models with additive noise: $Y=f(X,Z)+g(E)$.

  * •

Nonlinear causal models with multiplicative noise: $Y=f(X,Z)\cdot g(E)$.

  * •

Post-nonlinear causal models: $Y=h(f(X,Z)+g(E))$.

  * •

Heteroscedastic noise models: $Y=f(X,Z)+h(X,Z)\cdot g(E)$.

##  4 Quantile-Regression-based Counterfactual Inference: Practical Approaches

Building on the theoretical results presented in Section 3, we have
established that counterfactual inference can be reframed as an extended
quantile regression problem. This reformulation eliminates the requirement to
estimate the structural causal model and noise values for addressing
counterfactual inquiries. Accordingly, in this section, we introduce a
practical approach for counterfactual inference with quantile regression and
establish an upper bound on the generalization error.

In particular, in Section 4.1, we formulate counterfactual inference under a
bi-level optimization scheme. The upper level is dedicated to estimating the
targeted quantile level $\tau$, while the lower level endeavors to estimate
the quantile regression function at the specific quantile level. To solve this
bi-level optimization problem, Section 4.2 introduces a practical estimation
approach that employs neural networks capable of accommodating general causal
models and data distributions. Additionally, for more efficient inference of
multiple counterfactual samples, we adopt a compact representation to
encapsulate the variation in quantile regression functions for different
samples, eliminating the need to retrain the bilevel optimization for each new
sample. Furthermore, Section 4.3 provides an in-depth analysis of the
generalization ability of the proposed approach and establishes an upper bound
on the generalization error.

###  4.1 Counterfactual Inference under a Bi-Level Optimization Scheme

Suppose there are $N$ samples $\\{x_{i},y_{i},z_{i}\\}_{i=1}^{N}$ which are
realizations of random variables $X$, $Y$, and $Z$. We are interested in
finding the counterfactual outcome $y^{\prime}$ which is the realization of
$Y_{X=x^{\prime}}$ for a particular sample of interest $\langle x,y,z\rangle$.
An illustration of our method is provided in Fig. 1: we first estimate the
quantile $\tau=P(Y\leq y|x,z)$ and its corresponding quantile function
$\mu_{\tau}(X,z)=\min_{\mu}[P(Y\leq\mu|X,z)=\tau]$. Then we can infer the
counterfactual outcome with $\mu_{\tau}(x^{\prime},z)$.

However, estimating $\tau$ and $\mu_{\tau}$ remains challenging. A
straightforward way is to estimate $P(Y\leq y|X=x,Z=z)$ as $\tau$ first, e.g.,
with Monte Carlo, and then perform standard quantile regression to obtain the
corresponding quantile function. However, $P(Y\leq y|x,z)$ can be difficult to
estimate with finite training samples. For example, there may be few or even
only one training sample (itself) that have $X=x,Z=z$, leading to inaccurate
estimation of $\tau$, which has been demonstrated in the experiments in
Section 5.2.

To address this problem, we couple the estimations of $\tau$ and the quantile
function $\mu_{\tau}$, and formulate the counterfactual inference problem as a
bi-level optimization problem:

$\displaystyle\hat{\tau}$
$\displaystyle=\arg\min_{\tau}{L}(\hat{\mu}_{\tau}(x,z),y),\textrm{with}$ (1)
$\displaystyle\hat{\mu}_{\hat{\tau}}$
$\displaystyle=\arg\min_{\mu}R_{\hat{\tau}}[\mu],$

where the upper level is to estimate the quantile $\tau$ and $L(.,.)$ can be a
regular regression loss, such as $L_{1}$ and $L_{2}$ loss. As for the lower
level problem, it is a standard quantile regression problem for the particular
quantile level, with the quantile regression function:

$\displaystyle
R_{\tau}[\mu]=\frac{1}{N}\sum\nolimits_{i=1}^{N}l_{\tau}(y_{i}-\mu(x_{i},z_{i}))$
(2) $\displaystyle\text{and }l_{\tau}(\xi)=\begin{cases}\tau\xi,&\text{ if
}\xi\geq 0\\\ (\tau-1)\xi,&\text{ if }\xi<0.\end{cases}\vspace{-2mm}$ (3)

Figure 3: One training step of our proposed bi-level implementation.

We use the pin-ball loss $l_{\tau}$ as the objective as it has been shown that
the minimizer of the empirical pin-ball loss converges to the true quantile
$\mu_{\tau}$ under some mild assumptions (Takeuchi et al., 2006).

After the bi-level optimization, we obtain an accurate quantile regressor
$\hat{\mu}_{\tau}$. Since the factual outcome $y$ is the $\tau$-th quantile of
the $P(Y|x,z)$, $\tau$ is the minimizer of the objective
$L(\hat{\mu}_{\tau}^{*}(x,z),y)$. In other words, we have
$\hat{\tau}\\!\\!=\\!\\!\tau\\!\\!=\\!\\!P(Y\leq
y|x,z),\hat{\mu}_{\hat{\tau}}\\!=\\!\mu_{\tau}$. Through the lens of the above
bi-level optimization formulation, we avoid the direct estimation of $P(Y\leq
y|x,z)$ with finite samples. Next, we show how to efficiently solve this bi-
level optimization problem in practice.

###  4.2 An Efficient Neural-Network-Based Implementation

Although we have identifiability guarantees for the counterfactual outcome, it
still remains unclear how should we implement the framework to solve the bi-
level optimization problem. In this context, we present a scalable and
efficient neural network-based implementation, for both lower-level and upper-
level optimization. We employ a neural network $h$ to learn the quantile of
the sample and network $g$ to learn the quantile function. A training pipeline
is given in Fig. 3.

#### Lower-level optimization.

Since each sample of interest $x$ corresponds to a different quantile $\tau$
and quantile function $\mu_{\tau}$, the computational cost can be huge if we
learn the quantile regression function $\hat{\mu}_{\hat{\tau}}$ for every
interested sample separately. Hence, to achieve efficient quantile estimation
of counterfactual inference for multiple samples, we propose to learn a
conditional neural network $g$. Specifically, we use $\hat{\tau}$ to capture
the difference in different samples and concatenate $\hat{\tau}$ with each
training sample $\\{x_{i},z_{i}\\}_{i=1}^{N}$ as input to the network $g$ as
shown in the right part of Fig. 3. It eliminates the need to retrain the
bilevel optimization for each sample of interest:

$\vspace{-2mm}\hat{\mu}_{\hat{\tau}}(x_{i},z_{i})\Longrightarrow
g(x_{i},z_{i},\hat{\tau}).$

Then for each $\hat{\tau}$, we minimize the pin-ball loss
$\frac{1}{N}\sum_{i=1}^{N}l_{\hat{\tau}}(y_{i}-g(x_{i},z_{i},\hat{\tau}))$
following Eq. 4.1. Accordingly, after the optimization procedure, we have the
$\hat{\tau}$-th quantile regression output as $g(x,z,\hat{\tau})$ for every
sample of interest $\langle x,z,y\rangle$ with a shared neural network $g$.

#### Upper-level optimization.

To recover the quantile $P(Y\leq y|X=x,Z=z)$ from the factual observations
$\langle x,y,z\rangle$, we propose a data-dependent model to learn the
quantiles for each sample of interest automatically. Specifically, we employ a
neural network $h$ to infer $\tau$ from the observational data, i.e.,
$\hat{\tau}=h(x,z,y)$, and use $\hat{\tau}$ as the input of the network $f$ in
the lower-level problem to perform quantile regression. An important advantage
of the data-dependent implementation is that it allows inferring $\tau$ and
counterfactual outcomes for an unseen data sample $\langle x,z,y\rangle$. For
example, we may infer the $\tau$ for samples in the validation split even
though they have not been used in the training.

The optimization of the upper-level problem is more challenging as we also
need to consider the lower-level constraint besides the upper-level regression
loss $L(.,.)$. Thanks to the multi-level optimization library–Betty (Choe et
al., 2023b), we are able to solve this bi-level optimization problem without
considering the complex interactions between the lower and upper-level
problems. It uses hyper-gradient (Choe et al., 2023a) to update the network
$h$ automatically.

Summary of Training Pipeline. We present the training pipeline of one step in
Fig. 3. In each training step, we have four sub-steps: 1) fix the estimated
quantiles (i.e., network $h$); 2) train the neural network $g$ with the pin-
ball loss for each $\hat{\tau}$ by stochastic gradient descent (SGD) $m$ times
($m=30$ in our experiments); 3) send the optimized network $g$ to the upper-
level problem. 4) the bi-level optimization library updates the network $h$
one time by hyper-gradient to minimize the regression loss $L(.,.)$. Then we
continue the training steps until the models converge.

###  4.3 Generalization Bound of Empirical Estimator

Ideally, the counterfactual outcome $Y_{x^{\prime}}|X=x,Z=z,Y=y$ is
$\mu_{\tau}(x^{\prime},z)$. However, we are usually given a limited number of
training samples and the pair $(\tau,\mu_{\tau})$ is approximated by using the
bi-level optimization solution $(\hat{\tau},\hat{\mu}_{\hat{\tau}})$. For a
sample of interest $\langle x,y,z\rangle$, we have the estimated quantile
$\hat{\tau}=h(x,y,z)$ and counterfactual predictions
$y^{\prime}=\hat{\mu}_{\hat{\tau}}(x^{\prime},z)=g(x^{\prime},z,\hat{\tau})$
for any new value $x^{\prime}$. An essential problem is that we are unsure
about the generalization ability of the regressor $\hat{\mu}_{\hat{\tau}}$ on
the counterfactual input pair $\langle x^{\prime},z\rangle$ as the pair has
not been seen by $\hat{\mu}_{\hat{\tau}}$ during training.

|  |  |  |
---|---|---|---|---
Figure 4: Toy examples for counterfactual estimations. For the interested
sample, we traverse the value of $X$ with $x^{\prime}$ and compare against
with the true counterfactual $Y_{x=x^{\prime}}$. Our method is able to recover
the true trajectory from factual observations under the five causal models.
Figure 5: Box plots of the learned quantiles. The ground truth quantile is
$\Phi(0.5)=0.691$. Compared to the Monte Carlo (MC) method, our proposed
bilevel optimization method learns the quantile stably while the MC method has
a significantly larger variance, which makes it unsuitable for counterfactual
inference.

More formally, we are interested in analyzing the upper bound of the
generalization error
$\mathbb{E}_{x,z}[l_{\hat{\tau}}(\mu_{\hat{\tau}}(x,z)-\hat{\mu}_{\hat{\tau}}(x,z))]$.
Below, we employ the Rademacher complexity to upper bound the generalization
error.

######  Definition 2 (Rademacher complexity (Bartlett & Mendelson, 2002)).

Let $F$ be a hypothesis class mapping from $\mathcal{X}$ to $[0,b]$. Let
$\\{x_{i},z_{i}\\}_{i=1}^{N}$ be i.i.d. examples. Let
$\\{\sigma_{i}\\}_{i=1}^{N}$ be independent Rademacher variables taking values
in $\\{-1,+1\\}$ uniformly. The Rademacher complexity is defined as

$\displaystyle\mathfrak{R}(F)=\mathbb{E}_{x,z,\sigma}\left[\sup_{\mu\in
F}\frac{1}{N}\sum_{i=1}^{N}\sigma_{i}\mu(x_{i},z_{i})\right].$

Our main theoretical result is as follows.

######  Theorem 2.

Let $(\hat{\tau},\hat{\mu})\in(\mathfrak{T},F)$ by the optimization solution,
where $\mathfrak{T}$ is the parameter space. Let the loss function $l_{\tau}$
be upper bounded by $b$. Then, for any $\delta>0$, with probability at least
$1-\delta$, we have

$\displaystyle\mathbb{E}_{x,z}[l_{\hat{\tau}}(\mu_{\hat{\tau}}(x,z)-\hat{\mu}_{\hat{\tau}}(x,z))]$
$\displaystyle\leq\frac{1}{N}\sum_{i=1}^{N}l_{\hat{\tau}}(\mu_{\hat{\tau}}(x_{i},z_{i})-\hat{\mu}_{\hat{\tau}}(x_{i},z_{i}))$
$\displaystyle+4\mathfrak{R}(F)+\frac{4b}{\sqrt{N}}+b\sqrt{\frac{\log(1/\delta)}{2N}}.$

The Rademacher complexity has been widely used to derive generalization error
bounds in the statistical machine learning community (Mohri et al., 2018). Its
upper bound has also been widely studied. If $F$ is an RKHS and the hypotheses
are upper bounded, without any strong assumptions, $\mathfrak{R}(F)\leq
O(\sqrt{1/N})$ (Bartlett & Mendelson, 2002).

The upper bound of
$\mathbb{E}_{x,z}[l_{\hat{\tau}}(\mu_{\hat{\tau}}(x,z)-\hat{\mu}_{\hat{\tau}}(x,z))]$
heavily relies on the (i) empirical value
$\frac{1}{N}\sum_{i=1}^{N}l_{\hat{\tau}}(\mu_{\hat{\tau}}(x_{i},z_{i})-\hat{\mu}_{\hat{\tau}}(x_{i},z_{i}))$,
which can be minimized by utilizing the factual training samples and (ii) the
number of training samples $N$. Given a fixed number of training samples, our
theoretical results imply that the generalization error is bounded, i.e., the
counterfactual predictions will be very close to the ground truth
counterfactual outcome. We also empirically show that our method can achieve
much better performance given same number of samples compared with other
strong baselines in Section 5.4. In summary, the performance of our proposed
quantile regression method is guaranteed with finite samples.

##  5 Experimental Results

In this section, we begin by introducing the experimental setup including the
datasets, evaluation metrics, baseline methods, and implementation details.
Then we analyze the learned quantiles under different models and conduct a
comprehensive comparison with state-of-the-art approaches across diverse
datasets. Additionally, we also study the sample efficiency of our method, the
essential monotonicity assumption, and potential latent confounding cases.

###  5.1 Experiment Setup

Datasets. Due to the lack of observations in counterfactual scenarios in the
real world, we create the following datasets, which satisfy the monotonicity
assumption in our theorem.1, to evaluate the performance of our method. In
Cont-Dose and Dis-Dose, the covariate $Z$ is the age of patients, $X$
represents the (continuous and binary, respectively) dose of the medical
treatment, and $Y$ denotes the outcome. Both datasets contain 800 training and
200 testing samples. To further evaluate our method on high-dimensional data,
we create the Rotation-MNIST and Thick-Omniglot datasets. In Rotation-MNIST,
$Z$ is the original MNIST images (LeCun et al., 2010), $X$ is the rotation of
the digit, and noises are the changes to the RGB values of the rotated images.
In Thick-Omniglot, $Z$ is the original Omniglot images (Lake et al., 2019),
$X$ changes the thickness of the digits, and noises are the darkness of the
transformed images. We also adopt the semi-synthetic IHDP dataset (Hill,
2011).

Implementation. We provide the implementation details in the appendix and code
in the supplement.

Method | IHDP | Cont-Dose | Dis-Dose
---|---|---|---
Train | Test | Train | Test | Train | Test
DeepSCM (Pawlowski et al., 2020) | 2.37 $\pm$ 2. | 2.98 $\pm$ 4. | 0.38 $\pm$ .0 | 0.40 $\pm$ .0 | 0.33 $\pm$ .0 | 0.36 $\pm$ .0
CFQP-T (De Brouwer, 2022) | 1.81 $\pm$.1 | 1.80 $\pm$ .1 | 0.19 $\pm$ .0 | 0.19 $\pm$ .0 | 0.22 $\pm$ .0 | 0.22 $\pm$ .0
CFQP-U (De Brouwer, 2022) | 1.40 $\pm$ .1 | 1.30 $\pm$ .0 | 0.19 $\pm$ .0 | 0.18 $\pm$ .0 | 0.22 $\pm$ .0 | 0.22 $\pm$ .0
BGM (Nasr-Esfahany & Kiciman, 2023) | 4.35 $\pm$ .4 | 4.89 $\pm$ .5 | 0.31 $\pm$ .0 | 0.39 $\pm$ .1 | 0.27 $\pm$.0 | 0.29 $\pm$ .0
DCM (Chao et al., 2023) | 2.56 $\pm$ 2. | 2.76 $\pm$ 2. | 0.19 $\pm$ .0 | 0.16 $\pm$ .0 | 0.28 $\pm$ .0 | 0.29 $\pm$ .0
Ours | 1.29 $\pm$ .3 | 1.23 $\pm$ .2 | 0.06 $\pm$ .0 | 0.06 $\pm$ .0 | 0.20 $\pm$ .0 | 0.20 $\pm$ .0
Table 1: The RMSE performance for counterfactual inference on tabular
datasets. Note that we do not have access to the counterfactual outcomes for
training split during training. So, we measure the counterfactual inference
performances on both training and testing splits.

Method | Rotation-MNIST | Thick-Omniglot
---|---|---
Train | Test | Train | Test
DeepSCM | 6.61 $\pm$ .1 | 6.55 $\pm$ .1 | 11.31 $\pm$ .1 | 11.46 $\pm$ .1
CFQP-T | 2.82 $\pm$ .1 | 2.80 $\pm$ .0 | 4.36 $\pm$ .0 | 4.23 $\pm$ .0
CFQP-U | 1.79 $\pm$ .0 | 1.78 $\pm$ .0 | 3.28 $\pm$ .1 | 3.30 $\pm$ .1
BGM | 7.53 $\pm$ .2 | 7.52 $\pm$ .1 | 12.07 $\pm$ .2 | 12.27 $\pm$ .1
DCM | 5.30 $\pm$ .0 | 5.23 $\pm$ .0 | 8.80 $\pm$ .0 | 9.01 $\pm$ .1
Ours | 1.54 $\pm$ .1 | 1.52 $\pm$ .1 | 2.60 $\pm$ .0 | 2.96 $\pm$ .1

Table 2: The RMSE performance for counterfactual inference on image
transformation datasets.

Z | Y | CFQP | BGM | DCM | Ours | Truth
---|---|---|---|---|---|---
|  |  |  |  |  |
|  |  |  |  |  |
|  |  |  |  |  |
|  |  |  |  |  |
Figure 6: Example results on image transformation dataset.

Evaluation Metrics and Baseline Methods. We compute the root mean square error
(RMSE) between the counterfactual predictions and the ground truth. We compare
with state-of-the-art baselines: DeepSCM (Pawlowski et al., 2020), CFQP (De
Brouwer, 2022), BGM (Nasr-Esfahany & Kiciman, 2023) and DCM (Chao et al.,
2023). As CFQP provides two backbones, we use CFQP-T to denote the method with
a Transformer network and CFQP-U the method with a U-Net backbone. We run the
public code of baseline methods with their recommended hyper-parameters.

###  5.2 Analysis on Quantile Learning

We synthesize the data where the corresponding $\tau$ has closed-form
solutions, so we can compare the learned quantiles against the ground truth.
Specifically, we consider the following instantiations of the five causal
models mentioned in Section 3, including 1) linear causal models, $Y=X+Z+E$;
2) nonlinear additive noise model, $Y=\sin(2\pi X+Z)+E$; 3) nonlinear
multiplicative noise model, $Y=\exp(X-Z+0.5)\cdot E$; 4) post-nonlinear model,
$Y=\exp(\sin(\pi X+Z)+E)$; and 5) heteroscedastic model,
$Y=\exp(-5X+Z)+\exp(X+Z-0.5)\cdot E$. For these five causal models, we have
$F^{-1}(\tau)=E\Rightarrow\tau=F(E)$ where $F$ is the CDF of noise $E$. We set
$P(X)=P(Z)=\mathcal{U}[0,1]$. As for noise, we consider the isotropic Gaussian
distribution $\mathcal{N}(0,1)$. We sample 100,000 data points for each case.
For sample points of interest for counterfactual inference, we use
$X=0.5,Z=0.5,E=0.5$ and generate $Y$. Our goal is to learn a quantile for this
sample, where the ground-truth quantile $\tau$ is $\Phi(0.5)\approx 0.691$.

Figure 7: Analysis on the monotonicity assumption. (a)$Y=\exp(\cos(\pi
X+3Z)+\cos(E))$; (b)$Y=\exp(\cos(\pi X+3Z)+E^{2})$; (c)$Y=\exp(\cos(\pi
X+3Z)+\text{MLP}(E))$; (d)$Y=\exp(\text{MLP}(X,Z)+\text{MLP}(E))$;
(e)$Y=\exp(\text{MLP}(X,Z)+\text{MLP}(E))$. Although $Y$ is not monotonic
w.r.t $E$, the counterfactual outcome is still identifiable since we may have
$g(E)=E^{2},\cos(E),\text{MLP}(E)$ and $Y$ is monotonic w.r.t $g(E)$.
Empirically, we can see that the counterfactual predictions are very close to
the ground truth in the first four cases. As for the last case (e), it
exhibits deviation from the truth since the monotonicity assumption is
violated and the counterfactual outcome may not be identifiable.

Figure 8: The study of sample efficiency.

First, as shown in the upper row of Fig. 4, we compare with the linear
quantile regressor solved by linear programming with specified quantiles at
0.1, 0.5, and 0.9. We find that linear quantile regressors are limited by
linear parameterization and cannot handle nonlinear generation functions such
as the second causal model. Importantly, the baselines do not know the true
quantile as they did not leverage the factual observation of interest. By
contrast, we observe that our method can achieve accurate counterfactual
predictions under all five different causal models, no matter whether they are
linear or nonlinear. Second, one more straightforward way to estimate the
quantiles would be Monte Carlo (MC). Specifically, we use
$\frac{1}{m}\sum_{i=1}^{m}\mathbb{I}(y_{i}\leq y|x_{i}=x,z_{i}=z)$ to
approximate the quantiles $P(Y\leq y|X=x,Z=z)$, where $m$ is the number of
training samples that satisfy $x_{i}=x,z_{i}=z$; It is worth noting that MC
often failed to find samples with $X=0.5,Z=0.5$ exactly, so we relax it by
finding samples with $X$ and $Z$ in range $[0.49,0.51]$. The problem may be
more severe for discrete $X$ and $Z$ since such relaxation is not applicable
on discrete variables. The estimations are visualized in Fig. 5. We observe
that the variance of learned quantiles in different runs is large with the MC
method for all five simple causal models. By contrast, our bi-level
formulation learns the quantile and quantile function together and estimated
quantiles are much more accurate and stable. In particular, our method allows
learning the quantile functions for different samples simultaneously by
conditioning the estimated quantiles in the lower-level problem, which reduces
the training time significantly.

###  5.3 Comparisons with State-of-the-Art Approaches

Table 1 and Table 2 present the results of the counterfactual inference task
on tabular and image datasets, respectively. Our method achieves the lowest
RMSE across all datasets. We observe that the performance of BGM (Nasr-
Esfahany & Kiciman, 2023) on image transformation datasets is less than
satisfactory. BGM employs a conditional spline flow to transform the input to
Gaussian noise. However, its expressive power may be limited, especially on
high-dimensional image data, compared to our unconstrained neural networks. As
a result, it fails to generate accurate counterfactual predictions. In Fig. 6,
we present examples of counterfactual image generations. In the third and
fourth rows, the strong baseline CFQP (De Brouwer, 2022) fails to preserve the
background color of the outcome image, indicating potential challenges in
accurately identifying the noise value within the sample and consequently
generating unreliable counterfactual predictions. In contrast, our method
effectively learns the transformation and preserves the background color in
the prediction. This encouraging result further strengthens our
identifiability theory.

###  5.4 Generalization Bound and Sample Efficiency

We have theoretically shown that the expectation of the generalization error
is bounded by the empirical risk and the number of samples. To verify the
theorem, we perform an ablation study on the number of training samples. We
work on the Cont-Dose and Dis-Dose dataset and run experiments with different
numbers of samples $N\\!\\!=\\!\\!50,100,200,400,800$ and report the RMSE
averaged over different runs in Fig. 8. We find that the baselines BGM (Nasr-
Esfahany & Kiciman, 2023) and CFQP (Pawlowski et al., 2020) are very sensitive
to the number of training samples. In particular, CFQP exhibits a very large
variance when the number of samples is smaller than 200. A possible reason is
that it assigns samples to different clusters and samples in each cluster
become even smaller. In contrast, our method achieves the best performances
and the variance is very small, which supports our generalization bounds and
demonstrates the sample efficiency of our proposed approach.

###  5.5 Monotonicity Assumption and Latent Confounder Cases

An essential assumption of our method is that there exists a function $g(E)$
such that outcome $Y$ is strictly monotonic w.r.t $g(E)$. Since $g$ can be
arbitrary, our assumption covers a very wide spectrum of function classes, as
demonstrated by the superior empirical performance of our method. However, it
may still be violated in some cases. To investigate this problem, we consider
the five causal models as shown in Fig. 7. In the first four cases, although
$Y$ is not monotonic w.r.t $E$, when we set $g(E)=E^{2},\cos(E)$, or more
complex MLP($E$), $Y$ is monotonic w.r.t $g(E)$. Therefore, according to
Theorem 1, the counterfactual outcome should be identifiable, and the results
in the first four cases strongly support this conclusion. As for the last
case, the monotonicity assumption is violated and the counterfactual
predictions exhibit some deviation from the ground truth, which is expected.
The identifiability of the counterfactual outcome in non-monotonic generation
functions remains an important yet challenging task, which we defer to future
work.

Similar to (Nasr-Esfahany & Kiciman, 2023; Pawlowski et al., 2020; Chao et
al., 2023), we also assume no latent confounders. But it would be interesting
to testify about the performances under the confounder case. We considered
three different confounding cases. 1) $X\leftarrow C\rightarrow Z$; 2)
$Z\leftarrow C\rightarrow Y$ and 3) $X\leftarrow C\rightarrow Y$. $C$ is the
unknown confounder.

Scenario | RMSE
---|---
No C | 0.03 $\pm$ .0
$X\leftarrow C\rightarrow Z$ | 0.05 $\pm$ .0
$Z\leftarrow C\rightarrow Y$ | 0.06 $\pm$ .0
$X\leftarrow C\rightarrow Y$ | 0.22 $\pm$ .1
Table 3: Results with the presence of latent confounder $C$.

The detailed data-generating process is given in the supplementary and the
results are shown in Table 3. We observe that our method is quite robust under
the first three scenarios, which is reasonable as such $C$ does not affect the
causal effect from $X$ to $Y$. However, the performance becomes unsatisfactory
under the scenario $X\leftarrow C\rightarrow Y$, suggesting that the existence
of latent confounders behind $X$ and $Y$ potentially hamper the performances
of our method. Extending our approach to cover such latent confounders will be
our future work.

##  6 Conclusion

Counterfactual inference remains a significant challenge due to the lack of
counterfactual outcomes in the real world. In this paper, we advance
counterfactual inference from a novel perspective through quantile regression.
We first build a connection between counterfactual outcomes and quantile
regression and show that the counterfactual outcome corresponds to certain
quantities under rather mild conditions. Accordingly, we propose an effective
and efficient estimation approach based on neural-network implemented bi-level
optimization and show the generalization bound of the empirical estimator. We
verify our method on multiple simulated and semi-real-world datasets. The
superior performance over state-of-the-art baselines highlights the
effectiveness of our method.

## Impact Statement

This paper presents work whose goal is to advance the field of Machine
Learning. There are many potential societal consequences of our work, none
which we feel must be specifically highlighted here.

## References

  * Bartlett & Mendelson (2002) Bartlett, P. L. and Mendelson, S.  Rademacher and gaussian complexities: Risk bounds and structural results.  _Journal of Machine Learning Research_ , 3(Nov):463–482, 2002.
  * Bica et al. (2020) Bica, I., Jordon, J., and van der Schaar, M.  Estimating the effects of continuous-valued interventions using generative adversarial networks.  _Advances in Neural Information Processing Systems_ , 33:16434–16445, 2020.
  * Chao et al. (2023) Chao, P., Blöbaum, P., and Kasiviswanathan, S. P.  Interventional and counterfactual inference with diffusion models.  _arXiv preprint arXiv:2302.00860_ , 2023.
  * Choe et al. (2023a) Choe, S. K., Mehta, S. V., Ahn, H., Neiswanger, W., Xie, P., Strubell, E., and Xing, E.  Making scalable meta learning practical.  _arXiv preprint arXiv:2310.05674_ , 2023a.
  * Choe et al. (2023b) Choe, S. K., Neiswanger, W., Xie, P., and Xing, E.  Betty: An automatic differentiation library for multilevel optimization.  In _The Eleventh International Conference on Learning Representations_ , 2023b.  URL https://openreview.net/forum?id=LV_MeMS38Q9.
  * De Brouwer (2022) De Brouwer, E.  Deep counterfactual estimation with categorical background variables.  _Advances in Neural Information Processing Systems_ , 35:35213–35225, 2022.
  * De Lara et al. (2021) De Lara, L., González-Sanz, A., Asher, N., and Loubes, J.-M.  Transport-based counterfactual models.  _arXiv preprint arXiv:2108.13025_ , 2021.
  * Du et al. (2021) Du, X., Sun, L., Duivesteijn, W., Nikolaev, A., and Pechenizkiy, M.  Adversarial balancing-based representation learning for causal effect inference with observational data.  _Data Mining and Knowledge Discovery_ , 35(4):1713–1738, 2021.
  * Hill (2011) Hill, J. L.  Bayesian nonparametric modeling for causal inference.  _Journal of Computational and Graphical Statistics_ , 20(1), 2011.
  * Javaloy et al. (2023) Javaloy, A., Sánchez-Martín, P., and Valera, I.  Causal normalizing flows: from theory to practice.  _arXiv preprint arXiv:2306.05415_ , 2023.
  * Johansson et al. (2016) Johansson, F., Shalit, U., and Sontag, D.  Learning representations for counterfactual inference.  In _International conference on machine learning_ , pp. 3020–3029. PMLR, 2016.
  * Khemakhem et al. (2021) Khemakhem, I., Monti, R., Leech, R., and Hyvarinen, A.  Causal autoregressive flows.  In _International conference on artificial intelligence and statistics_ , pp. 3520–3528. PMLR, 2021.
  * Koenker & Hallock (2001) Koenker, R. and Hallock, K. F.  Quantile regression.  _Journal of economic perspectives_ , 15(4):143–156, 2001.
  * Lake et al. (2019) Lake, B. M., Salakhutdinov, R., and Tenenbaum, J. B.  The omniglot challenge: a 3-year progress report.  _Current Opinion in Behavioral Sciences_ , 29:97–104, 2019.
  * LeCun et al. (2010) LeCun, Y., Cortes, C., and Burges, C.  Mnist handwritten digit database.  _ATT Labs [Online]. Available: http://yann.lecun.com/exdb/mnist_ , 2, 2010.
  * Li & Yao (2022) Li, X. and Yao, L.  Contrastive individual treatment effects estimation.  In _2022 IEEE International Conference on Data Mining (ICDM)_ , pp. 1053–1058. IEEE, 2022.
  * Lu et al. (2020a) Lu, C., Huang, B., Wang, K., Hernández-Lobato, J. M., Zhang, K., and Schölkopf, B.  Sample-efficient reinforcement learning via counterfactual-based data augmentation.  _arXiv preprint arXiv:2012.09092_ , 2020a.
  * Lu et al. (2020b) Lu, D., Tao, C., Chen, J., Li, F., Guo, F., and Carin, L.  Reconsidering generative objectives for counterfactual reasoning.  _Advances in Neural Information Processing Systems_ , 33:21539–21553, 2020b.
  * Melnychuk et al. (2022) Melnychuk, V., Frauen, D., and Feuerriegel, S.  Causal transformer for estimating counterfactual outcomes.  In _International Conference on Machine Learning_ , pp. 15293–15329. PMLR, 2022.
  * Mohri et al. (2018) Mohri, M., Rostamizadeh, A., and Talwalkar, A.  _Foundations of machine learning_.  MIT press, 2018.
  * Nasr-Esfahany & Kiciman (2023) Nasr-Esfahany, A. and Kiciman, E.  Counterfactual (non-) identifiability of learned structural causal models.  _arXiv preprint arXiv:2301.09031_ , 2023.
  * Nasr-Esfahany et al. (2023) Nasr-Esfahany, A., Alizadeh, M., and Shah, D.  Counterfactual identifiability of bijective causal models.  _arXiv preprint arXiv:2302.02228_ , 2023.
  * Pawlowski et al. (2020) Pawlowski, N., Coelho de Castro, D., and Glocker, B.  Deep structural causal models for tractable counterfactual inference.  _Advances in Neural Information Processing Systems_ , 33:857–869, 2020.
  * Pearl (2000) Pearl, J.  _Causality: Models, Reasoning and Inference_.  Cambridge University Press, 2000.
  * Pearl et al. (2016) Pearl, J., Glymour, M., and Jewell, N. P.  _Causal inference in statistics: A primer_.  John Wiley Sons, 2016.
  * Powell (2020) Powell, D.  Quantile treatment effects in the presence of covariates.  _Review of Economics and Statistics_ , 102(5):994–1005, 2020.
  * Ribeiro et al. (2023) Ribeiro, F. D. S., Xia, T., Monteiro, M., Pawlowski, N., and Glocker, B.  High fidelity image counterfactuals with probabilistic causal models.  _arXiv preprint arXiv:2306.15764_ , 2023.
  * Sanchez & Tsaftaris (2022) Sanchez, P. and Tsaftaris, S. A.  Diffusion causal models for counterfactual estimation.  _arXiv preprint arXiv:2202.10166_ , 2022.
  * Song et al. (2020) Song, J., Meng, C., and Ermon, S.  Denoising diffusion implicit models.  _arXiv preprint arXiv:2010.02502_ , 2020.
  * Tagasovska & Lopez-Paz (2019) Tagasovska, N. and Lopez-Paz, D.  Single-model uncertainties for deep learning.  _Advances in Neural Information Processing Systems_ , 32, 2019.
  * Takeuchi et al. (2006) Takeuchi, I., Le, Q., Sears, T., Smola, A., et al.  Nonparametric quantile estimation.  2006\.
  * Takeuchi et al. (2009) Takeuchi, I., Nomura, K., and Kanamori, T.  Nonparametric conditional density estimation using piecewise-linear solution path of kernel quantile regression.  _Neural Computation_ , 21(2):533–559, 2009.
  * Xia et al. (2021) Xia, K., Lee, K.-Z., Bengio, Y., and Bareinboim, E.  The causal-neural connection: Expressiveness, learnability, and inference.  _Advances in Neural Information Processing Systems_ , 34:10823–10836, 2021.
  * Xia et al. (2022) Xia, K., Pan, Y., and Bareinboim, E.  Neural causal models for counterfactual identification and estimation.  _arXiv preprint arXiv:2210.00035_ , 2022.
  * Xie et al. (2020) Xie, Y., Cotton, C., and Zhu, Y.  Multiply robust estimation of causal quantile treatment effects.  _Statistics in Medicine_ , 39(28):4238–4251, 2020.
  * Yao et al. (2018) Yao, L., Li, S., Li, Y., Huai, M., Gao, J., and Zhang, A.  Representation learning for treatment effect estimation from observational data.  In _Neural Information Processing Systems_ , 2018.
  * Yoon et al. (2018) Yoon, J., Jordon, J., and Van Der Schaar, M.  Ganite: Estimation of individualized treatment effects using generative adversarial nets.  In _International conference on learning representations_ , 2018.
  * Zhang et al. (2015) Zhang, K., Wang, Z., Zhang, J., and Schölkopf, B.  On estimation of functional causal models: general results and application to the post-nonlinear causal model.  _ACM Transactions on Intelligent Systems and Technology (TIST)_ , 7(2):1–22, 2015.
  * Zhou et al. (2022) Zhou, G., Yao, L., Xu, X., Wang, C., and Zhu, L.  Cycle-balanced representation learning for counterfactual inference.  In _Proceedings of the 2022 SIAM International Conference on Data Mining (SDM)_ , pp. 442–450. SIAM, 2022.
  * Zhou et al. (2021a) Zhou, T., Carson IV, W. E., and Carlson, D.  Estimating potential outcome distributions with collaborating causal networks.  _arXiv preprint arXiv:2110.01664_ , 2021a.
  * Zhou et al. (2021b) Zhou, T., Li, Y., Wu, Y., and Carlson, D.  Estimating uncertainty intervals from collaborating networks.  _The Journal of Machine Learning Research_ , 22(1):11645–11691, 2021b.

##  Appendix A Other Related Work

Individual treatment effect (ITE) estimation is also closely related to the
counterfactual inference problem. A major difference is that ITE focuses on
the effect of treatment and the predictions of both factual and counterfactual
outcomes on unseen samples. GANITE (Yoon et al., 2018) first learns a
counterfactual generator with GAN by matching the joint distribution of
observed covariate and outcome variables, and then it generates a dataset by
feeding different treatment values and random noises and learns an ITE
generator to predict the factual and counterfactual outcomes. CFRNet
(Johansson et al., 2016) models ITE as a domain adaptation problem where there
is a distribution shift between effects under different treatment and match
the marginal distributions of representations under different treatments in
the representation space. ABCEI (Du et al., 2021) proposes to use adversarial
learning to balance the representations from treatment and control groups.
CBRE (Zhou et al., 2022) proposes to use cycle consistency to preserve the
semantics of the representations from two groups. Based on GANITE, SCIGAN
(Bica et al., 2020) further proposes a hierarchical discriminator to learn the
counterfactual generator when interventions are continuous, e.g., the dosage
of the treatment. SITE (Yao et al., 2018) uses propensity score to select
positive and negative pairs and proposes to minimize the middle point distance
to preserve the relationships in the representation space. Based on SITE, CITE
(Li & Yao, 2022) employs contrastive learning to preserve the relationships.
BV-NICE (Lu et al., 2020b) models the generation process as a latent variable
model where a confounder causes the treatment, covariate, and outcomes and
addresses the covariate imbalance with adversarial training. The goal of (Zhou
et al., 2021b) is to estimate uncertainty intervals by learning two networks
in an adversarial manner, where one is to estimate CDF and the other is to
estimate the quantile. Later, (Zhou et al., 2021a) extends (Zhou et al.,
2021b) to the ITE task. (Xie et al., 2020; Powell, 2020) propose ways to
estimate the quantile treatment effects unlike the average treatment effect.
The quantile treatment effect is measured on all samples with same value of
the quantile. On the contrary, our method learns different quantile for
different individuals and use the quantile to represent the property of the
conditional distribution.

##  Appendix B Proof of Theorem 1

###### Proof.

Note $Y=f(X,Z,E)$ can be equivalently represented as $f_{1}(f_{2}(X,Z),g(E))$,
and we denote $g(E)$ by $\tilde{E}$. We know that without further restrictions
on the function class of $f$, the causal model $f$ and the probabilistic
distribution $p(\tilde{E})$ are not identifiable (Zhang et al., 2015). Denote
by $f^{i}$ and $p^{i}(\tilde{E})$ as one solution, and we will see that the
counterfactual outcome actually does not depend on the index $i$; that is, it
is independent of which $f^{i}$ and $P^{i}(\tilde{E}_{t+1})$ we choose. Given
observed evidence $(X=x,Y=y,Z=z)$, because $f^{i}$ is strictly monotonic in
$\tilde{E}^{i}$, we can determine its value $\tilde{e}^{i}$, with
$\tilde{e}^{i}={f_{x,z}^{i}}^{-1}(y)$. Then, we can determine the value of the
cumulative distribution function of $\tilde{E}^{i}$ at $\tilde{e}^{i}$,
denoted by $\tau^{i}$.

Without loss of generality, we first show the case where $f^{i}$ is strictly
increasing w.r.t. $\tilde{E}^{i}$. Because $f$ is strictly increasing in
$\tilde{E}$ and $y=f^{i}(x,z,\tilde{e}^{i})$, $y$ is the $\tau^{i}$-th
quantile of $P(Y|X=x,Z=z)$. Then it is obvious that since $y$ and
$P(Y|X=x,Z=z)$ are determined, the value of $\tau^{i}$ is independent of the
index $i$, that is, it is identifiable. Thus, below, we will use $\tau$,
instead of $\tau^{i}$.

Since $g(E)\mathchoice{\mathrel{\hbox to0.0pt{$\displaystyle\perp$\hss}\mkern
2.0mu{\displaystyle\perp}}}{\mathrel{\hbox
to0.0pt{$\textstyle\perp$\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptstyle\perp$\hss}\mkern
2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptscriptstyle\perp$\hss}\mkern
2.0mu{\scriptscriptstyle\perp}}}(X;Z)$, when doing interventions on $X$, the
value $\tilde{e}^{i}$ will not change, as well as $e^{i}$. Hence, the
counterfactual outcome $Y_{X=x^{\prime}}|X=x,Y=y,Z=z$ can be calculated as
$f^{i}(X=x^{\prime},Z=z,E=e^{i})$, and such equivalence can be directly
derived following Pearl’s three-step procedure. Because $\tilde{e}^{i}$ does
not change after the intervention, the counterfactual outcome
$Y_{X=x^{\prime}}|X=x,Y=y,Z=z$ is the $\tau$-quantile of the conditional
distribution $P(Y|X=x^{\prime},Z=z)$. This quantile exists and it depends only
on the conditional distribution $P(Y|X=x^{\prime},Z=z)$, but not the chosen
function $f^{i}$ and $P^{i}(\tilde{E})$.

Therefore, the counterfactual outcome $Y_{X=x^{\prime}}|X=x,Y=y,Z=z$
corresponds to the $\tau$-th quantile of the conditional distribution
$P(Y|X=x^{\prime},Z=z)$, where $Y=y$ represents the $\tau$-th quantile of
$P(Y|X=x,Z=z)$. ∎

##  Appendix C Proof of Theorem 2

Theorem 2 shows that the generalization error is bounded if we minimize the
empirical loss, suggesting that our method learns a good quantile estimator
using finite training samples. As a consequence (of this generalization bound
and above identifiability theorem), our method is able to perform reliable
counterfactual inference given finite training samples.

In this section, we use $f$ to represent $\mu$. Further, for simplicity, we
ignore $z$ in the $f$ function, which doesn’t affect the generalization bound
since the concatenation of $z$ and $x$ can be treated as a single input.

We first give the following generalization error bound derived by Bartlett &
Mendelson (2002).

######  Theorem 3.

Let $F$ be a hypothesis class mapping from $\mathcal{X}$ to $[0,b]$. Let
$\\{x_{i}\\}_{i=1}^{N}$ be training samples with size $N$. Then, for any
$\delta>0$, with probability at least $1-\delta$, the following holds for all
$f\in F$:

$\displaystyle\mathbb{E}_{x}[f(x)]\leq\frac{1}{N}\sum_{i=1}^{N}f(x_{i})+2\mathfrak{R}(F)+b\sqrt{\frac{\log(1/\delta)}{2N}}.$

Inspired by the above theorem, we can derive a generalization error bound for
counterfactual inference by nonlinear quantile regression.

######  Theorem 4.

Let $(\hat{\tau},\hat{f}_{\hat{\tau}})\in(\mathfrak{T},F)$ be the optimization
solution. Let the loss function $l_{\tau}$ be upper bounded by $b$. Then, for
any $\delta>0$, with probability at least $1-\delta$, we have

$\displaystyle\mathbb{E}_{(x,y)}[l_{\hat{\tau}}(f^{*}(x)-\hat{f}(x))]\leq\frac{1}{N}\sum_{i=1}^{N}l_{\hat{\tau}}(f^{*}(x_{i})-\hat{f}(x_{i}))+2\mathfrak{R}(\mathfrak{T},F)+b\sqrt{\frac{\log(1/\delta)}{2N}},$

where

$\displaystyle\mathfrak{R}(\mathfrak{T},F)=\mathbb{E}_{x,\sigma}\left[\sup_{\tau\in\mathfrak{T},f\in
F}\frac{1}{N}\sum_{i=1}^{N}\sigma_{i}l_{\tau}(f^{*}(x_{i})-f(x_{i}))\right].$

Furthermore, we can derive the following result.

######  Theorem 5.

$\displaystyle\mathfrak{R}(\mathfrak{T},F)\leq
2\mathfrak{R}(F)+\frac{2b}{\sqrt{N}}.$

Proof. Let us rewrite
$l_{\tau}(f^{*}(x)-f(x))=\tau(f^{*}(x)-f(x))-1_{\\{f^{*}(x)-f(x)<0\\}}(f^{*}(x)-f(x))$,
where $1_{\\{A\\}}$ is the indicator function. We have

$\displaystyle\mathfrak{R}(\mathfrak{T},F)$
$\displaystyle=\mathbb{E}_{x,\sigma}\left[\sup_{\tau\in\mathfrak{T},f\in
F}\frac{1}{N}\sum_{i=1}^{N}\sigma_{i}l_{\tau}(f^{*}(x_{i})-f(x_{i}))\right]$
$\displaystyle=\mathbb{E}_{x,\sigma}\left[\sup_{\tau\in\mathfrak{T},f\in
F}\frac{1}{N}\sum_{i=1}^{N}\sigma_{i}(\tau(f^{*}(x_{i})-f(x_{i}))-1_{\\{f^{*}(x_{i})-f(x_{i})<0\\}}(f^{*}(x_{i})-f(x_{i})))\right]$
$\displaystyle\leq\mathbb{E}_{x,\sigma}\left[\sup_{\tau\in\mathfrak{T},f\in
F}\tau\cdot\frac{1}{N}\sum_{i=1}^{N}\sigma_{i}(f^{*}(x_{i})-f(x_{i}))\right]$
$\displaystyle+\mathbb{E}_{x,\sigma}\left[\sup_{f\in
F}\frac{1}{N}\sum_{i=1}^{N}\sigma_{i}1_{\\{f^{*}(x_{i})-f(x_{i})<0\\}}(f^{*}(x_{i})-f(x_{i}))\right]$
$\displaystyle=\mathbb{E}_{x,\sigma}\left[\sup_{f\in
F}\frac{1}{N}\sum_{i=1}^{N}\sigma_{i}(f^{*}(x_{i})-f(x_{i}))\right]$
$\displaystyle+\mathbb{E}_{x,\sigma}\left[\sup_{f\in
F}\frac{1}{N}\sum_{i=1}^{N}\sigma_{i}1_{\\{f^{*}(x_{i})-f(x_{i})<0\\}}(f^{*}(x_{i})-f(x_{i}))\right]$
$\displaystyle\leq 2\mathbb{E}_{x,\sigma}\left[\sup_{f\in
F}\frac{1}{N}\sum_{i=1}^{N}\sigma_{i}(f^{*}(x_{i})-f(x_{i}))\right]$
$\displaystyle=2\mathbb{E}_{x,\sigma}\left[\sup_{f\in
F}\frac{1}{N}\sum_{i=1}^{N}\sigma_{i}f(x_{i})\right]+2\mathbb{E}_{x,\sigma}\left[\frac{1}{N}\sum_{i=1}^{N}\sigma_{i}f^{*}(x_{i})\right]$
$\displaystyle\leq 2\mathfrak{R}(F)+\frac{2b}{\sqrt{N}}$

where the second inequality holds because that $\mathfrak{R}(\phi\circ F)\leq
L\mathfrak{R}(F)$ when $\phi:\mathbb{R}\rightarrow\mathbb{R}$ is an
$L$-Lipschitz and that $1_{x<0}x$ is $1$-Lipschitz w.r.t. $x$. The last
inequality holds because
$\mathbb{E}_{x,\sigma}\left[\frac{1}{N}\sum_{i=1}^{N}\sigma_{i}f^{*}(x_{i})\right]\leq\sqrt{\frac{b}{N}}$.
Specifically,

$\displaystyle\mathbb{E}_{x,\sigma}\left[\frac{1}{N}\sum_{i=1}^{N}\sigma_{i}f^{*}(x_{i})\right]$
$\displaystyle=\frac{1}{N}\mathbb{E}\left[\sum_{i=1}^{N}\sigma_{i}f^{*}(x_{i})\right]$
$\displaystyle\leq\frac{1}{N}\left[\mathbb{E}\left[\left(\sum_{i=1}^{N}\sigma_{i}f^{*}(x_{i})\right)^{2}\right]\right]^{1/2}$
$\displaystyle=\frac{1}{N}\left[\mathbb{E}\left[\sum_{i,j=1}^{N}\sigma_{i}\sigma_{j}f^{*}(x_{i})f^{*}(x_{j})\right]\right]^{1/2}$
$\displaystyle=\frac{1}{N}\left[\mathbb{E}\left[\sum_{i=1}^{N}(f^{*}(x_{i}))^{2}\right]\right]^{1/2}$
$\displaystyle\leq\frac{b}{\sqrt{N}}$ (4)

where the first inequality follows by Jensen’s inequality, the third equation
holds because $\\{\sigma_{i}\\}_{i=1}^{N}$ are i.i.d. and
$\mathbb{E}[\sigma_{i}\sigma_{j}]=\mathbb{E}[\sigma_{i}]\mathbb{E}[\sigma_{j}]=0$
when $i\neq j$ while $\mathbb{E}[\sigma_{i}\sigma_{i}]=1$.

Therefore, Theorem 2 can be derived by combining Theorem 4 and 5.

##  Appendix D Datasets

Datasets. Since our main result is based on the monotonicity assumption and
the observations in counterfactual scenarios are lacking in real-world, we
create following datasets.

Rotation-MNIST. We use MNIST images (LeCun et al., 2010) as $Z$ and rotation
angles as $X$. Then we randomly sample a noise value from $U[0,1]$ and adjust
the rotated images by adding values to the RGB channels. Therefore, the final
pixel values are strictly monotonic w.r.t the noise value when conditioning on
$X$ and $Z$. The training set consists of 60000 images while the testing set
consists of 10000 images. The dimension of $Z$ is $3\times 32\times 32=3072$
and the dimension of $X$ is 1, which determines the rotation angles within
[-45,45]. The dimension of $Y$ is also $3\times 32\times 32$. In
counterfactual inference, we are given a sample $X,Z,Y$, the question asks
what had the images would be if the rotation angle is set to another value. In
other words, we aim to rotate the digits in $Y$ while avoiding any other pixel
value changes.

Thick-Omniglot. We use Omniglot images (Lake et al., 2019) as $Z$ and the
thickness of characters as $X$. Then we randomly sample a noise value from
$U[0,1]$ and adjust the background darkness by multiplying the thicken
characters pixel values. Therefore, the final pixel values are strictly
monotonic w.r.t the noise value when conditioning on $X$ and $Z$. The number
of images is 19280 and We split the dataset as 80/20. We also resize the
images to 32$\times$32\. In this dataset, our goal is to adjust the thickness
by the new value of $X$ while preserving the brightness in $Y$.

Cont-Dose. We generate this dataset by setting $X$ as the treatment and $Z$ as
the age of patient and the outcome $Y$ denotes the effect. $X$ takes value
from [0,2] with space 0.1 uniformly. The numbers of training samples is 800
and the number of testing samples is 200.

Dis-Dose. We generate this dataset by setting $X$ as the treatment and $Z$ as
the age of patient and the outcome $Y$ denotes the effect. $X$ takes value
from $\\{0,1\\}$. The numbers of training samples is 800 and the number of
testing samples is 200.

We also use the semi-simulated dataset IHDP (Hill, 2011). It contains 100
splits. Each split consists of 675 training and 92 testing samples. The
dimension of $Z$ is 25 and the dimension of $X$ is 1 and the dimension of
output $Y$ is 1.

##  Appendix E Implementation

All codes are provided in the supplementary materials. Readers may refer to
the code for more implementation details.

We use betty (Choe et al., 2023b) as our bi-level optimization library. Then
we define the lower level loss as the pin-ball loss with the quantile passed
from the upper level. After optimizing the lower level net $g$ with 30
iterations, we fix the network $g$ and optimize the upper-level network $h$ on
the regression loss (We adopt MSE loss).

For toy datasets and cont-Dose dataset, we build the upper-level network as 3
layers and 1 output linear layer. The network $h$ takes as input the factual
observations $X$, $Z$ and $Y$. Then we concatenate the three variables and
feed into the model. The first layer is a linear layer transforming the input
into a hidden-dimension (200) and apply a SiLU activation function. Then we
use 2 residual blocks (with layernorm) and finally feed the output into the
linear layer with output size as 1. The final activation is Sigmoid since we
need to use $h$ to represent quantile, which is within range [0,1]. As for the
lower-level network $g$, it is mainly trained to get a regression at quantile
$\tau$, which is the output of network $h$. We first project the input
$X,Z,\tau$ into same dimensions respectively, then we feed into a linear layer
to match the dimension of $Y$.

For Dis-Dose and IHDP, $X$ is binary. Therefore, we use two networks with same
architecture as $h$. We feed the inputs $X,Z,Y$ into different sub-network
according to the value $X$. We find that a shallower network is better,
therefore, we do not use residual block here. For the network $h$, we use
3-layers (one input layer and 2 residual blocks).

For image transformation datasets, we use the convolutional neural network to
parameterize $h$ and $g$. The hidden dimension of $h$ is 32 while the hidden
dim of $g$ is 128. We downsample the images with Conv(4,2,1) and final
conv(4,1,0) to get the quantile estimations. In network $g$, we use a
symmetric network where the upsamplinmg network mimic the downsampling process
with ConvTranspose(4,2,1).

We use $n_{i}$ interested samples in training network $g$. For eacch
interested sample, we have one quantile $\tau$. Then we train the network $g$
with $n_{t}$ samples for each $\tau$. We set $n_{i}=256,n_{t}=64$ for most
cases. As for image datasets, we have to lower the values to reduce memory. So
we use $n_{i}=128,n_{t}=32$. In the upper-level problem, we use 64 (128)
samples to get the reconstruction loss. We use Adam optimizer with lr=2e-3 for
images while 1e-3 for the rest of them.

We choose the hyper-parameters based on the values on some toy datasets since
we can always generate them. Then we apply the best hyper-parameters to the
formal dataset, Dis-Dose, Cont-Dose, IHDP, Thick-Ominiglot and Rotation-MNIST.
During training, we use the reconstruction error on the training dataset as
metric to select the models.

##  Appendix F Confounding

In the main paper, we considered a simple case where $Y=X+Z+E$, where $Z\sim
U[0,1],X\sim U[0,1]$ and $e\sim N(0,1)$. We generate the latent confounder as
$C\sim U[-0.5,0.5]$. Then we add the confounder to the its children in
dfferent scenarios. For example, when we have $Z\leftarrow C\rightarrow Y$, we
add $Z=Z+C,Y=Y+C$. According to the table in the main paper, the performances
drop on the scenario $X\leftarrow C\rightarrow Y$ in this case, which is
kindly expected since we do not assume the existence of latent confounder and
$C$ influences $X$ and $Y$ now. (the code of data generation and training is
provided in the supp.)

##  Appendix G More Results on Image Transformation Datset

Learned Quantiles. Since we sample the noise values uniformly from [0,1]. The
target quantile $P(Y\leq y|X=x,Z=z)$ is the CDF of the noise. In other words,
if the noise is $a$, the target quantile should be $a$. We present the
examples of Thick-Omniglot in Fig. 10. Our estimated quantiles are not
affected by the thickness changes between the input $Z$ and the output $Y$.
The The learned quantiles are very close to the ground truth, further
demonstrating the effectiveness of our method. We also provide examples of
Rotation-MNIST in Fig. 9. Our method is able recover the quantiles of the true
noises accurately.

Z | E= 0.1 | 0.3 | 0.5 | 0.7 | 0.9 |
---|---|---|---|---|---|---
|  $\hat{\tau}=$ 0.1075 | 0.3088 | 0.5094 | 0.7032 | 0.9019 |
|  |  |  |  |  |
|  $\hat{\tau}$= 0.1034 | 0.2945 | 0.5002 | 0.6984 | 0.8972 |
|  |  |  |  |  |
|  $\hat{\tau}$= 0.0979 | 0.3069 | 0.5019 | 0.7060 | 0.9080 |
|  |  |  |  |  |
|  $\hat{\tau}$=0.1045 | 0.3026 | 0.4977 | 0.7003 | 0.8997 |
|  |  |  |  |  |
|  $\hat{\tau}$= 0.1057 | 0.3050 | 0.5024 | 0.7027 | 0.8970 |
|  |  |  |  |  |
|  $\hat{\tau}$= 0.1059 | 0.3032 | 0.5013 | 0.6956 | 0.8972 |
|  |  |  |  |  |
|  $\hat{\tau}$= 0.1059 | 0.3021 | 0.5012 | 0.7022 | 0.8967 |
|  |  |  |  |  |
|  $\hat{\tau}$= 0.1051 | 0.2991 | 0.5038 | 0.7017 | 0.8989 |
|  |  |  |  |  |
|  $\hat{\tau}$= 0.1007 | 0.2998 | 0.4979 | 0.7030 | 0.9027 |
|  |  |  |  |  |
|  $\hat{\tau}$= 0.1015 | 0.3015 | 0.4989 | 0.6986 | 0.8981 |
|  |  |  |  |  |
Figure 9: Examples of learned quantiles on Rotation-MNIST dataset. The learned quantiles $\hat{\tau}$ are very close to the true quantiles, i.e., 0.1, 0.3,0.5,0.7 and 0.9.  Z | E= 0.1 | 0.3 | 0.5 | 0.7 | 0.9
---|---|---|---|---|---
|  $\hat{\tau}=$ 0.0908 | 0.3034 | 0.5011 | 0.7062 | 0.8981
|  |  |  |  |
|  $\hat{\tau}$= 0.0934 | 0.3011 | 0.5014 | 0.7082 | 0.8929
|  |  |  |  |
|  $\hat{\tau}$= 0.0934 | 0.3055 | 0.5040 | 0.7089 | 0.8982
|  |  |  |  |
|  $\hat{\tau}$=0.0932 | 0.2901 | 0.4612 | 0.6431 | 0.8239
|  |  |  |  |
|  $\hat{\tau}$= 0.1028 | 0.2972 | 0.4956 | 0.7003 | 0.8961
|  |  |  |  |
|  $\hat{\tau}$= 0.0863 | 0.3079 | 0.4936 | 0.6877 | 0.8650
|  |  |  |  |
|  $\hat{\tau}$= 0.1076 | 0.3023 | 0.4953 | 0.6999 | 0.8960
|  |  |  |  |
|  $\hat{\tau}$= 0.0938 | 0.2995 | 0.4868 | 0.6837 | 0.8648
|  |  |  |  |
|  $\hat{\tau}$= 0.0991 | 0.3050 | 0.4993 | 0.7036 | 0.8947
|  |  |  |  |
|  $\hat{\tau}$= 0.0946 | 0.2982 | 0.5014 | 0.7107 | 0.8964
|  |  |  |  |
Figure 10: Examples of learned quantiles on Thick-Omniglot dataset. The
learned quantiles $\hat{\tau}$ are very close to the true quantiles, i.e.,
0.1, 0.3,0.5,0.7 and 0.9. Each row displays the images after adjusting
thickness by $X$ and the brightness with noise value. The noise determines the
brightness of the images. From the $Z$ and the output image, our method is not
affected by the changing thickness between $Z$ and $Y$ and able to capture the
brightness changes. The accurate estimation shows the effectiveness of our
approach even when $Z$ and $Y$ are high-dimensional.

Counterfactual Prediction Results We present more visual examples of Rotation-
MNIST in Fig. 4 and 5. We can see that the CFQP (De Brouwer, 2022) learns the
rotation when changing the value of $X$. However, the rotations results is
unsatisfactory. In particular, the background color of the rotated digits are
different from the ground truth $Y_{X=x^{\prime}}$ and the $Y$. This is
unwanted because $X$ only affects the rotation. It indicates that CFQP fails
to capture the uniqueness (noise) of the sample and leads to unnecessary
changes. As for BGM (Nasr-Esfahany & Kiciman, 2023), the images are almost
identical to the $Y$. The reason is that the BGM employs conditional spline
flow to mimic the strictly monotonic generation process. However, it is known
that the expressive power of flow, especially such strictly monotonic flow,
are limited compared to the unconstrained neural networks. It fails to
transform the high-dimensional input $Y$ into gaussian noise and therfore,
fails to utilize the conditions ($X$ and $Z$). As a consequence, when we
change the value of $X$, BGM ignores the condition and generates almost
identical images to the input $Y$. DCM (Chao et al., 2023) models the
generation process with conditional diffusion model and recovers the noise by
inverting the diffusion process. However, since the true noise is usually not
Gaussian, the transformed noise may not be meaningful to counterfactual
inference. Therefore, the counterfactual predictions are not accurate. In
contrast to above baseline methods, our approach learns the correct rotation
as well as preserving the correct color in the factual observation $Y$.

We also present examples of Thick-Omniglot in Fig. 11. The noises values
determines the darkness of the image. We obvserve that CFQP (De Brouwer, 2022)
fails to preserve the darkness of the image $Y$. Sometimes the predictions are
darker or much brighter. DCM (Chao et al., 2023) struggles to adjust the
thickness of the images. In contrast, our method learns to preserve the
brightness in the factual observation $Y$ while changing the thickness of the
digits.

Z | Y | CFQP | BGM | DCM | Ours | Truth
---|---|---|---|---|---|---
|  |  |  |  |  |
|  |  |  |  |  |
|  |  |  |  |  |
|  |  |  |  |  |
|  |  |  |  |  |
|  |  |  |  |  |
|  |  |  |  |  |
|  |  |  |  |  |
|  |  |  |  |  |
|  |  |  |  |  |
Table 4: Examples of counterfactual predictions on Rotation-MNIST.  Z | Y | CFQP | BGM | DCM | Ours | Truth
---|---|---|---|---|---|---
|  |  |  |  |  |
|  |  |  |  |  |
|  |  |  |  |  |
|  |  |  |  |  |
|  |  |  |  |  |
|  |  |  |  |  |
|  |  |  |  |  |
|  |  |  |  |  |
|  |  |  |  |  |
|  |  |  |  |  |
Table 5: Examples of counterfactual predictions on Rotation-MNIST. CFQP fails to preserve the brightness of image $Y$ in their counterfactual predictions. BGM fails to transform the image $Y$ given the limited expressive power of the used spline flow. DCM genrates unrealistic digits.  Z | Y | CFQP | BGM | DCM | Ours | Truth
---|---|---|---|---|---|---
|  |  |  |  |  |
|  |  |  |  |  |
|  |  |  |  |  |
|  |  |  |  |  |
|  |  |  |  |  |
|  |  |  |  |  |
|  |  |  |  |  |
|  |  |  |  |  |
|  |  |  |  |  |
|  |  |  |  |  |
Figure 11: Examples of counterfactual predictions on Thick-Omniglot dataset.