capbtabboxtable[][]

# VACA: Design of Variational Graph Autoencoders for Interventional and
Counterfactual Queries

Pablo Sanchez-Martin
Max Planck Institute for Intelligent Systems
Tübingen, Germany
<EMAIL_ADDRESS>
&Miriam Rateike*
Max Planck Institute for Intelligent Systems
Tübingen, Germany
<EMAIL_ADDRESS>
&Isabel Valera
Saarland University
Saarbücken, Germany
<EMAIL_ADDRESS>
Equal contribution

###### Abstract

In this paper, we introduce VACA, a novel class of variational graph
autoencoders for causal inference in the absence of hidden confounders, when
only observational data and the causal graph are available. Without making any
parametric assumptions, VACA mimics the necessary properties of a Structural
Causal Model (SCM) to provide a flexible and practical framework for
approximating interventions (_do-operator_) and abduction-action-prediction
steps. As a result, and as shown by our empirical results, VACA accurately
approximates the interventional and counterfactual distributions on diverse
SCMs. Finally, we apply VACA to evaluate counterfactual fairness in fair
classification problems, as well as to learn fair classifiers without
compromising performance.

_K_ eywords Graph Neural Network $\cdot$ Causality $\cdot$ Counterfactual
$\cdot$ Intervention $\cdot$ Variational Autoencoder

##  1 Introduction

Graph Neural Networks (GNNs) are a powerful tool for graph representation
learning and have been proven to excel in practical complex problems like
neural machine translation [1], traffic forecasting [5, 47], or drug discovery
[11].

In this work, we investigate to which extent the inductive bias of
GNNs–encoding the causal graph information–can be exploited to answer
interventional and counterfactual queries. More specifically, to approximate
the interventional and counterfactual distributions induced by interventions
on a casual model. To this end, we assume i) causal sufficiency–i.e., absence
of hidden confounders; and, ii) access to observational data and the true
causal graph. We stress that the causal graph can often be inferred from
expert knowledge [52] or via one of the approaches for causal discovery [12,
42]. With this analysis we aim to complement the concurrent line of research
that theoretically studies the use of Neural Networks (NN) [45], and more
recently GNNs [49], for causal inference.

To this end, we describe the architectural design conditions that a
variational graph autoencoder (VGAE)–as a density estimator that leverages a
priori graph structure–must fulfill so that it can approximate causal
interventions (_do-operator_) and abduction-action-prediction steps [33]. The
resulting Variational Causal Graph Autoencoder, referred to as VACA, enables
_approximating_ the observational, interventional and counterfactual
distributions induced by a causal model with unknown structural equations. We
remark that parametric assumptions on the structural causal equations are in
general not testable, may thus not hold in practice [34] and may lead to
inaccurate results, if misspecified. VACA addresses this limitation by
including uncertainty, i.e., a probabilistic model, in the estimation of the
causal-parent relationships.

We show in extensive synthetic experiments that VACA outperforms competing
methods [16, 17] on complex datasets at estimating not only the mean of the
interventional/counterfactual distribution (as in previous work), but also the
overall distribution (measured in terms of Maximum Mean Discrepancy [13]).
Finally, we show a practical use-case in which VACA is used to assess
counterfactual fairness of different classifiers trained on the real-world
German Credit dataset [7], as well as to learn counterfactually fair
classifiers without compromising performance.

###  1.1 Related Work

Deep generative models are enjoying increasing attention for causal queries in
complex data [26, 28]. Existing approaches for causal inference focus on i)
estimating the Average Treatment Effect (ATE)–a specific type of group-level
causal queries–by assuming a fixed causal graph that includes a treatment
variable [18, 25, 36, 39, 43, 50]; ii) discovering and intervening on the
causal latent structure of the (e.g., image) data [18, 28, 30, 40, 46]; or
iii) addressing interventional and/or counterfactual queries by fitting a
conditional model for each observed variable given its causal parents [10, 16,
21, 29, 31].

Within the scope of causality, GNNs have predominantly been used for causal
discovery [48, 51] and only very recently, concurrent with us, exploited to
answer interventional queries [49].

Khemakhem et al. [17] propose CAREFL, an autoregressive normalizing flow for
both causal discovery and inference. The authors focus on (multi-dimensional)
bi-variate graphs, but their approach can be extended to more general directed
acyclic graphs (DAGs) using e.g., neuronal spline flows [8]. However, causal
assumptions in a graph are modeled not only by the direction of edges, but
also the absence of edges [32]. For the task of causal inference CAREFL is
unable to exploit the absence of edges fully as it reduces a causal graph to
its causal ordering (which may not be unique). Further, the authors only
evaluate interventions in root nodes (which reduces to conditioning on the
intervened-upon variable).

Karimi et al. [16] answer interventional queries by fitting a conditional
variational autoencoder (CVAE) to each conditional in the Markov factorization
implied by the causal graph. As each observed variable is independently
fitted, the mismatch between the true and generated distribution can cause
errors that propagate to the distribution of its descendants. This can be
problematic, especially for long causal paths. Pawlowski et al. [31] propose
an approach similar to Karimi et al. [16], and additionally propose an
approach based on normalizing flows to approximate the causal parent-child
effect.

In contrast, VACA leverages i) GNNs to encode the causal graph information
(inductive bias), ii) the GNN message passing algorithm to approximate the
effect of interventions (do-operator [33]) in the causal graph, and iii)
jointly optimizes the observational distribution for all observed variables to
avoid error propagation along the Markov factorization. We thoroughly evaluate
the performance of VACA and compare it with related work, at approximating
both interventional and counterfactual distributions induced by interventions
on both root and non-root nodes in a wide variety of causal models.

##  2 Background

In this section, we first provide a brief overview on SCMs and then introduce
the main building block of VACA, i.e., variational graph autoencoders.

###  2.1 Structural causal models

An SCM ${\mathcal{M}}=(p({\mathbf{U}}),{{{\tilde{\mathbf{F}}}}})$ determines
how a set of $d$ endogenous (observed) random variables
${\mathbf{X}}:=\left\\{X_{1},\dots X_{d}\right\\}$ is generated from a set of
exogenous (unobserved) random variables ${\mathbf{U}}:=\left\\{U_{1},\dots
U_{d}\right\\}$ with prior distribution $p({\mathbf{U}})$ via the set of
_structural equations_
${{{\tilde{\mathbf{F}}}}}=\\{X_{i}:={{{\tilde{f}}}}_{i}\left({\mathbf{X}}_{{{{\mathrm{pa}(i)}}}},U_{i}\right)\\}_{i=1}^{d}$.
Here ${\mathbf{X}}_{{{{\mathrm{pa}(i)}}}}$ refers to the set of variables
directly causing $X_{i}$, i.e., parents of $i$. Similarly to [16, 17, 32], we
consider SCMs that are associated with a directed acyclic causal graph
(although Section 4.4 relaxes this assumption). We here denote the causal
graph by ${\mathcal{G}}:=({\mathcal{V}},{\mathcal{E}})$, where each node
$i\in{\mathcal{V}}$ corresponds to an endogenous variable $X_{i}$. The set of
directed edges $(j,i)\in{\mathcal{E}}$ represent the causal parent-child
relationship between endogenous variables [32], i.e. $X_{j}$ is a parent of
$X_{i}$. ${\mathcal{E}}$ can be represented by the adjacency matrix
${\mathbf{A}}\in\\{0,1\\}^{d\times d}$, such that $A_{ij}=1$ if
$(j,i)\in{\mathcal{E}}$ and $A_{ij}=0$, otherwise. We also define the set of
neighbors, a.k.a. parents, of node $i$ as
${{{\mathrm{pa}(i)}}}={\mathcal{N}}_{i}=\\{j\\}_{(j,i)\in{\mathcal{E}}}$ and
${{{\mathrm{pa}^{*}}(i)}}:={{{\mathrm{pa}(i)}}}\cup{i}$. Given an SCM, there
are two types of causal queries of general interest: interventional queries,
e.g., “What would happen to the population ${\mathbf{X}}$, if variable $X_{i}$
would be set to a fixed value $\alpha$?”, and counterfactual queries,
e.g.,“What would have happened to a specific factual sample $\mathbf{x}^{F}$,
had $X_{i}$ been set to a value $\alpha$?”. In more detail, _interventional
queries_ aim to evaluate the effect at the population level (rung 2) of a
specific intervention on, or equivalently manipulations of, a subset of the
endogenous variables $\mathcal{I}\subseteq[d]:=\\{1,\dots,d\\}$. Interventions
on an SCM ${\mathcal{M}}$ are often represented with the do-operator
$do(X_{i}=\alpha_{i})$ [33] and lead to a modified SCM
${\mathcal{M}}^{{{\mathcal{I}}}}$ which induces a new distribution over the
set of endogenous variables $p({\mathbf{X}}\mid do(X_{i}=\alpha_{i}))$, which
is referred to as the interventional distribution. In ${\mathcal{G}}$ an
intervention removes incoming edges to node $i$ and sets $X_{i}=\alpha$ (see
Figure 1(b)). A counterfactual query for a given factual instance
$\mathbf{x}^{F}$ aims to estimate what would have happened had
${{{\mathbf{X}_{{{\mathcal{I}}}}}}}$ instead taken value
${{\boldsymbol{\alpha}}}$. This effect is captured by the counterfactual
distribution
$p(\mathbf{x}^{CF}\mid\mathbf{x}^{F},do(X_{{\mathcal{I}}}={{\boldsymbol{\alpha}}}))$,
which can be computed using the _abduction-action-prediction_ procedure by
Pearl [33]. Refer to Section 3 for further details on the computation of the
interventional and counterfactual distributions within our framework.

$X_{1}$$X_{3}$$X_{2}$$U_{1}$$U_{2}$$U_{3}$

(a) ${\mathcal{G}}$ without intervention

$X_{1}$$X_{3}$$\alpha$ •$U_{1}$$U_{3}$

(b) ${\mathcal{G}}$ with intervention

Figure 1:  triangle SCM
${\mathcal{M}}=\\{p({\mathbf{U}}),{{{\tilde{\mathbf{F}}}}}\\}$,
${\mathbf{U}}\sim p({\mathbf{U}})$ with $d=|{\mathbf{X}}|=3$ endogenous
variables where $X_{1}:={{{\tilde{f}}}}_{1}(U_{1})$,
$X_{2}:={{{\tilde{f}}}}_{2}(X_{1},U_{2})$,
$X_{3}:={{{\tilde{f}}}}_{3}(X_{1},X_{2},U_{3})$ with (a) the corresponding
causal graph ${\mathcal{G}}$ and (b) the causal graph corresponding to
${\mathcal{M}}^{{{\mathcal{I}}}}$ after intervention $do(X_{2}=\alpha)$. Blue
(red) arrows highlight the direct (indirect) causal path from $X_{1}$ to
$X_{3}$ (via $X_{2}$).

###  2.2 Variational Graph Autoencoder and Graph Neural Networks

Variational Autoencoders (VAEs). VAEs [19] are powerful latent variable models
based on neural networks (NNs) for jointly i) learning expressive density
estimators $p({\mathbf{X}})\approx\int
p_{\theta}({\mathbf{X}}\mid{\mathbf{Z}})p({\mathbf{Z}})d{\mathbf{Z}}$, where
the likelihood function (a.k.a. _decoder_) is parameterized using a NN with
parameters $\theta$, and ii) performing approximate posterior inference over
the latent variables ${\mathbf{Z}}$ made possible via a variational
distribution (a.k.a. _encoder_)
$q_{\phi}\left({\mathbf{Z}}\mid{\mathbf{X}}\right)$ parameterized using a NN
with parameters ${\phi}$. The parameters $\theta$ and $\phi$ can be learned by
maximizing a lower bound on the log-evidence [2, 27, 35, 41].

Variational Graph Autoencoders (VGAEs). Kipf and Welling [20] extend VAEs to
account for prior graph structure information on the data [48]. VGAEs define a
(potentially multidimensional) latent variable $Z_{i}$ per observed variable
$X_{i}$, i.e., ${\mathbf{Z}}:=\\{Z_{1},\dots,Z_{d}\\}$. Additionally, VGAEs
rely on an adjacency matrix ${\mathbf{A}}$, which is used by two GNNs, one for
the encoder and one for the decoder, to enforce structure on the posterior
approximation $q_{\phi}({\mathbf{Z}}\mid{\mathbf{X}},{\mathbf{A}})$ and the
likelihood $p_{\theta}({\mathbf{X}}\mid{\mathbf{Z}},{\mathbf{A}})$. Hence,
${\mathbf{A}}$–given as prior–determines which variables $X_{i}$ influence
$Z_{j}$ $\forall i,j\in[d]$.

Graph Neural Networks (GNNs). In its most general form, a GNN is a composition
of message passing layers [11], where each layer updates the state of each
node in $\mathcal{G}$. In particular, the state of node $i$ at the output of
layer $l$, i.e., $\mathbf{h}_{i}^{l}$, is specified as:

$\displaystyle\mathbf{h}_{i}^{l}$
$\displaystyle=f^{u}\left(\mathbf{h}_{i}^{l-1},f^{a}\left(\\{\mathbf{m}_{ij}^{l}\\}_{j\in{{\mathcal{N}_{i}}}}\right);\theta^{l}_{u}\right).$
(1)

First, node $i$ receives a message
$\mathbf{m}_{ij}^{l}=f^{m}(\mathbf{h}_{i}^{l-1},\mathbf{h}_{j}^{l-1};\theta_{m}^{l})$
from each of its neighbors $j\in{\mathcal{N}}_{i}$. Then, these messages are
aggregated via $f^{a}$. Finally, $\mathbf{h}_{i}^{l}$ is computed as a
function $f^{u}$ of the node’s previous state $\mathbf{h}_{i}^{l-1}$ and the
aggregated message. Note, if a GNN has ${N_{h}}$ hidden layers, then the
output for node $i$ depends not only on its direct neighbors
${{\mathcal{N}_{i}}}$, but also on its neighbors up to order $N_{h}+1$ (hops).
For example, if $N_{h}=0$ ($N_{h}=1$) then the output for each node only
depend on its direct neighbors, i.e., _parents_ (2-hop neighbors, i.e.,
_grand-parents_). For a detailed description of GNNs, please refer to Appendix
A.

##  3 Observational, interventional and counterfactual distributions

In this section, we introduce the observational, interventional and
counterfactual distributions (triggered by any intervention of the form
$do({\mathbf{X}}_{{\mathcal{I}}}={{\boldsymbol{\alpha}}})$) that are induced
by an SCM $\mathcal{M}=\\{p({\mathbf{U}}),{{{\tilde{\mathbf{F}}}}}\\}$.
Specifically, we summarize the main properties of an SCM that will allow us to
propose a novel class of VGAEs, namely VACA, to compute accurate estimates of
these distributions using observational data and a known causal graph. To this
end, we assume the absence of hidden confounders, i.e., we assume that
$p({\mathbf{U}})=\prod_{i=1}^{d}p(U_{i})$.

Observational distribution. The SCM $\mathcal{M}$ determines the observational
distribution $p({\mathbf{X}})$ over the set of endogenous variables
${\mathbf{X}}=\left\\{X_{1},\dots X_{d}\right\\}$, which satisfies causal
factorization [38], i.e.,
$p({\mathbf{X}})=\prod_{i=1}^{d}p(X_{i}\mid{\mathbf{X}}_{\mathrm{pa}(i)}).$
That is, after marginalizing out the exogenous variables ${\mathbf{U}}$, the
distribution of each endogenous variable $X_{i}$ depends only on its parents,
i.e., ${\mathbf{X}}_{{{{\mathrm{pa}(i)}}}}$.

The observational distribution can alternatively be written only in terms of
the exogenous variables ${\mathbf{U}}$ as

$\displaystyle p({\mathbf{X}})={\mathbf{F}}_{\\#}[p({\mathbf{U}})],$ (2)

i.e., $p(\mathbf{X})$ is the pushforward of $P(\mathbf{U})$ through
$\mathbf{F}$. Here ${\mathbf{F}}:{\mathbf{U}}\rightarrow{\mathbf{X}}$
corresponds to the set of structural equations, which directly transform the
exogenous variables ${\mathbf{U}}$ into the endogenous variables
${\mathbf{X}}$. This is equivalent to ${{{\tilde{\mathbf{F}}}}}$, which takes
as input both the exogenous variable and the parent (endogenous) variables of
a target endogenous variables to compute its value.

Let us denote by ${{{\mathrm{an}(i)}}}$ the set of indexes of the ancestors of
$i$, and ${{{\mathrm{an}^{*}}(i)}}:={{{\mathrm{an}(i)}}}\cup\\{i\\}$. Then,
the causal factorization induced by ${\mathcal{M}}$ leads to the following
property of ${\mathbf{F}}({\mathbf{U}})$:

######  Property 1.

Each endogenous variable $X_{i}$ can be expressed as a function of its
exogenous variable $U_{i}$ and the ones of all its causal ancestors, i.e.,
${\mathbf{F}}({\mathbf{U}})$$=$$\\{X_{i}=f_{i}(\\{U_{j}\\}_{j\in{{{\mathrm{an}^{*}}(i)}}})\\}$.
This, together with the causal sufficiency assumption, implies that $X_{i}$ is
statistically independent of $U_{j},\forall j\notin{{{\mathrm{an}^{*}}(i)}}$.

Interventional distribution. As stated in Section 2.1, interventions on a set
of variables ${{\mathcal{I}}}$ can be performed using the do-operator, which
can be seen as a mapping
$do({{{\mathbf{X}_{{{\mathcal{I}}}}}}}={{\boldsymbol{\alpha}}}):{\mathcal{M}}\mapsto{\mathcal{M}}^{{{\mathcal{I}}}}=(p({\mathbf{U}}),{{{\tilde{\mathbf{F}}^{{{\mathcal{I}}}}}}})$
where
${{{\tilde{\mathbf{F}}^{{{\mathcal{I}}}}}}}=\\{{{{\tilde{f}}}}_{i}\\}_{i\not\in{{\mathcal{I}}}}\cup\\{\alpha_{i}\\}_{i\in{{\mathcal{I}}}}$.
As above, we can represent the resulting set of intervened structural
equations as
${{{\mathbf{F}^{{{\mathcal{I}}}}}}}=\\{f_{i}\\}_{i\not\in{{\mathcal{I}}}}\cup\\{\alpha_{i}\\}_{i\in{{\mathcal{I}}}}$,
and thus write the interventional distribution as:

$\displaystyle p({\mathbf{X}}\mid
do({{{\mathbf{X}_{{{\mathcal{I}}}}}}}={{\boldsymbol{\alpha}}}))={{{\mathbf{F}^{{{\mathcal{I}}}}}}}_{\\#}[p({\mathbf{U}})].$
(3)

######  Property 2.

After an intervention
$do({{{\mathbf{X}_{{{\mathcal{I}}}}}}}={{\boldsymbol{\alpha}}})$ on
${\mathcal{M}}$, all the causal paths from $U_{j}{{{{\text{
}{\forall}}}}}j\in{{{\mathrm{an}^{*}}(i)}}$ to $X_{i}$ that include an
intervened-upon variable in ${{{\mathbf{X}_{{{\mathcal{I}}}}}}}$ (i.e., the
causal paths where ${{{\mathbf{X}_{{{\mathcal{I}}}}}}}$ is a mediator) are
severed in ${{{\mathbf{F}^{{{\mathcal{I}}}}}}}$, while the rest of causal
paths remain untouched.

The above property is illustrated in Figure 1, where we can observe that after
an intervention $do(X_{2}=\alpha)$, the indirect causal path (in red) from
$X_{1}$, and thus from $U_{1}$, to $X_{3}$ via $X_{2}$ is severed, while the
direct path (in blue) remains.

Counterfactual distribution. Assuming the SCM
$\mathcal{M}=\\{p({\mathbf{U}}),{{{\tilde{\mathbf{F}}}}}\\}$ to be known, the
following three steps defined by Pearl [32] allow to compute counterfactuals
$\mathbf{x}^{CF}$: i) Abduction: infer the values of the exogenous variables
${\mathbf{U}}$ for a factual sample $\mathbf{x}^{F}$, i.e., compute
$p({\mathbf{U}}\mid\mathbf{x}^{F})$; ii) Action: intervene with
$do({{{\mathbf{X}_{{{\mathcal{I}}}}}}}={{\boldsymbol{\alpha}}}):{\mathcal{M}}\mapsto{\mathcal{M}}^{{{\mathcal{I}}}}=(p({\mathbf{U}}),{{{\tilde{\mathbf{F}}^{{{\mathcal{I}}}}}}})$;
and iii) Prediction: use the posterior distribution
$p({\mathbf{U}}\mid\mathbf{x}^{F})$ and the new structural equations
${{{\tilde{\mathbf{F}}^{{{\mathcal{I}}}}}}}$ to compute
$p(\mathbf{x}^{CF}\mid\mathbf{x}^{F})$. The prediction step can alternatively
be computed using the new set of structural equations
${{{\mathbf{F}^{{{\mathcal{I}}}}}}}$ defined in terms of the exogenous
variables ${\mathbf{U}}$, so that we can write the counterfactual distribution
as:

$\displaystyle{p(\mathbf{x}^{CF}\mid\mathbf{x}^{F},do({{{\mathbf{X}_{{{\mathcal{I}}}}}}}={{\boldsymbol{\alpha}}}))}={{{\mathbf{F}^{{{\mathcal{I}}}}}}}_{\\#}[p({\mathbf{U}}\mid\mathbf{x}^{F})].$
(4)

Importantly, the posterior distribution $p({\mathbf{U}}\mid\mathbf{x}^{F})$
satisfies:

######  Property 3.

In the abduction step, statistical independence implies that conditioned on
the endogenous variables of the factual sample $\mathbf{x}^{F}$, each
exogenous variable $U_{i}$ is independent of the factual value $x^{F}_{j}$ if
$j\neq i$ and the variable $X_{j}$ is not a parent of $X_{i}$, i.e.,
$j{{\not\in}}{{{\mathrm{pa}^{*}}(i)}}$.

##  4 Variational Causal Autoencoder (VACA)

In this section, we present a novel variational causal graph autoencoder
(VACA) to approximate the observational (2), interventional (3) and
counterfactual (4) distributions. While the underlying SCM ${\mathcal{M}}$ is
unknown, we assume access to the true causal graph ${\mathcal{G}}$ and
observational data $\\{\mathbf{x}_{n}\\}^{N}_{n=1}$, i.e., i.i.d. samples of
the observational distribution induced by ${\mathcal{M}}$ (in the absence of
hidden confounders).

######  Definition 4.1.

(VACA) . Given a causal graph ${\mathcal{G}}$ over a set of endogenous
variables ${\mathbf{X}}=\\{X_{1},\ldots,X_{d}\\}$, which establishes the set
of parents ${{{\mathrm{pa}(i)}}}$ for each variable $X_{i}$ (including the
$i$-th node), VACA is defined by:

  * •

A causal adjacency matrix ${\mathbf{A}}$, which is a $d\times d$ binary matrix
with elements $A_{ij}=1$ if $j\in{{{\mathrm{pa}^{*}}(i)}}$, i.e., when $i=j$
or $j$ is a parent of $i$. Otherwise, ${{{{A}_{ij}}}}=0$.

  * •

A prior distribution $p({\mathbf{Z}})=\prod_{i}p(Z_{i})$ over the set of
latent variables ${\mathbf{Z}}=\\{Z_{1},\ldots,Z_{d}\\}$.

  * •

A decoder ${p_{\theta}}({\mathbf{X}}\mid{\mathbf{Z}},{\mathbf{A}})$, which is
a GNN (parameterized by $\theta$) that takes as input the set of latent
variables ${\mathbf{Z}}$ and the causal adjacency matrix ${\mathbf{A}}$, and
outputs the parameters of the likelihood
${p_{\theta}}({\mathbf{X}}\mid{\mathbf{Z}},{\mathbf{A}})$.

  * •

An encoder $q_{\phi}({\mathbf{Z}}\mid{\mathbf{X}},{\mathbf{A}})$, which is a
GNN (parameterized by $\phi)$ that takes as input the endogenous variables
${\mathbf{X}}$ and the causal adjacency matrix ${\mathbf{A}}$, and outputs the
parameters of the posterior approximation
$q_{\phi}({\mathbf{Z}}\mid{\mathbf{X}},{\mathbf{A}})$.

Next, we discuss how to design VACA such that it is able to capture the
observational, interventional, and counterfactual distribution induced by an
unknown SCM. Importantly, we derive the necessary conditions on the design of
both the encoder and decoder GNNs such that VACA can mimic the SCM properties
introduced in Section 3.

###  4.1 Observational distribution

VACA approximates the observational distribution in (2) using the generative
model as

$\displaystyle\begin{split}&p({\mathbf{X}})\approx\int{p_{\theta}}({\mathbf{X}}\mid{\mathbf{Z}},{\mathbf{A}}){p({\mathbf{Z}})}d{\mathbf{Z}},\end{split}$
(5)

where
${p_{\theta}}({\mathbf{X}}\mid{\mathbf{Z}},{\mathbf{A}})=\prod_{i=1}^{d}{p_{\theta}}(X_{i}\mid{\mathbf{Z}},{\mathbf{A}})$.
Figure 3(a) depicts this generative process.

Relationship between $\mathbf{Z}$ and $\mathbf{U}$. When comparing (5) with
the true observational distribution in (2), we observe that the latent
variables ${\mathbf{Z}}$ play a similar role to the exogenous variables
${\mathbf{U}}$, and the decoder
${p_{\theta}}({\mathbf{X}}\mid{\mathbf{Z}},{\mathbf{A}})$ plays a similar role
to the structural equations ${\mathbf{F}}$. We remark that ${\mathbf{Z}}$ do
not need to correspond to the true exogenous variables (i.e.,
$p({\mathbf{U}})\neq p({\mathbf{Z}})$), and thus, the decoder does not aim to
approximate the causal structural equations. Yet, we assume that there is one
independent latent variable $Z_{i}$ for every observed variable $X_{i}$
capturing all the information of $X_{i}$ that cannot be explained by its
parents. Thus, since $X_{i}$ is in turn a (deterministic) function of its
parents ${{{\mathrm{pa}(i)}}}$ and its exogenous variable $U_{i}$, the
posterior $p\left(Z_{i}\mid X_{i},{{{\mathrm{pa}(i)}}}\right)$ aims to capture
the information that $U_{i}$ contributes to $X_{i}$ (i.e., the information of
$X_{i}$ not contributed by its parents). That is–similar to $p\left(U_{i}\mid
X_{i},{{{\mathrm{pa}(i)}}}\right)$, the (true) posterior
distribution–$p\left(Z_{i}\mid X_{i},{{{\mathrm{pa}(i)}}}\right)$ should
depend only on $X_{i}$ and parents ${{{\mathrm{pa}(i)}}}$.

Observational noise. VACA has observational noise that is not present in the
true SCM, where an observed variable is assumed to be a deterministic
transformation of its exogenous variables and parents via the structural
equations (SEs). As VACA does not have access to the true SEs (nor to the true
distribution of the exogenous variables), the noise of the likelihood
$p_{\theta}({\mathbf{X}}\mid{\mathbf{Z}},{\mathbf{A}})$ can be interpreted as
an estimate of the uncertainty on the estimated observational distribution
(due to the uncertainty on the true SCM).

Here, we seek to ensure that the $p({\mathbf{X}})$ induced by VACA complies
with causal factorization (Property 1 in Section 3). To that end, the design
of the decoder GNN must assure that
${p_{\theta}}(X_{i}\mid{\mathbf{Z}},{\mathbf{A}})={p_{\theta}}(X_{i}\mid{\mathbf{Z}}_{{{{\mathrm{an}^{*}}(i)}}})$.
That is, that $X_{i}$ depends only on $Z_{j}$ if $j=i$ or $X_{j}$ is an
ancestor of $X_{i}$ in the causal graph.

######  Proposition 1.

(Causal factorization). VACA satisfies causal factorization,
${p_{\theta}}({\mathbf{X}}\mid{\mathbf{Z}},{\mathbf{A}})=\prod_{i}p_{\theta_{i}}(X_{i}\mid{\mathbf{Z}}_{{{{\mathrm{an}^{*}}(i)}}})$,
if and only if the number of hidden layers in the decoder is greater or equal
than ${\delta}-1$, with ${\delta}$ being the length of the longest shortest
path between any two endogenous nodes.

The above proposition (proved in Appendix B) is based on the fact that, in a
GNN with $N_{h}$ hidden layers (and $N_{h}+1$ layers in total), the output for
the $i$-th node depends on its neighbors of up to $N_{h}+1$ hops. As an
example, consider the following chain causal graph: $X_{1}\rightarrow
X_{2}\rightarrow X_{3}$, such that $\delta=2$. 4.1 In the decoder, the first
layer yields a hidden representation for the $3$-rd node
$h_{3}^{1}:=f(f(Z_{2}),Z_{3})$ that only depends on $Z_{2}$ and $Z_{3}$. Thus,
we need a second layer for its output
$h_{3}^{2}:=f(h_{2},Z_{3})=f(f(f(Z_{1}),Z_{2}),Z_{3})$ to depend on $Z_{1}$
(note that $X_{1}$ is an ancestor of $X_{3}$).

$X_{1}$$X_{2}$$X_{3}$$Z_{1}$$Z_{2}$$Z_{3}$•$h_{1}$•$h_{2}$•$h_{3}$

(a) Original

$X_{1}$$X_{2}$$X_{3}$$Z_{1}$$Z_{2}$$Z_{3}$•$h_{1}$•$h_{2}$•$h_{3}$

(b) Intervened

Figure 2: VACA decoder (a) with and (b) without intervening on $X_{2}$.
Message passing in the GNN correspond to direct (blue) and indirect (red)
causal paths in Figure 1.

(a) Observational

(b) Interventional

(c) Counterfactual

Figure 3: VACA generation of (a) observational, (b) interventional, and (c)
counterfactual samples. The ‘hat’ in $\hat{{\mathbf{X}}}$ and
$\hat{{\mathbf{x}}}^{CF}$ indicates that they are sample estimates of the true
random variables.

###  4.2 Interventional distribution

VACA approximates the interventional distribution in (3) as (illustrated in
Figure 2):

$\displaystyle p({\mathbf{X}}\mid do({{{\mathbf{X}_{{{\mathcal{I}}}}}}}$
$\displaystyle={{\boldsymbol{\alpha}}}))\approx\int\int{p_{\theta}}({\mathbf{X}}\mid{\tilde{{\mathbf{Z}}}},{\tilde{{\mathbf{Z}}}^{{{\mathcal{I}}}}},{{{\mathbf{A}^{{{\mathcal{I}}}}}}})p({\tilde{{\mathbf{Z}}}}){q_{\phi}}({\tilde{{\mathbf{Z}}}^{{{\mathcal{I}}}}}\mid{{{\mathbf{A}^{{{\mathcal{I}}}}}}},{{{\mathbf{X}_{{{\mathcal{I}}}}}}})d{\tilde{{\mathbf{Z}}}}d{\tilde{{\mathbf{Z}}}^{{{\mathcal{I}}}}},$
(6)

where
${\tilde{{\mathbf{Z}}}^{{{\mathcal{I}}}}}=\\{{Z}^{{{\mathcal{I}}}}_{i}\\}_{i\in{{\mathcal{I}}}}$
is the subset of latent variables associated with the intervened-upon
variables ${{{\mathbf{X}_{{{\mathcal{I}}}}}}}$, and
${\tilde{{\mathbf{Z}}}}=\\{Z_{i}\\}_{i{{\not\in}}{{\mathcal{I}}}}$ denotes the
subset of latent variables associated with the rest of the observed variables.
Importantly, here the do-operator is performed on the causal adjacency matrix
as
$do(X_{{\mathcal{I}}}={{\boldsymbol{\alpha}}}):{\mathbf{A}}\mapsto{{{\mathbf{A}^{{{\mathcal{I}}}}}}}=\\{{{{{A}_{ij}}}}\\}_{{{{{\text{
}{\forall}}}}}i\not\in{{\mathcal{I}}},j}\cup\\{{{{{A}_{ij}}}}=0\\}_{{{{{\text{
}{\forall}}}}}i\in{{\mathcal{I}}},j}$. This ensures that $X_{i}$ for
$i\in{{\mathcal{I}}}$ is independent of $Z_{j}$ for all $j\neq i$. Note that
in order for (6) to be able to approximate the interventional distribution in
(3), an intervention on VACA should satisfy Property 2, i.e.:

######  Proposition 2.

(Causal interventions). VACA captures causal interventions if and only if the
number of hidden layers in its decoder is greater than or equal to $\gamma-1$,
with $\gamma$ being the length of the longest path between any two endogenous
nodes in ${\mathcal{G}}$.

To illustrate this, Figure 2 depicts how messages are exchanged in a one-
hidden-layer decoder GNN corresponding to the causal graph ${\mathcal{G}}$ in
Figure 1 (_triangle_ with ${\gamma}=2$), both (a) without and (b) with an
intervention on $X_{2}$. We highlight in blue the direct messages (sent via
direct causal path in ${\mathcal{G}}$), and in red the indirect messages (sent
via indirect causal path in ${\mathcal{G}}$) from $Z_{1}$ to $X_{3}$. Observe
that, similarly to Figure 1, in (a) there is an indirect path (via $h_{2}$)
from $Z_{1}$ to $X_{3}$; while in (b) this path is severed. Hence, the hidden
layer $(h_{1},h_{2},h_{3})$ allows to distinguish between direct and indirect
paths and thus to capture interventional effects. As the condition in
Proposition 2 is more restrictive than the one in Proposition 1, VACA is able
to approximate the observational and interventional distributions (as
empirically validated in Appendix D) if:

Design condition 1 (necessary condition) _The decoder GNN of VACA has at least
as many hidden layers as ${\gamma}-1$, with ${\gamma}$ being the longest
directed path in the causal graph ${\mathcal{G}}$._

###  4.3 Counterfactual distribution

VACA approximates the counterfactual distribution in (4) as (illustrated in
Figure 3(c)):

$\displaystyle\begin{split}p(\mathbf{x}^{CF}&\mid
do({{{\mathbf{X}_{{{\mathcal{I}}}}}}}={{\boldsymbol{\alpha}}}),\mathbf{x}^{F})\approx\underbrace{\int\int\underbrace{{p_{\theta}}({\mathbf{X}}\mid{\tilde{{\mathbf{Z}}}^{{\mathbf{F}}}},{\tilde{{\mathbf{Z}}}^{{{\mathcal{I}}}}},{{{\mathbf{A}^{{{\mathcal{I}}}}}}}){q_{\phi}}({\tilde{{\mathbf{Z}}}^{{{\mathcal{I}}}}}\mid{{\mathbf{x^{{{\mathcal{I}}}}}}},{{{\mathbf{A}^{{{\mathcal{I}}}}}}})}_{action}\underbrace{{q_{\phi}}({\tilde{{\mathbf{Z}}}^{{\mathbf{F}}}}\mid\mathbf{x}^{F},{\mathbf{A}})}_{abduction}d{\tilde{{\mathbf{Z}}}^{{{\mathcal{I}}}}}d{\tilde{{\mathbf{Z}}}^{{\mathbf{F}}}}}_{prediction},\end{split}$
(7)

where $\mathbf{x}^{F}$ represents a sample from ${\mathbf{X}}$ for which we
seek to compute the distribution over counterfactual $\mathbf{x}^{CF}$ and
${\tilde{{\mathbf{Z}}}^{{\mathbf{F}}}}=\\{{Z}^{F}_{i}\\}_{i{{\not\in}}{{\mathcal{I}}}}$.
Note that two different passes of the encoder are necessary: one for the
_abduction_ step of the factual instance
${q_{\phi}}({\tilde{{\mathbf{Z}}}^{{\mathbf{F}}}}\mid\mathbf{x}^{F},{\mathbf{A}})$;
and another one for the _action_ step (intervention)
${q_{\phi}}({\tilde{{\mathbf{Z}}}^{{{\mathcal{I}}}}}\mid{{\mathbf{x^{{{\mathcal{I}}}}}}},{{{\mathbf{A}^{{{\mathcal{I}}}}}}})$
with $x^{{\mathcal{I}}}_{i}=\alpha_{i}{{{{\text{
}{\forall}}}}}i\in{{\mathcal{I}}}$ (we remark that the rest of the values in
${{\mathbf{x^{{{\mathcal{I}}}}}}}$ do not affect the overall counterfactual
computation).

We then evaluate the likelihood making sure that the resulting counterfactual
sample $\mathbf{x}^{CF}$ only depends on the
${\tilde{{\mathbf{Z}}}^{{\mathbf{F}}}}$ and
${\tilde{{\mathbf{Z}}}^{{{\mathcal{I}}}}}$. Importantly, in order for VACA to
be able to approximate the counterfactual distribution, we need its abduction
(and action) step(s) to comply with Property 3, i.e.:

######  Proposition 3.

(Abduction). The abduction step of an observed sample
${\mathbf{x}}=\\{x_{1},\ldots,x_{d}\\}$ in VACA satisfies that for all $i$ the
posterior of $Z_{i}$ is independent on the subset
$\\{x_{j}\\}_{j{{\not\in}}{{{\mathrm{pa}^{*}}(i)}}}\subseteq{\mathbf{x}}$, if
and only if the encoder GNN has no hidden layers.

The above result (proved in Appendix B) can be shown by the message passing
algorithm computed by the encoder GNN, and leads to the second design
condition of VACA:

Design condition 2 (necessary condition): _The encoder GNN of VACA has no
hidden layers._

In other words, from the definition of the SCM, the posterior distribution of
$U_{i}$ only depends on the parents, i.e $U_{i}\mid X_{{{{\mathrm{pa}(i)}}}}$.
In order for VACA to mimic this property, the GNN that parameterized the
encoder contains no hidden layers: in the message passing algorithm, in the
$k$-th iteration (layer), a node depends on its $k$-hop ancestors; requiring
$k=1$ in the encoder refers to a GNN without hidden layers. Note that while
the above condition may look restrictive and limiting the capacity of our
encoder, we may choose arbitrarily complex NNs for the message $f^{m}$ and
update $f^{u}$ functions, as well as one or more aggregation functions
$f^{a}$, e.g., sum or max, to model the encoder [4].

###  4.4 Practical considerations

Next, we briefly discuss practical implementation considerations to handle
complex causal models, which often appear in real world applications [7, 6]
For further details on VACA implementation, refer to Appendix C.

Heterogeneous causal nodes. So far, we have modeled each endogenous variable
$X_{i}$ as a node in the causal graph ${\mathcal{G}}$, and thus in the VACA
GNNs. Yet, in some application domains the relationships between a subset of
$k_{i}$ variables may be unknown, or they may be affected by hidden
confounders. In such cases, we assume that set of $k_{i}$ variables to be
correlated and model them as one multidimensional and potentially
heterogeneous node ${\mathbf{X}}_{i}=\\{X_{i1},\dots,X_{ik_{i}}\\}$ that share
the same latent random variable $Z_{i}$. This allows us to deal with a large
variety of graphs in practice.

Heterogeneous endogenous variables. Heterogeneous causal nodes require us to
model different functions for each node, i.e. nodes may now contain a mix of
continuous/discrete variables. In general GNNs are parametrized such that the
parameters of the message function $f^{m}$ and update function $f^{u}$ are
shared for all the nodes and edges in the graph. However, similar to the
structural equations ${\mathbf{F}}$, we can define a unique set of parameters
$\theta_{mij}$ for each $f^{m}_{ij}$ (see (1)), so that we can model a
different function for every edge in the causal graph. Further, we can also
assume different update functions $f^{u}_{i}$ for each node $i$, by
introducing different update parameters $\theta_{ui}$. Note, however, that
VACA can be used with any GNN architecture [4, 15].

##  5 Evaluation

In this section, we evaluate the potential of VACA in approximating the
outcomes of causal queries and compare it to two competing methods in
synthetic experiments. The synthetic setting allows us to have access the true
SEs, which is necessary to evaluate interventional distributions and
especially counterfactuals. We consider interventions of the form
$do(x_{i}=\alpha_{i})$ for several values of $\alpha_{i}$ on both root and
non-root nodes. We compute all results over the same 10 random seeds and
report mean and standard deviation. Refer to Appendix E for a complete
description of the experimental setup. Moreover, our code is publicly
available at GitHub111https://github.com/psanch21/VACA.

Datasets. We consider 6 different synthetic causal graphs that differ in the
number of nodes $d$, diameter ${\delta}$, and longest path ${\gamma}$. Here,
we report the results for i) the collider ($d=3$, ${\delta}=1$, ${\gamma}=1$)
with linear (LIN) and non-linear (NLIN) additive noise SEs, ii) the loan from
Karimi et al. [16] ($d=7$, ${\delta}=2$, ${\gamma}=3$), and iii) the adult
($d=11$, ${\delta}=2$, ${\gamma}=3$) graphs. Note that the two latter ones are
synthetic versions of the German Credit dataset [7] and the Adult datasets
[6], respectively. See Appendix E.1 for further details on the graphs and
Appendix E.4 for the results with the remaining graphs.

Metrics. We evaluate the observational distribution using the Maximum Mean
Discrepancy (MMD) [13] as distance-measure between the true and estimated
distributions, i.e., the lower the MMD the better the distributions match. For
the interventional distribution, we additionally report the average estimation
squared error of the mean (MeanE) and of the standard deviation (StdE) over
all descendants of the intervened-upon variables. For the counterfactual
distribution we report the mean squared error (MSE) as well as the standard
deviation of the squared error (SSE) between the true and the estimated
counterfactual value. More details in Appendix E.3.

Baselines. We compare VACA with MultiCVAE [16] and CAREFL [17] described in
Section 1.1. For a fair comparison, all model hyperparameters have been cross-
validated using a similar computational budget (see Appendix E.2). In Table
LABEL:tab:results, we report for each model and SCM the best configuration
according to observational MMD. We also include a time-complexity analysis in
Appendix E.5.

Results. Table LABEL:tab:results summarizes the results. We observe that, in
general, MultiCVAE underperforms the other methods. This may be explained by
the fact that MultiCVAE trains each node independently, and thus the
discrepancy between the true and generated distributions in one node may be
amplified in its descendants. Comparing VACA to CAREFL, we first observe that
VACA performs consistently better in terms of observational MMD, i.e., VACA is
able to generate observational samples that better resemble the true ones.
Second, regarding the interventional distribution, CAREFL does a good job at
fitting the mean (i.e., low MeanE). However,VACA performs consistently better
both at approximating the standard deviation (i.e., low StdE) and the true
samples (i.e., low MMD). This can be explained by VACA i) leveraging the
causal graph (contrary to CAREFL that relies on causal ordering), and ii)
optimizing the log-evidence in (5) jointly (contrary to the sequential
optimization of MultiCVAE). Thus, VACA approximates the distribution as a
whole better, which is a desirable property for studying interventions on a
population-level rather than just on average. Lastly, in the approximation of
counterfactuals we observe that both CAREFL and VACA exhibit similar
performance in terms of MSE and SSE. Note however that CAREFL performs exact
inference while VACA is built on approximate inference and is trained on a
lower bound on the log-evidence. Finding tighter bounds could boost VACA
performance.

##  6 Use case: counterfactual fairness

We finally show two practical use-cases of our method: assessing
counterfactual fairness and training counterfactually fair classifiers. We use
the public German Credit dataset [7] and rely on the causal model proposed by
Chiappa [3] with the following random variables ${\mathbf{X}}$: sensitive
feature $S=\\{\textit{sex}\\}$, and non-sensitive features
$C=~{}\\{\textit{age}\\}$, $R=\\{\text{{credit amount}, {repayment
history}\\}}$ and $H=\\{\textit{checking
account},\textit{savings},\textit{housing}\\}$. Then, we aim to predict the
binary feature $Y=\\{\textit{credit risk}\\}$ from ${\mathbf{X}}$. See
Appendix F for further details.

Counterfactual fairness. Let $S\subset{\mathbf{X}}$ be a sensitive attribute
(e.g., gender), then the counterfactual unfairness [23] of a classifier
${{h}}:{\mathbf{X}}\rightarrow Y$ is measured
$\forall\mathbf{x}^{CF},\alpha^{\prime}\neq\alpha,y$ as:

$\displaystyle{\textit{uf}}=\mid P({{h}}(\mathbf{x}^{CF})=y\mid
do(S=\alpha),\mathbf{x}^{F})-P({{h}}(\mathbf{x}^{CF})=y\mid
do(S=\alpha^{\prime}),\mathbf{x}^{F})\mid$ (8)

A classifier is counterfactually fair (${\textit{uf}}=0$), if, given a factual
$\mathbf{x}^{F}$ with sensitive attribute $S=\alpha$, had its sensitive
attribute been different $S=\alpha^{\prime}$, the classifier prediction would
remain the same. We can use VACA to generate counterfactual estimates to
_audit_ the fairness level of a classifier. Following [23], we audit: i) a
full model ${{h}_{\text{full}}}:{\mathbf{X}}\rightarrow Y$ that takes as input
the complete variable set; ii) an unaware model
${{h}_{\text{unaw}}}:{\mathbf{X}}\backslash S\rightarrow Y$ that takes as
input all variables but the sensitive one; iii) and a fair model
${{h}_{\text{fair-x}}}:\\{X_{i}|S{{\not\in}}{{{\mathrm{an}^{*}}(i)}}\\}\rightarrow
Y$ that takes as input all non-descendant variables of the sensitive
attribute. Moreover, we show that we can _learn a fair classifier_
${{h}_{\text{fair-z}}}:{\mathbf{Z}}\backslash Z_{S}\rightarrow Y$, which takes
as input the latent variables generated by the VACA encoder without the one of
the sensitive attribute $Z_{S}$.

Fairness Auditing. Table LABEL:tab:fairness summarizes the unfairness level
and f1-score for a support vector machine (SVM) classifier. See Appendix F for
results of a logistic regression classifier. As we do not have access to the
true data generation process, we evaluate the auditing task by the resulting
ranking of the different classifiers according to their unfairness level.
Based on the counterfactual generation by VACA the _full_ classifier is
consistently less fair than the _unaware_ and the _fair-x_ classifier,
respectively. This ranking is consistent with the one in [23].

Fairness Classification Table LABEL:tab:fairness shows that for the fair-x
classifier fairness comes at the expense of accuracy compared to the full
classifier. On the contrary, even though VACA has been trained for
representation learning without access to classification labels, fair-z is a
fair classifier (with comparable fairness level to the fair-x one) while
keeping the performance comparable to the unfair full classifier. VACA,
therefore, does not only allow us to audit counterfactual fairness but also
provides a practical approach to train accurate and fair classifiers.

##  7 Conclusion, Limitations and Impact

In this work, we have proposed VACA, a variational causal autoencoder based on
GNNs that: i) is specially designed to capture the properties of SCMs; ii)
inherently handles heterogeneous data; and iii) provides good approximations
of interventional and counterfactual distributions as a whole for SCMs of
different complexities. As demonstrated by extensive synthetic experiments,
VACA provides accurate results for a wide range of interventions in diverse
SCMs leading to more consistent results than competing methods [16, 17].
Finally, we have applied VACA for counterfactually fair classification.

Practical limitations. The expressive power of VACA to model complex
structural equations, e.g., in domains such as biology [37], is limited by the
GNN architectures of the encoder and the decoder. As discussed in the GNN
literature [4], especially aggregation functions may limit expressiveness. We
expect VACA to benefit from advances in the field. Second, long causal paths
would require VACA to increase the number of layers in the decoder (see Design
condition 1). However, the GNNs performance is known to deteriorate with depth
[9, 14, 24].

Social impact. Trusting counterfactuals is of great importance for decision
making, e.g. in the political or medical domain. We thus encourage anyone who
uses VACA (or any other ML method for causal inference) to i) fully understand
the model assumptions and to verify (up to the possible extend) that they are
fulfilled; as well as ii) to be aware of the identifiability problem in
counterfactual queries [49, 45].

Future work. First, it would be important to evaluate the sensitivity of VACA
to errors in the assumed causal graph, as well as to the presence of hidden
confounders. We plan to extend VACA to handle more complex causal models
including, e.g., hidden confounders and non-DAG causal graphs. Second, it
would be interesting to perform ablation studies on the limitations of
available GNNs architectures [44] for the VACA encoder and decoder; as well as
on how the performance deteriorates as we increase the length of the causal
path and thus the required number of hidden layers [24]. Finally, it would be
intriguing to apply VACA to other causal questions recently discussed in the
machine learning literature, such as privacy-preserving causal inference [22]
or explainable machine learning [16].

## Acknowledgments

We would like to thank Adrián Javaloy Bornás, Amir Hossein-Karimi, Jonas
Kleesen and Maryam Meghdadi Esfahani for helpful feedback and discussions.
Moreover, a special thanks to Diego Baptista Theuerkauf for invaluable help
with formalizing proofs. Moreover, the authors would like to thank Ilyes
Khemakhem for helpful insights in how to generalize their CAREFL approach to
arbitrary graphs.

Pablo Sánchez Martín is supported by the German Research Foundation through
the Cluster of Excellence “Machine Learning – New Perspectives for Science”,
EXC 2064/1, project number 390727645. Miriam Rateike is supported by the
German Federal Ministry of Education and Research (BMBF): Tübingen AI Center,
FKZ: 01IS18039B. The authors thank for the generous funding support. The
authors also thank the International Max Planck Research School for
Intelligent Systems (IMPRS-IS) for supporting Pablo Sańchez Martín.

## References

  * Bastings et al. [2017] Jasmijn Bastings, Ivan Titov, Wilker Aziz, Diego Marcheggiani, and Khalil Sima’an.  Graph convolutional encoders for syntax-aware neural machine translation.  In _Proceedings of the Conference on Empirical Methods in Natural Language Processing (EMNLP)_ , volume 3, 2017.
  * Burda et al. [2016] Yuri Burda, Roger Grosse, and Ruslan Salakhutdinov.  Importance weighted autoencoders.  In _Proceedings of the International Conference onLearning Representations (ICLR)_ , volume 4, 2016.
  * Chiappa [2019] Silvia Chiappa.  Path-specific counterfactual fairness.  In _Proceedings of the AAAI Conference on Artificial Intelligence_ , volume 33, 2019.
  * Corso et al. [2020] Gabriele Corso, Luca Cavalleri, Dominique Beaini, Pietro Liò, and Petar Veličković.  Principal neighbourhood aggregation for graph nets.  In _Advances in Neural Information Processing Systems (NeurIPS)_ , volume 33, 2020.
  * Derrow-Pinion et al. [2021] Austin Derrow-Pinion, Jennifer She, David Wong, Oliver Lange, Todd Hester, Luis Perez, Marc Nunkesser, Seongjae Lee, Xueying Guo, Brett Wiltshire, et al.  Eta prediction with graph neural networks in google maps.  _arXiv preprint arXiv:2108.11482_ , 2021.
  * Dua and Graff [2017a] Dheeru Dua and Casey Graff.  UCI machine learning repository, 2017a.  URL https://archive.ics.uci.edu/ml/datasets/adult.
  * Dua and Graff [2017b] Dheeru Dua and Casey Graff.  UCI machine learning repository, 2017b.  URL https://archive.ics.uci.edu/ml/datasets/statlog+(german+credit+data).
  * Durkan et al. [2019] Conor Durkan, Artur Bekasov, Iain Murray, and George Papamakarios.  Neural spline flows.  In _Advances in Neural Information Processing Systems_ , volume 32, 2019.
  * Gallicchio and Micheli [2020] Claudio Gallicchio and Alessio Micheli.  Fast and deep graph neural networks.  In _Proceedings of the AAAI Conference on Artificial Intelligence_ , volume 34, 2020.
  * Garrido et al. [2021] Sergio Garrido, Stanislav Borysov, Jeppe Rich, and Francisco Pereira.  Estimating causal effects with the neural autoregressive density estimator.  _Journal of Causal Inference_ , 9, 2021.
  * Gilmer et al. [2017] Justin Gilmer, Samuel S Schoenholz, Patrick F Riley, Oriol Vinyals, and George E Dahl.  Neural message passing for quantum chemistry.  In _Proceedings of the International Conference on Machine Learning (ICML)_ , volume 34. PMLR, 2017.
  * Glymour et al. [2019] Clark Glymour, Kun Zhang, and Peter Spirtes.  Review of causal discovery methods based on graphical models.  _Frontiers in genetics_ , 10:524, 2019.
  * Gretton et al. [2012] Arthur Gretton, Karsten M Borgwardt, Malte J Rasch, Bernhard Schölkopf, and Alexander Smola.  A kernel two-sample test.  _The Journal of Machine Learning Research (JMLR)_ , 13, 2012.
  * Gu et al. [2020] Fangda Gu, Heng Chang, Wenwu Zhu, Somayeh Sojoudi, and Laurent El Ghaoui.  Implicit graph neural networks.  In _Advances in Neural Information Processing Systems_ , volume 33, 2020.
  * Hamilton et al. [2017] Will Hamilton, Zhitao Ying, and Jure Leskovec.  Inductive representation learning on large graphs.  In _Advances in Neural Information Processing Systems (NeurIPS)_ , volume 30, 2017.
  * Karimi et al. [2020] Amir-Hossein Karimi, Julius von Kügelgen, Bernhard Schölkopf, and Isabel Valera.  Algorithmic recourse under imperfect causal knowledge: a probabilistic approach.  In _Advances in Neural Information Processing Systems_ , volume 33, pages 265–277, 2020.
  * Khemakhem et al. [2021] Ilyes Khemakhem, Ricardo Monti, Robert Leech, and Aapo Hyvarinen.  Causal autoregressive flows.  In _Proceedings of International Conference on Artificial Intelligence and Statistics (AISTATS)_ , volume 24. PMLR, 2021.
  * Kim et al. [2021] Hyemi Kim, Seungjae Shin, JoonHo Jang, Kyungwoo Song, Weonyoung Joo, Wanmo Kang, and Il-Chul Moon.  Counterfactual fairness with disentangled causal effect variational autoencoder.  In _Proceedings of the AAAI Conference on Artificial Intelligence_ , volume 35, 2021.
  * Kingma and Welling [2014] Diederik P Kingma and Max Welling.  Auto-encoding variational bayes.  In _Proceedings of the International Conference on Learning Representations (ICLR)_ , volume 2, 2014.
  * Kipf and Welling [2016] Thomas N Kipf and Max Welling.  Variational graph auto-encoders.  _arXiv preprint arXiv:1611.07308_ , 2016.
  * Kocaoglu et al. [2018] Murat Kocaoglu, Christopher Snyder, Alexandros G. Dimakis, and Sriram Vishwanath.  Causalgan: Learning causal implicit generative models with adversarial training.  In _Proceedings of the International Conference on Learning Representations (ICLR)_ , volume 6, 2018.
  * Kusner et al. [2016] Matt J Kusner, Yu Sun, Karthik Sridharan, and Kilian Q Weinberger.  Private causal inference.  In _Proceedings of the Conference on Artificial Intelligence and Statistics (AISTATS)_ , volume 19. PMLR, 2016.
  * Kusner et al. [2017] Matt J Kusner, Joshua Loftus, Chris Russell, and Ricardo Silva.  Counterfactual fairness.  In _Advances in Neural Information Processing Systems (NeurIPS)_ , volume 30, 2017.
  * Li et al. [2018] Qimai Li, Zhichao Han, and Xiao-Ming Wu.  Deeper insights into graph convolutional networks for semi-supervised learning.  In _Proceedings of the AAAI Conference on Artificial Intelligence_ , volume 32, 2018.
  * Louizos et al. [2017] Christos Louizos, Uri Shalit, Joris M Mooij, David Sontag, Richard Zemel, and Max Welling.  Causal effect inference with deep latent-variable models.  In _Advances in Neural Information Processing Systems (NeurIPS)_ , volume 30, 2017.
  * Moraffah et al. [2020] Raha Moraffah, Bahman Moraffah, Mansooreh Karami, Adrienne Raglin, and Huan Liu.  Can: A causal adversarial network for learning observational and interventional distributions.  _arXiv preprint arXiv:2008.11376_ , 2020.
  * Nowozin [2018] Sebastian Nowozin.  Debiasing evidence approximations: On importance-weighted autoencoders and jackknife variational inference.  In _Proceedings of the International Conference on Learning Representations (ICML)_ , volume 35. PMLR, 2018.
  * Parafita and Vitrià [2019] Álvaro Parafita and Jordi Vitrià.  Explaining visual models by causal attribution.  In _2019 IEEE/CVF International Conference on Computer Vision Workshop (ICCVW)_ , pages 4167–4175. IEEE, 2019.
  * Parafita and Vitrià [2020] Álvaro Parafita and Jordi Vitrià.  Causal inference with deep causal graphs.  _arXiv preprint arXiv:2006.08380_ , 2020.
  * Parafita and Vitrià [2019] Álvaro Parafita and Jordi Vitrià.  Explaining visual models by causal attribution.  In _International Conference on Computer Vision Workshop (ICCVW)_ , 2019.
  * Pawlowski et al. [2020] Nick Pawlowski, Daniel Coelho de Castro, and Ben Glocker.  Deep structural causal models for tractable counterfactual inference.  In _Advances in Neural Information Processing Systems (NeurIPS)_ , volume 33, 2020.
  * Pearl [2009a] Judea Pearl.  Causal inference in statistics: An overview.  _Statistics surveys_ , 3, 2009a.
  * Pearl [2009b] Judea Pearl.  _Causality_.  Cambridge university press, 2009b.
  * Peters et al. [2017] Jonas Peters, Dominik Janzing, and Bernhard Schölkopf.  _Elements of causal inference: Foundations and learning algorithms_.  The MIT Press, 2017.
  * Rainforth et al. [2018] Tom Rainforth, Adam Kosiorek, Tuan Anh Le, Chris Maddison, Maximilian Igl, Frank Wood, and Yee Whye Teh.  Tighter variational bounds are not necessarily better.  In _Proceedings of the International Conference on Machine Learning (ICML)_ , volume 35. PMLR, 2018.
  * Rakesh et al. [2018] Vineeth Rakesh, Ruocheng Guo, Raha Moraffah, Nitin Agarwal, and Huan Liu.  Linked causal variational autoencoder for inferring paired spillover effects.  In _Proceedings of the International Conference on Information and Knowledge Management (CIKM)_. ACM, 2018.
  * Sachs et al. [2005] Karen Sachs, Omar Perez, Dana Pe’er, Douglas A Lauffenburger, and Garry P Nolan.  Causal protein-signaling networks derived from multiparameter single-cell data.  _Science_ , 308(5721):523–529, 2005.
  * Schölkopf [2019] Bernhard Schölkopf.  Causality for machine learning.  _arXiv preprint arXiv:1911.10500_ , 2019.
  * Schwab et al. [2018] Patrick Schwab, Lorenz Linhardt, and Walter Karlen.  Perfect match: A simple method for learning representations for counterfactual inference with neural networks.  _arXiv preprint arXiv:1810.00656_ , 2018.
  * Shen et al. [2020] Xinwei Shen, Furui Liu, Hanze Dong, Qing Lian, Zhitang Chen, and Tong Zhang.  Disentangled generative causal representation learning.  _arXiv preprint arXiv:2010.02637_ , 2020.
  * Tucker et al. [2018] George Tucker, Dieterich Lawson, Shixiang Gu, and Chris J Maddison.  Doubly reparameterized gradient estimators for monte carlo objectives.  In _International Conference on Learning Representations_ , 2018.
  * Vowels et al. [2021] Matthew J Vowels, Necati Cihan Camgoz, and Richard Bowden.  D’ya like dags? a survey on structure learning and causal discovery.  _arXiv preprint arXiv:2103.02582_ , 2021.
  * Vowels et al. [2020] Matthew James Vowels, Necati Cihan Camgoz, and Richard Bowden.  Targeted vae: Structured inference and targeted learning for causal parameter estimation.  _arXiv preprint arXiv:2009.13472_ , 2020.
  * Wu et al. [2020] Zonghan Wu, Shirui Pan, Fengwen Chen, Guodong Long, Chengqi Zhang, and S Yu Philip.  A comprehensive survey on graph neural networks.  _IEEE Transactions on Neural Networks and Learning Systems_ , 2020\.
  * Xia et al. [2021] Kevin Xia, Kai-Zhan Lee, Yoshua Bengio, and Elias Bareinboim.  The causal-neural connection: Expressiveness, learnability, and inference.  _arXiv preprint arXiv:2107.00793_ , 2021.
  * Yang et al. [2020] Mengyue Yang, Furui Liu, Zhitang Chen, Xinwei Shen, Jianye Hao, and Jun Wang.  Causalvae: Disentangled representation learning via neural structural causal models.  _arXiv preprint arXiv:2004.08697_ , 2020.
  * Yu et al. [2018] Bing Yu, Haoteng Yin, and Zhanxing Zhu.  Spatio-temporal graph convolutional networks: a deep learning framework for traffic forecasting.  In _Proceedings of the International Joint Conference on Artificial Intelligence (IJCAI)_ , volume 27, 2018.
  * Yu et al. [2019] Yue Yu, Jie Chen, Tian Gao, and Mo Yu.  Dag-gnn: Dag structure learning with graph neural networks.  In _Proceedings of the International Conference on Machine Learning (ICML)_ , volume 36. PMLR, 2019.
  * Zečević et al. [2021] Matej Zečević, Devendra Singh Dhami, Petar Veličković, and Kristian Kersting.  Relating graph neural networks to structural causal models.  _arXiv preprint arXiv:2109.04173_ , 2021.
  * Zhang et al. [2020] Cheng Zhang, Kun Zhang, and Yingzhen Li.  A causal view on robustness of neural networks.  _Advances in Neural Information Processing Systems_ , 33, 2020.
  * Zhang et al. [2019] Muhan Zhang, Shali Jiang, Zhicheng Cui, Roman Garnett, and Yixin Chen.  D-vae: A variational autoencoder for directed acyclic graphs.  In _Advances in Neural Information Processing Systems (NeurIPS)_ , volume 32, 2019.
  * Zheng and Kleinberg [2019] Min Zheng and Samantha Kleinberg.  Using domain knowledge to overcome latent variables in causal inference from time series.  In _Proceedings of the Machine Learning for Healthcare Conference (MLHC)_. PMLR, 2019.

##  Appendix A Background details on message passing Graph Neural Networks

A directed graph with $|{\mathcal{V}}|=d$ nodes can be represented as
${\mathcal{G}}:=({\mathcal{V}},{\mathcal{E}})$, where ${\mathcal{V}}$ denotes
the set of nodes and $(j,i)\in{\mathcal{E}}$ is the set of directed edges from
$j$ to $i$. Additionally, we denote with ${\mathbf{X}}\in{\mathbb{R}}^{d\times
F_{X}}$ the features of the nodes (the row index identifies the node, i.e.,
the $i$-th row contains the $F_{X}$-dimensional features of the node $i$)
Also, the adjacency matrix ${\mathbf{A}}\in\\{0,1\\}^{d\times d}$ of
${\mathcal{G}}$ is defined as $A_{ij}=1$ if there is an edge from $j$ to $i$
and $A_{ij}=0$ otherwise. Then, a directed graph can be alternatively
represented as ${\mathcal{G}}=({\mathbf{X}},{\mathbf{A}})$. Given a graph
${\mathcal{G}}$, a Graph Neural Network (GNN) with parameters $\theta$ is a
function $f_{\theta}:{\mathbb{R}}^{d\times F_{X}}\times\\{0,1\\}^{d\times
d}\rightarrow{\mathbb{R}}^{d\times F_{H}}$ that takes into account the graph
structure contained in the adjacency matrix ${\mathbf{A}}\in\\{0,1\\}^{d\times
d}$ and transforms the node features ${\mathbf{X}}$ into different features
$\mathbf{H}\in{\mathbb{R}}^{d\times F_{H}}$, i.e.,
$\mathbf{H}=f({\mathbf{X}},{\mathbf{A}};\theta)$ (for readability we consider
$f_{\theta}(\cdot)\equiv f(\cdot;\theta)$). Importantly, at the output of the
GNN we have a graph $(\mathbf{H},{\mathbf{A}})$ that preserves the structure
of the input graph $({\mathbf{X}},{\mathbf{A}})$.

A GNN based on message passing [11] is a type of spatial convolution GNNs [44]
in which information is passed following the message passing process: In each
layer $l$ of the GNN each node $i$ receives information from its neighbors
${\mathcal{N}}_{i}$, a.k.a. the parents of $i$. In a message passing GNN, the
feature vector of the $i$-th node at the output of a layer $l$—i.e.,
$h_{i}^{l}$—is computed in three steps:

  1. 1.

Message. The message from node $j$ to node $i$ is defined as
$m_{ij}^{l}=f^{m}(h_{i}^{l-1},h_{j}^{l-1};\theta^{l}_{m})$, where
$h_{i}^{l-1}$ are the features of node $i$ at layer $l-1$, $h_{j}^{l-1}$ are
the features of node $j$ at the previous layer $l-1$, and $f^{m}$ is a neural
network (usually a linear layer) parametrized by $\theta^{l}_{m}$.

  2. 2.

Aggregator. The aggregator is a function in charge of combining all the
incoming messages at each node $i$ into a single message, a.k.a. the
aggregated message $M_{i}^{l}=f^{a}(\\{m_{ij}^{l}\mid
j\in{\mathcal{N}}_{i}\\})$. Notice that $f^{a}$ does not have any parameters.
Choices of $f^{a}$ may be the sum, mean, standard deviation, max or min over
the inputs, i.e., messages [4].

  3. 3.

Update. The update function
$h_{i}^{l}=f^{u}(h_{i}^{l-1},M_{i}^{l};\theta_{u}^{l})$ takes the aggregated
message and the representation of node $i$ at layer $l-1$ and outputs the new
representation for node $i$ at layer $l$. The function $f^{u}$ is defined as a
neural network (usually a linear layer) with parameters $\theta^{l}_{u}$.

Putting the three steps together, we obtain the general form of a message
passing based GNN layer as
$h_{i}^{l}=f^{u}\left(h_{i}^{l-1},f^{a}\left(\\{f^{m}(h_{i}^{l-1},h_{j}^{l-1};\theta^{l}_{m})\mid
j\in{{\mathcal{N}_{i}}}\\}\right);\theta^{l}_{u}\right)$.

Algorithm 1 describes the propagation of information (i.e., messages ) in a
GNN with $L$ layers.

Input:  A directed graph ${\mathcal{G}}$ with $d$ nodes, adjacency matrix
${\mathbf{A}}$ and node features ${\mathbf{X}}$.

Output:  $\mathbf{H}=\\{h_{i}^{L}\\}_{i=1}^{d}$.

$h^{0}_{i}=x_{i}$ $\forall i$

for  _$l=1,\dots,L$ ;_

__ // For each layer $l$

do

_ _ for  _$i=1,\dots,d$ ;_

__ // For node $i$

do

$m_{ij}^{l}=f^{m}(h_{i}^{l-1},h_{j}^{l-1};\theta^{l}_{m})$ $\forall
j\in{\mathcal{N}}_{i}$;

// Compute the messages

$M_{i}^{l}=f^{a}(\\{m_{ij}^{l}\\}_{j\in{\mathcal{N}}_{i}})$ ;

// Compute the aggregated message

$h_{i}^{l}=f^{u}(h_{i}^{l-1},M_{i}^{l};\theta_{u}^{l})$ ;

// Compute the node features at layer $l$

end for

end for

Algorithm 1 Message passing GNN with $L$ layers decomposed in three operations

##  Appendix B Proofs

For completeness, this section first formalizes the meaning of causal
factorization, interventions and the abduction step in VACA.

VACA causal factorization refers to the factorization of the joint
distribution as

$\displaystyle{p_{\theta}}({\mathbf{X}}\mid{\mathbf{Z}},{\mathbf{A}})=\prod_{i}p_{\theta_{i}}(X_{i};{\boldsymbol{\eta}}_{i}),$

where the likelihood parameters
${\boldsymbol{\eta}}_{i}={\boldsymbol{\eta}}_{i}({\mathbf{Z}}_{{{{\mathrm{an}^{*}}(i)}}})$
are a function of all (and only) the features of the ancestors of $i$ and the
features of $i$.

A VACA intervention is performed by removing all the edges towards the
intervened node $i$, such that ${\mathcal{N}}_{i}=\varnothing$, while the rest
of the edges remains untouched.

In a VACA abduction step, the posterior distribution factorizes as

$\displaystyle
q_{\phi}({\mathbf{Z}}\mid{\mathbf{X}})=\prod_{i}q_{\phi_{i}}(Z_{i};{\boldsymbol{\eta}}_{i}^{enc}),$

where the distribution parameters
${\boldsymbol{\eta}}_{i}^{enc}={\boldsymbol{\eta}}_{i}^{enc}({\mathbf{X}}_{{{{\mathrm{pa}^{*}}(i)}}})$
are a function of all (and only) the features of node $i$ and the features of
its the parents.

#### Notation.

Consider a causal graph ${\mathcal{G}}:=({\mathbf{X}},{\mathbf{A}})$, which is
a directed acyclic graph (DAG). Let us define a path of length $n$ from node
$u$ to node $v$ in ${\mathcal{G}}$ as
$p(u,v)=(u,w_{1},w_{2},\dots,w_{n-1},v)$, which is an ordered sequence of
unique nodes such that i) there exists an edge in ${\mathcal{G}}$ between
concurrent nodes, ii) the first node is $u$, iii) and the last node is $v$. We
refer to the length of the path as $|p(u,v)|$, i.e., the number of edges in
the path, or alternatively, the number of nodes minus one. Let us define
$P(u,v)$ as the set of unique paths connecting $u$ to $v$. Let us define the
shortest path from $u$ to $v$ as $p^{-}(u,v)$ (i.e. the path with the minimum
number of edges to go from $u$ to $v$) and its length as
$d^{-}(u,v)=|p^{-}(u,v)|$. Let us define the longest path from $u$ to $v$ as
$p^{+}(u,v)$ (i.e. the path with the maximum number of edges to go from $u$ to
$v$) and its length as $d^{+}(u,v)=|p^{+}(u,v)|$. Let us define the set of
ancestors of node $i$ (i.e., ${{{\mathrm{an}(i)}}}$) as the set of nodes with
paths to $i$, i.e., $\\{j\mid|P(j,i)|>0\\}$. As for a GNN, we define the
number of hidden layers (total number of layers minus one) as ${N_{h}}$.

Then, we define the diameter ${\delta}$ of the graph ${\mathcal{G}}$ to be the
length of the longest shortest path and ${\gamma}$ to be the length of the
longest path of the graph, which we compute as

$\displaystyle{\delta}=\max_{u,v\in{\mathcal{G}}}d^{-}(u,v)\quad\text{ and
}\quad{\gamma}=\max_{u,v\in{\mathcal{G}}}d^{+}(u,v).$

######  Lemma 1.

A message passing Graph Neural Network (GNN) has at least ${N_{h}}$ hidden
layers if and only if the output feature of every node $i$ (i.e.,
$h_{i}^{{N_{h}}+1}$) receives information from any other node $j$ via paths
$p(j,i)$ such that $|p(j,i)|\leq{N_{h}}+1$.

###### Proof.

Step 1. The statement is that the feature of every node $i$ at the output of a
GNN with ${N_{h}}$ hidden layers (i.e., $h_{i}^{{N_{h}}+1}$) receives
information from any other node $j$ via paths $p(j,i)$ such that
$|p(j,i)|\leq{N_{h}}+1$. We give a proof by induction on ${N_{h}}$ for an
arbitrary node $i$, with input feature to the GNN $h_{i}^{0}$ and output
feature $h_{i}^{{N_{h}}+1}$.

Base case: The statement holds for ${N_{h}}=0$. By definition, a message
passing GNN with one layer only exchanges messages between neighboring nodes.
Hence, i) the output feature of node $i$ is
$h_{i}^{1}=f(\\{h_{i}^{0}\\}\cup\\{h_{j}^{0}|j\in{\mathcal{N}}_{i}\\};\theta)$,
which is only a function of the $1$-hop ancestors (i.e., parents); ii)
information is exchanged via paths that fulfill $|p(j,i)|\leq 1$.

Inductive step: We assume the statement holds for ${N_{h}}=k-1$. In this case,
i) the output feature of node $i$ is
$h_{i}^{k}=f(\\{h_{i}^{k-1}\\}\cup\\{h_{j}^{k-1}|j\in{\mathcal{N}}_{i}\\};\theta)$,
which is a function of the $k$-hop ancestors; ii) information is exchanged via
paths that fulfill $|p(j,i)|\leq k$. For ${N_{h}}=k$, the output feature of
node $i$ is
$h_{i}^{k+1}=f(\\{h_{i}^{k}\\}\cup\\{h_{j}^{k}|j\in{\mathcal{N}}_{i}\\};\theta)$.
Since $h_{j}^{k}$ is a function of $k$-hop ancestors of node $j$, it follows
that $h_{i}^{k+1}$ is a function of the ancestors of $h_{j}^{k}$ and on
parents of node $i$, i.e., the output feature of node $i$ is a function of its
$(k+1)$-ancestors. Then, it follows that for ${N_{h}}=k$, information is
exchanged via paths that fulfill $|p(j,i)|\leq k+1$.

Step 2. We assume that the output feature of every node $i$ receives
information from any other node $j$ via paths $p(j,i)$ such that
$|p(j,i)|\leq{N_{h}}+1$. Then, there exist node $i$ and $j$ such that
$|p(j,i)|={N_{h}}+1$. Then, by definition, a message passing GNN needs at
least ${N_{h}}$ hidden layers to capture paths $p(j,i)$ with length
$|p(j,i)|\leq{N_{h}}+1$. ∎

Remark. Lemma 1 implies that, for a given GNN with ${N_{h}}$ hidden layers,
the set of paths through which the output feature of node $i$ is a function of
node $j$ is

$\displaystyle\begin{split}P_{\text{GNN}}(j,i)=&\\{p(j,i)\mid p(j,i)\in
P(j,i)\\\ &\text{ and }|p(j,i)|\leq{N_{h}}+1\\}.\end{split}$ (9)

Additionally, if $|P_{\text{GNN}}(j,i)|=\varnothing$ then the output feature
of node $i$ is not a function of node $j$.

#### Proposition 1 (Causal factorization).

VACA satisfies causal factorization,
${p_{\theta}}({\mathbf{X}}\mid{\mathbf{Z}},{\mathbf{A}})=\prod_{i}p_{\theta_{i}}(X_{i}\mid{\mathbf{Z}}_{{{{\mathrm{an}^{*}}(i)}}})$,
if and only if the number of hidden layers in the decoder is greater or equal
than ${\delta}-1$, with ${\delta}$ being the length of the longest shortest
path between any two endogenous nodes.

###### Proof.

Consider a causal graph ${\mathcal{G}}:=({\mathbf{X}},{\mathbf{E}})$ with
diameter ${\delta}$ and a GNN decoder with ${N_{h}}$ hidden layers. We assume
VACA to satisfy causal factorization, i.e., ${\boldsymbol{\eta}}_{i}$ is a
function ${\mathbf{Z}}_{{{{\mathrm{an}^{*}}(i)}}}$ for all $i$. Therefore,
there exist node $i$ and $j$ such that $d^{-}(j,i)={\delta}$ (notice that this
implies that $j$ is an ancestor of $i$). Thus, by Lemma 1, the GNN decoder has
${N_{h}}\geq{\delta}-1$ hidden layers. The converse is true because Lemma 1 is
a bi-conditional statement. ∎

#### Proposition 2 (Causal interventions).

VACA captures causal interventions if and only if the number of hidden layers
in its decoder is greater than or equal to $\gamma-1$, with $\gamma$ being the
length of the longest path between any two endogenous nodes in
${\mathcal{G}}$.

A causal intervention involves to severe all the incoming edges to the
intervened nodes. Thus, VACA can only capture causal interventions, if it can
model all the causal paths, i.e., $P_{\text{GNN}}(j,i)=P(j,i){{{{\text{
}{\forall}}}}}i,j$. Otherwise, severing some paths will have no effect in the
resulting intervention, as we prove next.

###### Proof.

Consider a causal graph ${\mathcal{G}}:=({\mathbf{X}},{\mathbf{E}})$ with
length of the longest path between two nodes ${\gamma}$ and a GNN decoder with
${N_{h}}$ hidden layers. We assume that the GNN decoder models all the causal
paths, i.e., $P_{\text{GNN}}(j,i)=P(j,i){{{{\text{ }{\forall}}}}}i,j$. By
definition of ${\gamma}$, there exists at least one node $i$ with an ancestor
$j$ such that $d^{+}(j,i)={\gamma}$. Thus, by Lemma 1, the GNN decoder has
${N_{h}}\geq{\gamma}-1$ hidden layers. The converse is true because Lemma 1 is
a bi-conditional statement. ∎

#### Proposition 3 (Abduction).

The abduction step of an observed sample
${\mathbf{x}}=\\{x_{1},\ldots,x_{d}\\}$ in VACA satisfies that for all $i$ the
posterior of $Z_{i}$ is independent on the subset
$\\{x_{j}\\}_{j{{\not\in}}{{{\mathrm{pa}^{*}}(i)}}}\subseteq{\mathbf{x}}$, if
and only if the encoder GNN has no hidden layers.

###### Proof.

Consider a causal graph ${\mathcal{G}}:=({\mathbf{X}},{\mathbf{E}})$ and a GNN
encoder with ${N_{h}}$ hidden layers. We assume the posterior of $Z_{i}$ is
independent on the subset
$\\{x_{j}\\}_{j{{\not\in}}{{{\mathrm{pa}^{*}}(i)}}}\subseteq{\mathbf{x}}$,
i.e., the parameters ${\boldsymbol{\eta}}^{enc}_{i}$ (the output of the GNN)
is a function of $\\{x_{j}\\}_{j\in{{{\mathrm{pa}^{*}}(i)}}}$. Then, the GNN
only models paths $p(j,i)$ such that $d^{+}(j,i)=1$. It follows, by Lemma 1,
that the number of hidden layers of the encoder GNN is ${N_{h}}=0$. The
converse is true because Lemma 1 is a bi-conditional statement. ∎

${N_{h}}$ |  collider (${\delta}=1$, ${\gamma}=1$) |  triangle (${\delta}=1$, ${\gamma}=2$) |  chain (${\delta}=2$, ${\gamma}=2$)
---|---|---|---
MMD Obs. (%) | MMD Inter.(%) | MMD Obs.(%) | MMD Inter.(%) | MMD Obs.(%) | MMD Inter.(%)
0 | 1.86 $\pm$ 0.84 | 1.16 $\pm$ 0.86 | 21.76 $\pm$ 5.80 | 48.47 $\pm$ 12.57 | 9.69 $\pm$ 1.92 | 17.24 $\pm$ 2.82
1 | 1.31 $\pm$ 0.36 | 0.75 $\pm$ 0.17 | 9.17 $\pm$ 2.39 | 19.91 $\pm$ 4.75 | 6.44 $\pm$ 1.52 | 10.36 $\pm$ 2.67
2 | 1.32 $\pm$ 0.41 | 1.09 $\pm$ 1.08 | 7.19 $\pm$ 3.14 | 14.10 $\pm$ 6.45 | 4.35 $\pm$ 1.77 | 6.55 $\pm$ 1.64
Table 1:  Evaluation of the observational and interventional distributions
generated by VACA with different numbers of hidden layers ${N_{h}}$ in the
decoder. All metrics are multiplied by 100 (%).

##  Appendix C VACA implementation details

In this section, we extend Section 4.4 and provide further details about the
implementation of VACA for complex real-world datasets and causal graphs.

###  C.1 Heterogeneous endogenous variables

As described in Appendix A, each layer $l$ of a GNN uses the same parameters
$\theta=\\{\theta_{m}^{l},\theta_{u}^{l}\\}$ (corresponding to $f^{m}$ and
$f^{u}$) to update the features of every node, i.e.,
$h^{l}_{i}=f(\\{h_{i}^{l-1}\\}\cup\\{h_{j}^{l-1}\mid
j\in{{\mathcal{N}_{i}}}\\};\theta)$ with $f_{\theta}$ being reused for all
$i$. Nonetheless, the structural equations of an SCM define a unique function
$f_{i}$ for each node (see Property 1). To mimic this behavior, we will rely
on _port numbering_. In particular, for a given causal graph ${\mathcal{G}}$,
we uniquely identify each node with an index $i$ and each edge with the pair
of indexes of the nodes it connects.

Then, we define a _disjoint GNN layer_ by the following characteristics:

  * •

The node indexes define unique update functions $f^{u}_{i}$ for each node,
with parameters $\theta_{ui}$.

  * •

The edge indexes define unique message functions $f^{m}_{ij}$ for each edge,
with parameters $\theta_{mij}$.

Consequently, parameters are not shared among nodes and we can mimic the
diversity of the structural equations of an SCM and model heterogeneous
endogenous variables. In our GitHub
repository222https://github.com/psanch21/VACA, we present a PyTorch Geometric
implementation of the _disjoint GNN layer_.

###  C.2 Heterogeneous causal nodes

Assume an SCM with $d$ endogenous variables. As described in Section 4.4, is
it possible to model an endogenous variable ${\mathbf{X}}_{i}$ of the SCM as a
heterogeneous node, i.e., ${\mathbf{X}}_{i}=\\{X_{i1},\dots,X_{ik_{i}}\\}$,
where $k_{i}$ is the number of random variables in node $i$. In this section
we describe the implications this has on the design of the encoder and decoder
of VACA.

#### Implications for the decoder.

Given the heterogeneous nature of the nodes, the likelihood of VACA factorizes
as follows

$\displaystyle\begin{split}&p_{\theta}({\mathbf{X}}\mid{\mathbf{Z}})=\prod_{i=1}^{d}p_{\theta_{i}}({\mathbf{X}}_{i};{\boldsymbol{\eta}}_{i})=\prod_{i=1}^{d}\prod_{j=1}^{k_{i}}p_{\theta_{i}}(X_{ij};{\boldsymbol{\eta}}_{ij})\\\
&\text{ where
}{\boldsymbol{\eta}}_{i}={\boldsymbol{\eta}}_{i}({\mathbf{Z}}_{{{{\mathrm{an}^{*}}(i)}}})\text{
and
}{\boldsymbol{\eta}}_{ij}={\boldsymbol{\eta}}_{ij}({\mathbf{Z}}_{{{{\mathrm{an}^{*}}(i)}}}).\end{split}$

Figure 4: Heterogeneous VACA decoder architecture. Figure 5: Heterogeneous
VACA encoder architecture.

Note, each $p(X_{ij};{\boldsymbol{\eta}}_{ij})$ can be model with a different
distribution, e.g., Gaussian or categorical. This means that the likelihood
parameters ${\boldsymbol{\eta}}_{ij}$ for each random variable $X_{ij}$ may be
different depending on its distribution type. However, the decoder GNN
transforms the latent features into different features
$\mathbf{H}\in{\mathbb{R}}^{d\times F_{h}}$, where the feature vector of each
node $i$ has the same dimensionality $F_{h}$. As a consequence, $\mathbf{H}$
cannot model the diversity in the likelihood parameters
${\boldsymbol{\eta}}_{i}$. To overcome such limitation, we add at the output
of the GNN decoder a neural network (NN) per node $i$ with parameters
$\theta_{i}^{dec}$. Such a NN transforms $h_{i}$, i.e. the output features of
node $i$, into the set of likelihood parameters of each node
${\boldsymbol{\eta}}_{i}=\\{{\boldsymbol{\eta}}_{ij}\\}^{k_{i}}_{j=1}$, such
that the likelihood parameters of each random variable
${\boldsymbol{\eta}}_{ij}$ satisfy the constraints of the corresponding
likelihood $p(X_{ij};{\boldsymbol{\eta}}_{ij})$ (e.g., non-negativity of
variance for a Gaussian distribution). See Figure 4 for an illustration.

#### Implications for the encoder.

Due to the heterogeneous nature of nodes, each endogenous variable
${\mathbf{X}}_{i}=\\{X_{i1},\dots,X_{ik_{i}}\\}$, can have a different number
of random variables $k_{i}$ and thus the node $i$ corresponding to it in the
GNN will have features of different dimensions. However, as described in
Appendix A, a GNN takes as input in general a matrix feature
${\mathbf{X}}\in\mathbb{R}^{d\times F_{x^{enc}}}$. This implies, the features
of every node share the same dimensionality $F_{x^{enc}}$. To overcome this
limitation, we include for each node $i$ a neural network (NN) with parameters
$\theta^{{enc}}_{i}$ that transforms the corresponding heterogeneous random
variable ${\mathbf{X}}_{i}$ into a feature vector with the dimension
$F_{x^{enc}}$. See Figure 5 for an illustration.

##  Appendix D Validating & understanding VACA

Here, we empirically validate the design conditions of VACA and provide an
analysis of the effect of the latent space dimension $\dim\mathbf{z}$ on the
performance.

#### Validating VACA design conditions.

In a first step we empirically validate our design choices for the VACA
encoder and decoder. We show how the number of hidden layers $N_{h}$ in the
decoder affects the quality of the estimation of the observational and
interventional distribution. We do so for three SCMs, with different values of
longest shortest directed path $\delta$ and longest directed path $\gamma$.
Our observations in Table 1 match our expectations. The collider
($\delta=\gamma=1$) does not need any hidden layer to provide accurate
estimate of both the observational and interventional distributions. In
contrast, the triangle ($\delta=1,\gamma=2$) – in accordance with Proposition
2 – needs at least one hidden layer to get a more accurate estimate of the
interventional distribution. Finally, as stated by Propositions 1 and 2, the
chain ($\delta=\gamma=2$) requires at least one hidden layer to accurately
approximate both the observational and interventional distributions.

#### Analysis of the latent space dimension.

Here we present an analysis of the performance of VACA with respect to the
dimension of the latent space $\dim\mathbf{z}$. Without loss of generality, we
focus the analysis on the collider, triangle and chain graphs, and the NLIN
and NADD structural equations. The results are depicted in Figure 6. All
results are averaged over 10 different initializations. Regarding the
observational distribution (left column), we observe that VACA overfits as we
increase the $\dim\mathbf{z}$ for all SEMs and graphs under consideration.
Specially, the performance degrades notoriously with $\dim\mathbf{z}=32$. As
shown in Table 2, the number of parameters of VACA increases linearly with
$\dim\mathbf{z}$. As for the interventional distribution (middle column), we
observe similar behavior: for large values of $\dim\mathbf{z}$ performance
decreases and variance increases. For counterfactuals (right column), we
observe a similar behavior but considerably less pronounced. In summary, we
encourage practitioners to keep $\dim\mathbf{z}$ small to avoid overfitting
and obtain better and more consistent performance.

Figure 6: Performance of VACA varying the dimension of the latent space $\dim\boldsymbol{z}\in\\{2,4,8,16,32\\}$ for the non-linear (top row) and non-additive (bottom column) SEMs of the collider, triangle and chain graphs. We evaluate the quality of the observational (MMD Obs.) (left column), interventional (MMD Inter.) (middle column), and counterfactual (MSE CF) (right column) distributions. | collider | triangle | chain
---|---|---|---
$\dim\mathbf{z}$ | NLIN | NADD | NLIN | NADD | NLIN | NADD
2 | 3785 | 3785 | 4542 | 4542 | 3785 | 3785
4 | 4445 | 4445 | 5334 | 5334 | 4445 | 4445
8 | 5765 | 5765 | 6918 | 6918 | 5765 | 5765
16 | 8405 | 8405 | 10086 | 10086 | 8405 | 8405
32 | 13685 | 13685 | 16422 | 16422 | 13685 | 13685
Table 2: Number of parameters of VACA for different $\dim\mathbf{z}$ and
datasets

##  Appendix E Experiments: setting, metrics and further results

This section provides a complete description of the experimental set-up,
including the (semi-)synthetic datasets (Section E.1), training of VACA,
MultiCVAE [16] and CAREFL [17] (Section E.2), metrics reported in the
experiments (Section E.3), additional results (Section E.4), complexity of the
algorithm (Section E.5), and computing infrastructure (Section E.6).

###  E.1 Datasets

The following (semi-)synthetic datasets are taken from or inspired by [16].
The distribution of exogenous variables $p({\mathbf{U}})$ for triangle, chain
and collider follows Table 3 with with $\operatorname{MoG}$ denoting a mixture
of Gaussian distributions.

SCM | $p(U_{1})$ | $p(U_{2})$ | $p(U_{3})$
---|---|---|---
LIN | $\operatorname{MoG}(0.5\mathcal{N}(-2,1.5)+0.5\mathcal{N}(1.5,1))$ | $\mathcal{N}(0,1)$ | $\mathcal{N}(0,1)$
NLIN | $\operatorname{MoG}(0.5\mathcal{N}(-2,1.5)+0.5\mathcal{N}(1.5,1))$ | $\mathcal{N}(0,0.1)$ | $\mathcal{N}(0,1)$
NADD | $\operatorname{MoG}(0.5\mathcal{N}(-2.5,1)+0.5\mathcal{N}(2.5,1))$ | $\mathcal{N}(0,0.25)$ | $\mathcal{N}(0,0.0625)$
Table 3: Distribution of exogenous variables $p({\mathbf{U}})$ for SCM triangle, chain, collider. | SCM | ${{{\tilde{f}}}}_{1}:=X_{1}$ | ${{{\tilde{f}}}}_{2}:=X_{2}$ | ${{{\tilde{f}}}}_{3}:=X_{3}$
---|---|---|---|---
collider | LIN | $U_{1}$ | $U_{2}$ | $0.05X_{1}+0.25X_{2}+U_{3}$
NLIN | $U_{1}$ | $U_{2}$ | $0.05X_{1}+0.25(X_{2})^{2}+U_{3}$
NADD | $U_{1}$ | $U_{2}$ | $-1+0.1\operatorname{sgn}(U_{3})((X_{1})^{2}+(X_{2})^{2})U_{3}$
triangle | LIN | $U_{1}$ | $-X_{1}+U_{2}$ | $X_{1}+0.25X_{2}+U_{3}$
NLIN | $U_{1}$ | $-1+\frac{3}{(1+\operatorname{exp}(-2X_{1}))}+U_{2}$ | $X_{1}+0.25(X_{2})^{2}+U_{3}$
NADD | $U_{1}$ | $0.25\operatorname{sgn}(U_{2})*(X_{1})^{2}(1+(U_{2})^{2})$ | $-1+0.1\operatorname{sgn}(U_{3})((X_{1})^{2}+(X_{2})^{2})+U_{3}$
chain | LIN | $U_{1}$ | $-X_{1}+U_{2}$ | $0.25*X_{2}+U_{3}$
NLIN | $U_{1}$ | $-1+\frac{3}{(1+\operatorname{exp}(-2X_{1}))}+U_{2}$ | $0.25*(X_{2})^{2}+U_{3}$
NADD | $U_{1}$ | $0.25\operatorname{sgn}(U_{2})(X_{1})^{2}(1+(U_{2})^{2})$ | $-1+0.1\operatorname{sgn}(U_{3})((X_{2})^{2})+U_{3}$
Table 4: Structural equations ${{{\tilde{\mathbf{F}}}}}$ for different SCMs
with ${\mathbf{U}}\sim p({\mathbf{U}})$ in Table 3. Function
$\operatorname{sgn}(x)$ returns an element-wise indication of the sign of $x$.

#### Collider.

The collider is a synthetic dataset, which consists of 3 endogenous variables.
The structural equations are shown in Table 4. Figure 7 illustrates the
corresponding causal graph with $d=|{\mathbf{X}}|=3$ nodes, diameter
${\delta}=1$ and longest path ${\gamma}=1$.

$X_{3}$$X_{2}$$X_{1}$ Figure 7: Causal graph for variables ${\mathbf{X}}$ of
SCM collider.

#### Triangle.

The triangle is a synthetic dataset, which consists of 3 endogenous variables.
The structural equations are shown in Table 4. Figure 8 illustrates the
corresponding causal graph with $d=|{\mathbf{X}}|=3$ nodes, diameter
${\delta}=1$ and longest path ${\gamma}=2$.

$X_{1}$$X_{3}$$X_{2}$ Figure 8: Causal graph for variables ${\mathbf{X}}$ of
SCM triangle.

#### Chain.

The chain is a synthetic dataset, which consists of 3 endogenous variables.
The structural equations are shown in Table 4. Figure 9 illustrates the
corresponding causal graph with $d=|{\mathbf{X}}|=3$ nodes, diameter
${\delta}=2$ and longest path ${\gamma}=2$.

$X_{1}$$X_{2}$$X_{3}$ Figure 9: Causal graph for variables ${\mathbf{X}}$ of
SCM chain.

#### M-graph.

The M-graph is a synthetic dataset, which consists of 5 endogenous variables.
Here, the distributions of exogenous variables follow $U_{i}\sim
p(U_{i})=\mathcal{N}(0,1){{{{\text{ }{\forall}}}}}i\in{1\dots 5}$. The
structural equations are shown in Table 5 and Figure 10 illustrates the
corresponding causal graph with $d=|{\mathbf{X}}|=5$ nodes, diameter
${\delta}=1$ and longest path ${\gamma}=1$.

$X_{4}$$X_{2}$$X_{1}$$X_{5}$$X_{3}$ Figure 10: Causal graph for variables ${\mathbf{X}}$ of SCM M-graph. SCM | ${{{\tilde{f}}}}_{1}:=X_{1}$ | ${{{\tilde{f}}}}_{2}:=X_{2}$ | ${{{\tilde{f}}}}_{3}:=X_{3}$ | ${{{\tilde{f}}}}_{4}:=X_{4}$ | ${{{\tilde{f}}}}_{5}:=X_{5}$
---|---|---|---|---|---
LIN | $U_{1}$ | $U_{2}$ | $X_{1}+U_{3}$ | $-X_{2}+0.5X_{1}+U_{4}$ | $-1.5X_{2}+U_{5}$
NLIN | $U_{1}$ | $U_{2}$ | $X_{1}+0.5(X_{1})^{2}+U_{3}$ | $-X_{2}+0.5(X_{1})^{2}+U_{4}$ | $-1.5(X_{2})^{2}+U_{5}$
NADD | $U_{1}$ | $U_{2}$ | $X_{1}*U_{3}$ | $(-X_{2}+0.5*(X_{1})^{2})U_{4}$ | $(-1.5(X_{2})^{2})U_{5}$
Table 5: Structural Equations ${{{\tilde{\mathbf{F}}}}}$ for SCM M-graph with
$U_{i}\sim p(U_{i})=\mathcal{N}(0,1)\,{{{{\text{ }{\forall}}}}}i\in{1\dots
5}$.

(a) collider LIN

(b) collider NLIN

(c) collider NADD

(d) triangle LIN

(e) triangle NLIN

(f) triangle NADD

(g) chain LIN

(h) chain NLIN

(i) chain NADD

(j) M-graph LIN

(k) M-graph NLIN

(l) M-graph NADD

Figure 11: Histograms and scatter plots of pairwise feature relations for the
collider, triangle, chain and M-graph datasets.

#### Loan.

The loan is a semi-synthetic dataset from [16], which reflects a loan approval
setting in the real-world inspired by German Credit dataset [7]. It consists
of 7 endogenous variables: gender $G$, age $A$, education $E$, loan amount
$L$, loan duration $D$, income $I$ and savings $S$ with the following
structural equations and distributions of exogenous variables:

$\displaystyle\begin{split}f_{G}:&G=U_{G}\\\ f_{A}:&A=-35+U_{A}\\\
f_{E}:&E=-0.5+\left(1+e^{+1-0.5G-\left(1+e^{-0.1A}\right)^{-1}-U_{E}}\right)^{-1}\\\
f_{L}:&L=1+0.01(A-5)(5-A)+G+U_{L}\\\ f_{D}:&D=-1+0.1A+2G+L+U_{D}\\\
f_{I}:&I=-4+0.1(A+35)+2G+GE+U_{I}\\\
f_{S}:&S=-4+1.5\mathbb{I}_{\\{I>0\\}}I+U_{S}\end{split}$

with $U_{G}~{}\sim\operatorname{Bernoulli}(0.5)$,
$U_{A}~{}\sim\operatorname{Gamma}(10,3.5)$, $U_{E}~{}\sim\mathcal{N}(0,0.25)$,
$U_{L}~{}\sim~{}\mathcal{N}(0,4)$, $U_{D}~{}\sim~{}\mathcal{N}(0,9)$,
$U_{S}~{}\sim~{}\mathcal{N}(0,25)$, $U_{I}~{}\sim~{}\mathcal{N}(0,4)$.

Note, the authors model variables w.r.t. their relative meaning in terms of
deviation from the mean. See [16] for further details. Figure 12 illustrates
the corresponding causal graph with $d=|{\mathbf{X}}|=7$ nodes, diameter
${\delta}=2$ and longest path ${\gamma}=3$.

$G$$A$$L$$D$$E$$I$$S$ Figure 12: Causal graph for variables ${\mathbf{X}}$ of
SCM loan. Figure 13: Histograms and scatter plots of pairwise feature
relations for the loan dataset.

#### Adult.

We introduce the new semi-synthetic adult dataset which aims to reflect the
relationship between the variables that affect the annual income of a person
inspired by the real-world Adult dataset [6]. Figure 15 shows the
corresponding causal graph with $d=|{\mathbf{X}}|=11$ nodes, diameter
${\delta}=2$ and longest path ${\gamma}=6$. We base the causal graph on [3].
Our dataset consists of the following 11 endogenous variables:

  * •

Race $R$ is an independent value.

  * •

Age $A$ is an independent value.

  * •

Native country ($N$) is an independent value.

  * •

Sex $S$ is an independent value.

  * •

Education level ($E$) depends on the sex, the native county and the race.

  * •

Hours (worked) per week $H$ depends additionally on the race, the native
country and the sex. We consider this continuous variable to be between 0 and
80 hours.

  * •

Work status $W$, i.e., being unemployed or self-employed, is a discrete
variable that depends on the number of hours worked per day, age, native
country, and level of education.

  * •

Martial status $M$, i.e., married never married or separated, depends on age,
race, work status, hours per week, native country and sex. Individuals under a
certain age are never married.

  * •

Occupation sector $O$ is a discrete variable with values technology , social
sciences and medicine. It depends on the race, age, education, martial status
and sex.

  * •

Relationship status $L$ is a discrete variable, with values wife, own-child,
husband, not-in-family, unmarried. It depends on martial status, education,
age, native country, and sex.

  * •

Income $I$ depends on race, age, education, occupation, work status, martial
status, hours per week, relationship status, native country, and sex.

These aspects are reflected in the following structural equations :

$\displaystyle\begin{split}&f_{R}:R=U_{R}\\\ &f_{A}:A=U_{A}+17\\\
&f_{N}:N=U_{N}\\\ &f_{S}:S=U_{S}\\\
&f_{E}:E=\exp\left(2\mathbb{I}_{\\{R=0\\}})+\mathbb{I}_{\\{R=1\\}})+\sigma(A-30)\right)\\\
&\quad+(0.5\mathbb{I}_{\\{S=0\\}}+\mathbb{I}_{\\{S=1\\}})(2\mathbb{I}_{\\{N=1\\}}+5\mathbb{I}_{\\{N=2\\}}\\\
&\quad+\mathbb{I}_{\\{N=3\\}})+U_{E}\\\
&f_{H}:H=((40\mathbb{I}_{\\{N=0\\}}+36\mathbb{I}_{\\{N=1\\}}+50\mathbb{I}_{\\{N=2\\}}\\\
&\quad+30\mathbb{I}_{\\{N=3\\}})*\left(0.5\mathbb{I}_{\\{R=0\\}}+\mathbb{I}_{\\{R=1\\}}+1.3\mathbb{I}_{\\{R=2\\}}\right)\\\
&\quad+2\exp\left(-(A-30)^{2}\right)+5\mid\operatorname{tanh}(E-2)\mid\\\
&\quad+2\mathbb{I}_{\\{S=0\\}}+U_{H})\mathbb{I}_{\\{A<70\\}}\\\
&f_{W}:W=w_{2}\mathbb{I}_{\\{w_{2}>=0\\}},\\\ &\quad
w_{2}=w_{1}\mathbb{I}_{\\{w_{1}<=3\\}}+3\mathbb{I}_{\\{w_{1}>3\\}}\\\ &\quad
w_{1}=\mathbb{I}_{\\{5\mid\operatorname{tanh}(E-2)\mid+\sigma(H-30+U_{W})>0.3\\}}\\\
&\quad+\mathbb{I}_{\\{\sigma(H-30+U_{W}>0.3\\}}\mathbb{I}_{\\{(A+1.5U_{W})>50\\}}\\\
&\quad-\mathbb{I}_{\\{N=0\\}}+\mathbb{I}_{\\{N=2\\}}+3\mathbb{I}_{\\{N=3\\}})\\\
&f_{M}:M=\operatorname{mode}(r_{2},a_{2},W,h_{2},H,g_{3})\\\ &\quad
r_{2}=2*\mathbb{I}_{\\{r_{1}=1\\}}+\mathbb{I}_{\\{r_{1}=2\\}}\\\ &\quad
r_{1}=\operatorname{int}(R+0.2U_{M})\mathbb{I}_{\\{R\in[0,2]\\}}+2\mathbb{I}_{\\{R>2\\}}\\\
&\quad g_{3}=\mathbb{I}_{\\{g_{2}=0\\}}+2*\mathbb{I}_{\\{g_{2}=1\\}}\\\ &\quad
g_{2}=0*\mathbb{I}_{\\{g_{1}<0\\}}+\mathbb{I}_{\\{g_{1}>1\\}}+g_{1}\mathbb{I}_{\\{g_{1}\in[0,1]\\}}\\\
&\quad g_{1}=\operatorname{int}(G+0.5U_{M})\\\ \end{split}$

$\displaystyle\begin{split}&\quad
a_{2}=2\mathbb{I}_{\\{a_{1}\in(20,40]\\}}+\mathbb{I}_{\\{a_{1}\in(40,50]\\}}+2\mathbb{I}_{\\{a_{1}>=50\\}}\\\
&\quad a_{1}=A+2U_{M}\\\ &\quad
h_{2}=h_{1}\mathbb{I}_{\\{h_{1}<=2\\}}+2\mathbb{I}_{\\{h_{1}>2\\}}\\\ &\quad
h_{1}=3\operatorname{int}(\sigma(H-30))\\\
&f_{L}:L=0\mathbb{I}_{\\{(M=1)\wedge(c<-1)\\}}+1\mathbb{I}_{\\{(M=1)\wedge(c\geq-1)\\}}\\\
&\quad+2\mathbb{I}_{\\{(M\neq 1)\wedge(c\geq-1)\\}}+1\mathbb{I}_{\\{(M\neq
1)\wedge(c<-1)\\}}\\\ &\quad
c=c_{n}+c_{e}+2\mathbb{I}_{\\{A<20\\}}-2\mathbb{I}_{\\{S=0\\}}\\\ &\quad
c_{n}=U_{O}\mathbb{I}_{\\{N=0\\}}-U_{O}\mathbb{I}_{\\{N=1\\}}+2U_{O}\mathbb{I}_{\\{N=2\\}}+2\mathbb{I}_{\\{N=3\\}}\\\
&\quad c_{e}=\sigma(E-30)\\\
&f_{O}:O=0\mathbb{I}_{\\{k<1\\}}+1\mathbb{I}_{\\{1\leq k\leq
4\\}}+2\mathbb{I}_{\\{4<k\\}}\\\ &\quad k=R+k_{a}+k_{e}+W+3M+4S\\\ &\quad
k_{a}=2e^{-(A+U_{O}-20)^{2}}\\\ &\quad k_{e}=-\sigma(E*U_{O}-30)\\\
&f_{I}:I=U_{I}\\\
&\quad+10,000\mathbb{I}_{\\{R>1.5\\}}+20,000\mathbb{I}_{\\{R<1.5\\}}\\\
&\quad+3,000\mathbb{I}_{\\{21\leq A<30\\}}+8,000\mathbb{I}_{\\{30\leq A\\}}\\\
&\quad+5,000\mathbb{I}_{\\{E<2\\}}+10,000\mathbb{I}_{\\{2\leq
E<10\\}}+30,000\mathbb{I}_{\\{10\leq E\\}}\\\
&\quad+5,000\mathbb{I}_{\\{O=1\\}}+15,000\mathbb{I}_{\\{O=2\\}}\\\
&\quad+5,000\mathbb{I}_{\\{W=0\\}}+7,000\mathbb{I}_{\\{W=1\\}}\\\
&\quad+1,000\mathbb{I}_{\\{M=0\\}}+4,000\mathbb{I}_{\\{M=1\\}}-2,000\mathbb{I}_{\\{M=2\\}}\\\
&\quad+15,000\mathbb{I}_{\\{H>45\\}}+10,000\mathbb{I}_{\\{N\geq 2\\}}\\\
&\quad+4,000\mathbb{I}_{\\{S=1\\}}+3,000\mathbb{I}_{\\{R\leq 1\\}}\\\
\end{split}$

$R$$A$$N$$S$$H$$E$$W$$M$$L$$O$$I$ Figure 15: Causal graph for variables
${\mathbf{X}}$ of SCM adult. Figure 16:  Histograms and scatter plots of
pairwise feature relations for the adult dataset.

###  E.2 Training and cross-validation

This section details the hyperparameter cross-validation of VACA, MultiCVAE
[16] and CAREFL [17] for the experiments in Tables LABEL:tab:results and Table
10. Across experiments and models we generate synthetic datasets consisting of
5000 training samples, 2500 test samples and 2500 validation samples and we
use a batch size of 1000. Al the models have been cross-validated using a
similar computational budget, i.e. the number of combinations of
hyperparameters cross-validated is the same for all the models (see Table 6).
Additionally, we run each configuration over 10 different random
initializations. Regarding the number of iterations, we assumed a large enough
number of epochs for training loss to converge, and implemented early stopping
to prevent overfitting the training data. The best configuration has been
chosen in terms of the best observational MMD.

SCM | MultiCVAE | CAREFL | VACA |
---|---|---|---|---
chain | 72 | 72 | 72 |
collider | 72 | 72 | 72 |
triangle | 72 | 72 | 72 |
M-graph | 72 | 72 | 72 |
loan | 40 | 24 | 24 |
Table 6: Number of combinations of hyperparameters cross-validated for each of
the models and datasets. Note that for each combination we run 10 different
seeds.

#### VACA.

We optimize the ELBO [19] and use the IWAE [2] with $K=100$ as the objective
metric for the early stopping procedure. We use the Rectified Linear Unit
(ReLU) as activation function. We trained with a learning rate $\eta=0.005$
for a maximum of 500 epochs, or alternatively until the objective metric does
not improve in 50 epochs. Also, we regularize the training of VACA using a
novel _parents dropout:_ randomly removing all incoming edges to the nodes
with probability $p\in[0,1)$. In our experiments, we observe that adding this
regularization improves overall performance. We cross-validated the parents
dropout rate with values $\\{0.1,0.2\\}$, the number of hidden layers of the
decoder with values $\\{0,1,2,3,4,5\\}$ (the specific values depend on the
diameter of the graph) with 16 neurons each, and weather to use or not
residual connections. The best models are reported in Table 7. We use a latent
variable dimension of 4 and a Gaussian likelihood with a small variance
$\sigma^{2}=\lambda_{KLD}/2$ with $\lambda_{KLD}=0.05$.

SCM |  | Encoder Arch. | Decoder Arch. | Parents Dropout | Residual
---|---|---|---|---|---
chain | LIN | $1\times 16\times 4$ | $4\times 16\times 1$ | 0.1 | 1
NLIN | $1\times 16\times 4$ | $4\times 16\times 16\times 1$ | 0.1 | 0
NADD | $1\times 16\times 4$ | $4\times 64\times 1$ | 0.1 | 1
collider | LIN | $1\times 16\times 4$ | $4\times 16\times 16\times 1$ | 0.2 | 1
NLIN | $1\times 16\times 4$ | $4\times 16\times 1$ | 0.1 | 0
NADD | $1\times 16\times 4$ | $4\times 16\times 16\times 1$ | 0.1 | 1
triangle | LIN | $1\times 16\times 4$ | $4\times 16\times 16\times 1$ | 0.1 | 0
NLIN | $1\times 16\times 4$ | $4\times 16\times 16\times 1$ | 0.2 | 0
NADD | $1\times 16\times 4$ | $4\times 16\times 16\times 1$ | 0.1 | 0
M-graph | LIN | $1\times 16\times 4$ | $4\times 16\times 1$ | 0.1 | 0
NLIN | $1\times 16\times 4$ | $4\times 16\times 16\times 1$ | 0.1 | 0
NADD | $1\times 16\times 4$ | $4\times 16\times 16\times 1$ | 0.1 | 0
loan | - | $1\times 16\times 16\times 4$ | $4\times 16\times 16\times 16\times 1$ | 0.1 2 | 0
adult | - | $1\times 16\times 16\times 4$ | $4\times 8\times 8\times 8\times 8\times 1$ | 0.2 | 1
Table 7: Hyperparemeter selection for our VACA training for the SCMs on the synthetic datasets triangle, collider, chain and M-graph and on the semi-synthetic dataset loan. Note, the encoder architecture refers to the layers in function $f^{m}$, while the decoder architecture refers to the different GNN layers.  SCM |  | $\lambda_{\text{KLD}}$ | Encoder Arch. | Decoder Arch. | Dropout rate
---|---|---|---|---|---
chain | LIN | 0.05 | $1\times 16\times 4$ | $4\times 128\times 1$ | 0.0
NLIN | 0.05 | $1\times 16\times 1$ | $1\times 128\times 1$ | 0.0
NADD | 0.05 | $1\times 16\times 4$ | $4\times 128\times 1$ | 0.0
collider | LIN | 0.05 | $1\times 16\times 4$ | $4\times 1$ | 0.0
NLIN | 0.05 | $1\times 16\times 4$ | $4\times 1$ | 0.0
NADD | 0.05 | $1\times 16\times 4$ | $4\times 1$ | 0.0
triangle | LIN | 0.05 | $1\times 16\times 4$ | $4\times 128\times 1$ | 0.2
NLIN | 0.05 | $1\times 16\times 1$ | $1\times 128\times 1$ | 0.1
NADD | 0.05 | $1\times 16\times 4$ | $4\times 128\times 1$ | 0.2
M-graph | LIN | 0.05 | $1\times 16\times 4$ | $4\times 1$ | 0.2
NLIN | 0.05 | $1\times 16\times 4$ | $4\times 1$ | 0.1
NADD | 0.05 | $1\times 16\times 4$ | $4\times 32\times 32\times 1$ | 0.2
loan | - | 0.05 | $1\times 16\times 16\times 1$ | $1\times 64\times 64\times 1$ | 0.0
adult | - | 0.05 | $1\times 16\times 16\times 2$ | $2\times 16\times 1$ | 0.0
Table 8: Hyperparemeter selection for MultiCVAE [16] training for the SCMs on
the synthetic datasets with three nodes (i.e., triangle, collider and chain),
M-graph and for the semi-synthetic dataset loan.

#### MultiCVAE.

Karimi et al. [16] propose to train a conditional variational autoencoder
(CVAE) for each endogenous variable that is not a root node in the causal
graph, we refer to as MultiCVAE. Different from [16], our implementation also
models non-root nodes as CVAEs, since our goal is to model the joint
distribution, while [16] target interventional and counterfactual
distributions for algorithmic recourse only. Additionally, we perform the
necessary modifications for training on normalized data.

We cross-validated the number and size of hidden layers of the decoder, the
dropout rate, and the dimension of the latent space. The best models
(according to the observational MMD) are reported in Table 8. As with VACA, we
assume $\lambda_{KLD}=0.05$ for all SCMs and CVAEs.

#### CAREFL.

Khemakhem et al. [17] propose CAREFL, an autoregressive causal flows model for
causal discovery, which also allows to answer interventional and
counterfactual queries. The authors rely on real-valued non-volume preserving
(real NVP) transformations, since they mainly focus on the multivariate bi-
variate case. As this flow architecture is not suited for general graphs, we
use their framework with Neural Spline Autoregressive Flows [8]. We have cross
validated the number of flows $\\{2,3,4,5,6\\}$ and the number of hidden units
of the neural networks $\\{5,10,16,32,64,96\\}$. The final configuration–i.e.,
the configuration with the lowest observational MMD–is displayed in Table 9.

SCM |  | Flows | Hidden Units
---|---|---|---
chain | LIN | 4 | 5
NLIN | 3 | 64
NADD | 5 | 5
collider | LIN | 2 | 96
NLIN | 2 | 64
NADD | 3 | 16
triangle | LIN | 4 | 64
NLIN | 4 | 5
NADD | 3 | 96
M-graph | LIN | 4 | 64
NLIN | 4 | 96
NADD | 5 | 96
loan | - | 4 | 64
adult | - | 2 | 96
Table 9: Hyperparemeter selection for CAREFL [17] training for different SCMs.

###  E.3 Performance metrics

In the following we describe the metrics used to evaluate the performance of
VACA in Section 5. In all experiments we use (semi-)synthetic datasets with
access to samples from the ground truth distribution
$\\{{\mathbf{x}}_{i}\\}^{n}_{i=0}\sim P$ as well as from the estimated
distribution $\\{{\hat{\mathbf{x}}}_{i}\\}^{n}_{i=0}\sim Q$.

For the interventional and counterfactual distribution, we perform a set of
interventions
${{\mathcal{I}}}=\\{do({\mathbf{X}}_{{{\mathcal{I}}}_{j}}=\alpha_{j})\\}_{j}$,
where ${{\mathcal{I}}}_{j}\in[d]$ and
$\alpha_{j}\in\\{-1.0,-0.5,0.0,0.5,1.0\\}\times\sigma_{{{\mathcal{I}}}_{j}}$
with $\sigma_{{{\mathcal{I}}}_{j}}$ as the empirical standard deviation of the
intervened variable ${\mathbf{X}}_{{{\mathcal{I}}}_{j}}$ prior to intervention
(i.e., in the observational distribution). Note that we only intervene on one
variable at a time. For each intervention in ${{\mathcal{I}}}$, we are
interested in the estimated distribution of variables causally affected by the
intervention $\\{{\mathbf{X}}_{i}|i\in des({{\mathcal{I}}}_{j})\\}$, i.e., the
set of descendants of the variable intervened. Note that
$des({{\mathcal{I}}}_{j})$ refers to the set of indexes of the descendants. It
follows that we do not intervene on leaf nodes.

#### Mean Maximum Discrepancy (MMD).

The Mean Maximum Discrepancy (MMD) [13] is a kernel-based distance-measure
between two distributions $P$ and $Q$ on the basis of samples from both
distributions. The smaller the MMD, the more likely it is that the sets of
samples are drawn from the same distributions, i.e. the better distributions
match.

Without access to underlying distribution, we can compute an unbiased
empirical squared MMD estimate using a kernel function $k$ as[13]:

$\displaystyle{\widehat{\mathrm{MMD}}}^{2}({\mathbf{X}},\hat{{\mathbf{X}}})=\frac{1}{n(n-1)}\left(\sum_{i=1}^{n}\sum_{j=1}^{n}k\left({\mathbf{x}}_{i},{\mathbf{x}}_{j}\right)+\sum_{i=1}^{n}\sum_{j=1}^{n}k\left({\hat{\mathbf{x}}}_{i},{\hat{\mathbf{x}}}_{j}\right)-2\sum_{i=1}^{n}\sum_{j=1}^{n}k\left({\mathbf{x}}_{i},{\hat{\mathbf{x}}}_{j}\right)\right).$
(10)

In our implementation we use as kernel a mixture of RBF (Gaussian) kernels
with different bandwidths and sample size $n=1000$.

#### Estimation squared error for the mean (MeanE).

For the interventional distribution, we compute the estimation squared error
for the mean (MeanE) as the average (across interventions) of the squared
difference between the empirical means of the true and estimated
interventional distributions (for the descendants of the intervened
variables):

$\displaystyle\operatorname{\text{MeanE}}$
$\displaystyle=\frac{1}{|{{\mathcal{I}}}|}\sum_{{{\mathcal{I}}}_{j}\in{{\mathcal{I}}}}\frac{1}{|des({{\mathcal{I}}}_{j})|}\sum_{i\in
des({{\mathcal{I}}}_{j})}\left(E[X_{i}^{{{\mathcal{I}}}_{j}}]-E[{\hat{X}}_{i}^{{{\mathcal{I}}}_{j}}]\right)^{2}$

#### Estimation squared error for the standard deviation (StdE).

For the interventional distribution, we compute the estimation squared error
for the standard deviation (StdE) as the average (across interventions) of the
squared difference between the empirical standard deviation of the true
$\tilde{\sigma}(X_{i}^{{{\mathcal{I}}}_{j}})$ and estimated
$\tilde{\sigma}({\hat{X}}_{i}^{{{\mathcal{I}}}_{j}})$ interventional
distributions (for the descendants of the intervened variables):

$\displaystyle\operatorname{\text{StdE}}$
$\displaystyle=\frac{1}{|{{\mathcal{I}}}|}\sum_{{{\mathcal{I}}}_{j}\in{{\mathcal{I}}}}\frac{1}{|des({{\mathcal{I}}}_{j})|}\sum_{i\in
des({{\mathcal{I}}}_{j})}\left(\tilde{\sigma}(X_{i}^{{{\mathcal{I}}}_{j}})-\tilde{\sigma}({\hat{X}}_{i}^{{{\mathcal{I}}}_{j}})\right)^{2}$

#### Mean squared error (MSE).

For the counterfactual distribution, we compute the mean squared error (MSE)
as the average (across interventions) of the pairwise squared difference
between true and estimated counterfactual values for the descendants of the
intervened variable. More in detail, let us define the random variable
$T^{{{\mathcal{I}}}_{j}}$ as the Frobenius norm of the difference between true
$\mathbf{x}^{CF}_{des({{\mathcal{I}}}_{j})}$ and estimated
$\hat{\mathbf{x}}^{CF}_{des({{\mathcal{I}}}_{j})}$ counterfactual values for
the descendants of the intervened variable, i.e.,

$\displaystyle
T^{{{\mathcal{I}}}_{j}}=||\mathbf{x}^{CF}_{des({{\mathcal{I}}}_{j})}-\hat{\mathbf{x}}^{CF}_{des({{\mathcal{I}}}_{j})}||^{2}_{2},$
(11)

Thus, we can compute the counterfactual MSE as:

$\displaystyle\operatorname{\text{MSE}}=\frac{1}{|{{\mathcal{I}}}|}\sum_{{{\mathcal{I}}}_{j}\in{{\mathcal{I}}}}\frac{1}{|des({{\mathcal{I}}}_{j})|}E\left[T^{{{\mathcal{I}}}_{j}}\right]$
(12)

#### Standard deviation of the squared error (SSE).

Similarly, we can compute the average (across interventions) of the standard
deviation of the counterfactual squared error as:

$\displaystyle\operatorname{\text{SSE}}=\frac{1}{|{{\mathcal{I}}}|}\sum_{{{\mathcal{I}}}_{j}\in{{\mathcal{I}}}}\frac{1}{|des({{\mathcal{I}}}_{j})|}\sum_{i\in
des({{\mathcal{I}}})}\tilde{\sigma}\left(T^{{{\mathcal{I}}}_{j}}\right),$ (13)

where $\tilde{\sigma}\left(T^{{{\mathcal{I}}}_{j}}\right)$ denotes the
empirical standard deviation of $T^{{{\mathcal{I}}}_{j}}$.

###  E.4 Additional results

|  |  | Obs. | Interventional | Counterfactuals |
---|---|---|---|---|---|---
SCM | Model | MMD | MMD | MeanE | StdE | MSE | SSE | Num. params
collider |  LIN | MultiCVAE | 30.37$\pm$8.16 | 44.70$\pm$12.25 | 13.29$\pm$4.78 | 46.56$\pm$2.40 | 87.41$\pm$3.64 | 65.15$\pm$2.83 | 553
CAREFL | 9.27$\pm$1.49 | 4.86$\pm$0.45 | 0.35$\pm$0.08 | 81.89$\pm$1.78 | 8.11$\pm$0.58 | 7.83$\pm$0.55 | 6420
VACA | 1.50$\pm$0.67 | 1.57$\pm$0.41 | 0.75$\pm$0.31 | 41.99$\pm$0.30 | 9.86$\pm$0.74 | 7.06$\pm$0.38 | 5600
NLIN | MultiCVAE | 28.03$\pm$9.12 | 41.60$\pm$12.62 | 10.49$\pm$4.12 | 46.48$\pm$2.43 | 82.32$\pm$2.61 | 62.05$\pm$1.87 | 553
CAREFL | 10.38$\pm$2.00 | 4.69$\pm$0.38 | 0.19$\pm$0.07 | 80.68$\pm$2.08 | 6.93$\pm$0.40 | 7.15$\pm$0.64 | 4308
VACA | 0.95$\pm$0.27 | 0.97$\pm$0.23 | 0.26$\pm$0.12 | 42.20$\pm$0.24 | 5.01$\pm$0.73 | 4.08$\pm$0.54 | 1805
NADD | MultiCVAE | 29.72$\pm$8.90 | 117.67$\pm$46.20 | 26.23$\pm$11.09 | 39.75$\pm$1.34 | 73.93$\pm$10.44 | 51.44$\pm$6.06 | 553
CAREFL | 9.79$\pm$1.29 | 9.13$\pm$1.45 | 1.58$\pm$0.94 | 102.13$\pm$3.16 | 9.92$\pm$0.51 | 32.66$\pm$0.86 | 1710
VACA | 1.71$\pm$0.48 | 29.87$\pm$1.94 | 10.71$\pm$1.05 | 50.25$\pm$0.87 | 31.99$\pm$0.84 | 39.10$\pm$0.67 | 5600
triangle |  LIN | MultiCVAE | 33.12$\pm$3.89 | 157.50$\pm$15.08 | 12.10$\pm$2.03 | 44.66$\pm$2.57 | 65.95$\pm$3.05 | 42.39$\pm$1.75 | 3243
CAREFL | 13.50$\pm$1.83 | 8.88$\pm$0.54 | 0.62$\pm$0.25 | 80.11$\pm$2.44 | 7.05$\pm$1.09 | 8.91$\pm$1.43 | 8616
VACA | 3.64$\pm$1.64 | 13.82$\pm$1.87 | 2.90$\pm$0.49 | 25.82$\pm$0.33 | 23.52$\pm$1.83 | 16.02$\pm$1.32 | 5334
NLIN | MultiCVAE | 46.65$\pm$9.05 | 218.23$\pm$30.38 | 17.50$\pm$4.85 | 27.55$\pm$1.98 | 52.34$\pm$4.50 | 32.45$\pm$2.22 | 1785
CAREFL | 13.55$\pm$2.25 | 19.01$\pm$1.05 | 1.61$\pm$0.44 | 103.03$\pm$2.26 | 9.38$\pm$4.21 | 10.38$\pm$4.20 | 828
VACA | 6.04$\pm$1.87 | 9.53$\pm$3.60 | 1.44$\pm$0.59 | 17.93$\pm$0.19 | 13.24$\pm$1.93 | 8.26$\pm$1.29 | 5334
NADD | MultiCVAE | 16.99$\pm$7.22 | 133.83$\pm$16.96 | 2.59$\pm$3.38 | 18.25$\pm$0.83 | 42.43$\pm$0.87 | 20.55$\pm$1.15 | 3243
CAREFL | 13.58$\pm$1.69 | 74.51$\pm$14.18 | 1.31$\pm$0.84 | 148.03$\pm$6.96 | 11.62$\pm$1.50 | 53.71$\pm$0.98 | 9630
VACA | 1.72$\pm$0.67 | 30.10$\pm$4.41 | 0.19$\pm$0.14 | 26.64$\pm$1.30 | 22.50$\pm$2.61 | 41.16$\pm$3.91 | 5334
M-graph |  LIN | MultiCVAE | 38.60$\pm$4.06 | 80.89$\pm$20.73 | 25.82$\pm$6.26 | 17.73$\pm$3.39 | 54.72$\pm$4.71 | 27.47$\pm$2.28 | 933
CAREFL | 19.55$\pm$3.48 | 15.38$\pm$0.75 | 0.77$\pm$0.23 | 68.32$\pm$1.35 | 6.94$\pm$0.23 | 9.97$\pm$0.15 | 18200
VACA | 1.76$\pm$0.34 | 4.60$\pm$0.81 | 1.86$\pm$0.23 | 11.79$\pm$0.11 | 12.83$\pm$0.69 | 9.10$\pm$0.51 | 3249
NLIN | MultiCVAE | 37.95$\pm$6.21 | 84.79$\pm$18.15 | 24.97$\pm$9.10 | 16.66$\pm$3.28 | 51.30$\pm$5.24 | 25.90$\pm$1.23 | 933
CAREFL | 20.72$\pm$2.78 | 18.63$\pm$1.36 | 2.57$\pm$0.69 | 74.41$\pm$1.02 | 9.43$\pm$0.34 | 24.94$\pm$0.31 | 27160
VACA | 1.95$\pm$0.40 | 6.78$\pm$0.95 | 2.27$\pm$0.41 | 12.55$\pm$0.43 | 16.26$\pm$0.95 | 17.65$\pm$1.12 | 8001
NADD | MultiCVAE | 19.97$\pm$5.03 | 80.25$\pm$14.12 | 0.92$\pm$0.60 | 32.94$\pm$1.77 | 38.28$\pm$1.00 | 32.61$\pm$1.11 | 7277
CAREFL | 19.36$\pm$1.59 | 18.37$\pm$0.78 | 0.28$\pm$0.09 | 79.71$\pm$2.30 | 32.66$\pm$0.31 | 49.78$\pm$0.19 | 33950
VACA | 3.06$\pm$0.84 | 17.80$\pm$1.78 | 0.09$\pm$0.05 | 21.22$\pm$0.56 | 41.06$\pm$0.31 | 48.39$\pm$0.35 | 8001
chain |  LIN | MultiCVAE | 29.38$\pm$2.96 | 131.24$\pm$12.04 | 6.75$\pm$1.66 | 38.96$\pm$2.17 | 59.15$\pm$1.61 | 41.37$\pm$1.07 | 3099
CAREFL | 11.70$\pm$1.46 | 12.00$\pm$1.23 | 0.98$\pm$0.38 | 81.15$\pm$2.55 | 9.90$\pm$1.42 | 11.72$\pm$1.80 | 828
VACA | 4.47$\pm$0.72 | 5.73$\pm$1.11 | 3.60$\pm$0.62 | 22.91$\pm$0.23 | 23.09$\pm$1.18 | 15.77$\pm$0.77 | 5648
NLIN | MultiCVAE | 30.33$\pm$4.53 | 159.00$\pm$11.81 | 8.68$\pm$4.51 | 43.49$\pm$2.91 | 60.34$\pm$1.87 | 40.88$\pm$1.08 | 1641
CAREFL | 12.35$\pm$2.03 | 11.08$\pm$1.02 | 1.04$\pm$0.43 | 79.06$\pm$2.16 | 10.56$\pm$0.74 | 12.05$\pm$1.17 | 6462
VACA | 4.66$\pm$1.06 | 6.17$\pm$1.24 | 3.69$\pm$0.46 | 22.41$\pm$0.45 | 23.17$\pm$2.04 | 15.69$\pm$1.25 | 4445
NADD | MultiCVAE | 29.65$\pm$2.89 | 132.74$\pm$11.58 | 7.10$\pm$1.75 | 40.59$\pm$2.23 | 59.21$\pm$1.66 | 41.26$\pm$1.08 | 3099
CAREFL | 11.24$\pm$1.41 | 10.86$\pm$1.04 | 1.15$\pm$0.32 | 80.26$\pm$2.52 | 10.18$\pm$2.11 | 12.92$\pm$2.86 | 1035
VACA | 4.67$\pm$0.91 | 6.09$\pm$1.16 | 3.98$\pm$0.72 | 23.14$\pm$0.20 | 23.58$\pm$1.33 | 15.96$\pm$0.81 | 5648
loan |  | MultiCVAE | 90.38$\pm$11.31 | 213.65$\pm$5.38 | 12.24$\pm$1.33 | 65.78$\pm$1.13 | 40.98$\pm$0.35 | 15.12$\pm$0.16 | 33717
- | CAREFL | 22.10$\pm$1.64 | 27.38$\pm$4.07 | 6.74$\pm$4.25 | 50.13$\pm$2.47 | 11.15$\pm$2.57 | 6.59$\pm$0.38 | 2880
| VACA | 2.22$\pm$0.25 | 6.87$\pm$0.66 | 4.35$\pm$0.35 | 3.83$\pm$0.08 | 10.30$\pm$0.40 | 6.41$\pm$0.11 | 30402
adult |  | MultiCVAE | 140.15$\pm$6.37 | 155.52$\pm$5.93 | 12.18$\pm$2.36 | 63.52$\pm$4.05 | 39.96$\pm$0.36 | 16.37$\pm$0.65 | 6549
| CAREFL | 31.31$\pm$1.58 | 34.31$\pm$5.77 | 12.54$\pm$3.17 | 41.26$\pm$3.44 | 1.23$\pm$0.17 | 3.55$\pm$0.90 | 127420
| VACA | 4.51$\pm$0.45 | 12.68$\pm$1.95 | 1.65$\pm$0.23 | 3.37$\pm$0.09 | 5.33$\pm$0.27 | 5.67$\pm$0.20 | 63432
Table 10: Performance of different methods at estimating the observational,
interventional and counterfactual of different SCMs. All metrics are
multiplied by 100.

Figure 17: Training time (in minutes) for MultiCVAE, CAREFL and VACA for the
different graphs.

Figure 18: Number of parameters of MultiCVAE, CAREFL and VACA for the
different graphs.

In the following we present additional results that empirically show the
potential of VACA to model interventional and counterfactual queries. In
particular, we report the results for the triangle, M-graph, and chain graphs.
We remark that the following results are consistent with the ones reported in
the main manuscript for the collider, loan and adult.

#### Results for interventional distributions.

Table 10 (middle columns) reports the MMD, MeanE, and StdE for the
interventional distribution. In accordance with the results shown in the main
manuscript, we can observe that i) VACA consistently outperforms other methods
in terms of MMD; ii) the three methods provide comparable results in capturing
the mean of the interventional distribution (MeanE); and iii) CAREFL and
MultiCVAE often fail to capture the standard deviation of the interventional
distribution (StdE), while VACA provides a more accurate estimate of the
overall interventional distribution (as can be easily seen in the MMD).

#### Results for the counterfactuals.

Table 10 also reports the results for the counterfactual distribution, in
terms of MSE and SSE. As reported in the main text, we observe that CAREFL
provides more accurate estimates than VACA and MultiCVAE in terms of MSE,
which may be explained by the fact that CAREFL performs exact inference as
opposed to the approximated inference of the other two approaches. However,
CAREFL presents high variance in its results (see SSE). In contrast, VACA
leads to regularly lower values of SSE, which suggests more consistent
counterfactual estimations across factual samples and interventions.

###  E.5 Complexity analysis

In this section we compare the amount of time that it takes VACA and the two
baselines to converge during training. Figure 17 shows the time (in minutes)
it takes to train the models for the configurations of hyperparameters cross-
validated (see Section E.2). Note we also show the results averaged over 10
different initializations. We can extract three main points from Figure 17.
Firstly, we observe that most of the experiments take less than 60 minutes to
converge, only some experiments with the loan graph take longer. Secondly,
MultiCVAE takes longer to converge than CAREFL and VACA on average. This can
be explained by the fact that nodes of the graph are trained sequentially.
Thirdly, we observe that CAREFL and VACA take similar amount of time to
converge on average. We remark that experiments were run in a computer
cluster, whose performance depends on its congestion, i.e. the number of
people using the cluster and the amount of experiments. Thus, it is possible
that part of the variance in the times is due to different congestion
situations in the cluster.

Figure 18 shows a box-plot of the number of parameters of the configuration of
hyperparameters cross-validated for the models and the datasets under study.
Firstly, we observe that the configurations chosen for MultiCVAE contain less
number of parameters than the configurations of the other two methods.
However, increasing the number of parameters leads to an increase in training
time. MultiCVAE already has the longest training time compared to the other
methods, so the comparison is unfair in terms of time complexity. Regarding
CAREFL and VACA, we observe that the number of parameters cross-validated
overlaps in all the cases.

###  E.6 Computing infrastructure

All the experiments conducted in this work were executed in the same computer
cluster based on the Linux OS. Each experiment was assigned to 2 CPUs (they
could be assigned only to 1 CPU but the more CPUs, the faster the data-loading
process is) and 8GB of memory. For information about the required software
packages please refer to the official GitHub repository of VACA
https://github.com/psanch21/VACA.

##  Appendix F Further details on the counterfactual fairness use-case

In this section we provide further details on dataset, training, metrics and
additional results for the use-case of counterfactual fairness in Section 6.

###  F.1 German Credit Dataset

The German Credit dataset from the UCI repository [7] contains 20 attributes
from 1000 loan applicants. We rely on the causal model proposed by in [3] for
the following subset of features as exogenous variables ${\mathbf{X}}$ (see
Figure 19): sensitive feature $S=\\{\textit{sex}\\}$, and non-sensitive
features $C=~{}\\{\textit{age}\\}$, $R=\\{\text{{credit amount}, {repayment
history}\\}}$ and $H=\\{\textit{checking
account},\textit{savings},\textit{housing}\\}$. The causal graph in Figure 19
has a diameter ${\delta}=1$ and longest path ${\gamma}=1$. The goal of a
classifier $h$ is to predict $Y=\\{creditrisk\\}$ from ${\mathbf{X}}$. We load
and pre-process the data using the $\operatorname{aif360}$ library such that
the dataset contains binary outcome variable $Y$ (0-bad, 1-good) and a binary
sensitive attribute $S$ (0-female, 1-male). Note that the dataset contains 700
labels $Y=1$ and 300 labels $Y=0$, i.e., it is imbalanced. It also contains
690 males $S=1$ and 310 females $Y=0$. Note also that the causal model
contains heterogeneous causal nodes ($R$ and $S$), as addressed in Section
4.4. For example, [3] assume that the relationship between credit amount and
repayment history is unknown, or that it may be affected by hidden
confounders. This leads to an undirected path between the random variables and
they are grouped together in one multidimensional causal node $R$. This
applies similarly to node $S$.

$S$$R$$H$$C$ Figure 19: Causal graph for variables ${\mathbf{X}}$ of the
German Credit dataset [3]

###  F.2 Training

In this section we provide further information on training VACA on the German
Credit dataset [7] and detail the different classifiers in Section 6. We use a
80% training, 10% validation, 10% training data split.

#### VACA.

Training for VACA was performed on normalized data—performing normalization
only on the continuous variables, i.e. r.v. $C$ and $R$ in Fig 19. We trained
a heterogeneous VACA as described in Section with a message passing function
$f^{m}$ with one hidden layer of $16$ neurons, a decoder with one hidden layer
of $16$ neurons and a latent variable with dimension $4$. We trained the model
using the PIWAE [35] approach with $\lambda_{KLD}=0.05$, specifically, the
encoder with the IWAE [2] objective with $K=5$ and the decoder with a
$\beta$-ELBO with $\beta=0.5$. We use a parents dropout rate (see Appendix
E.2) of $0.2$, learning rate of $0.005$ and batch size $100$.

#### Classifiers.

Logistic Regresion ($\operatorname{LR}$) and Support Vector Machine
($\operatorname{SVM}$) classifiers are taken from the scikit-learn library and
trained with default parameters as well as
$\operatorname{class\\_weight}=\operatorname{balanced}$ due to the class
imbalance of the dataset.

###  F.3 Metrics

In this section we detail the measures f1-score ($\operatorname{f1}$),
unfairness ($\operatorname{uf}$) and accuracy ($\operatorname{acc}$) in Table
11.

#### f1-score.

Due to class imbalance, we measure classifier performance with the f1-score.
The f1-score is the weighted average of the precision and recall and can
assume values between 0 and 1; the higher the values the better. Our
implementation relies on the $\operatorname{f1\\_score}$ from the scikit-learn
library. We compute in expectation over our training dataset:

$\displaystyle\operatorname{f1}={{\mathbb{E}}}\left[2\times\frac{\operatorname{precision}\times\operatorname{recall}}{\operatorname{precision}+\operatorname{recall}}\right].$
(14)

where $\operatorname{precision}$ is the ratio $\frac{TP}{TP+FP}$ with the
number of true positives $TP$ and the number of false positives $FP$ and
$\operatorname{recall}$ is the ratio $\frac{TP}{TP+FN}$ with the number of
false positives $FP$.

#### Counterfactual (un)fairness.

We measure counterfactual unfairness [23] with counterfactual instances
$\mathbf{x}^{CF}$ and classifier prediction
$h(\mathbf{x}^{CF})={\hat{y}^{CF}}$ as expectation over our training dataset:

$\displaystyle\operatorname{uf}={{\mathbb{E}}}\left[\left|p({y^{F}})=1|\mathbf{x}^{F})-p({\hat{y}^{CF}}=1|do(S=a^{\prime}),\mathbf{x}^{F})\right|\right]$
(15)

where $a^{\prime}=1-a$ as $S\in\\{0,1\\}$.

#### Accuracy.

In Table 11, we report additionally the prediction accuracy as performance
measure of classifier $h$ with respect to factuals (samples)
$(\mathbf{x}^{F},{y^{F}})$ and prediction $h(\mathbf{x}^{F})={\hat{y}^{F}}$ in
expectation over our training dataset:

$\displaystyle\operatorname{acc}={{\mathbb{E}}}\left[\mathds{1}\left({y^{F}}_{i}={\hat{y}^{F}}_{i}\right)\right].$
(16)

Our implementation relies the $\operatorname{accuracy\\_score}$ from the
scikit-learn library.

Metr. | CLF | full | unaware | fair-x | fair-z
---|---|---|---|---|---
$\uparrow$ acc | LR | 61.00 | 62.00 | 43.00 | 63.60 $\pm$ 4.74
| SVM | 66.00 | 64.00 | 51.00 | 64.20 $\pm$ 4.52
$\uparrow$ f1 | LR | 67.23 | 68.33 | 47.71 | 70.89 $\pm$ 4.11
| SVM | 71.67 | 69.49 | 59.50 | 70.79 $\pm$ 5.15
$\downarrow$ uf | LR | 18.30 $\pm$ 3.22 | 17.65 $\pm$ 3.20 | 0.16 $\pm$ 0.02 | 0.44 $\pm$ 0.24
| SVM | 14.01 $\pm$ 2.26 | 13.27 $\pm$ 2.28 | 0.14 $\pm$ 0.02 | 0.51 $\pm$ 0.19
Table 11: Counterfactual unfairness (uf), accuracy (acc) and f1-score (f1) of
a LR and SVM classifier over 10 VACA seeds. Values multiplied by 100.