As with the GAN-NCM, for estimation experiments, only one model is trained,
the query loss $\mathbb{D}_{Q}$ is not computed, and optimization is purely to
fit the datasets in $\mathbb{Z}$. For each setting in Fig. 6, we ran 5
different settings of $n$, increasing exponentially from $10^{3}$ to $10^{5}$,
but only for $d=1$ due to tractability issues with higher dimensions. We ran
10 trials for each setting, with means and confidence intervals computed
across these 10 trials.

We use the full batch per epoch. We use a learning rate of $\eta=4\cdot
10^{-3}$, and run for up to 3000 epochs with early stopping. We use $m=10^{6}$
Monte Carlo samples to evaluate the query from $\widehat{M}$ via the Gumbel-
max trick.

###  B.7 Runtime Experiment Procedure

For the runtime experiments procedure, we conducted NCM-GAN and NCM-MLE on the
same dataset generated in the first graph of Fig. 4 as the dimensionality $d$
of $Z$ scales from 1 to 7. Similar to estimation experiments, only one model
is trained, the query loss is not computed, and optimization is purely to fit
the datasets. We record the first 100 epochs of training time for both
methods. For each setting in Fig. 5, positivity is enforced when generating
data and $n=10000$. We ran 4 trials for each setting, with means and
confidence intervals computed across these 4 trials.

For GAN-NCM, we use a learning rate of $\eta=2\cdot 10^{-5}$ and a batch size
of 1000. For MLE-NCM, we use a learning rate of $\eta=4\cdot 10^{-3}$ and the
full batch per epoch.

###  B.8 Hardware

The MLE-NCM was trained on NVIDIA Tesla V100 GPUs provided by Amazon Web
Services. All other experiments are run via jobs on a cluster, and each thread
may be assigned either an NVIDIA K80 or NVIDIA P100 GPU. In total,
approximately 10000 GPU hours were used for experiments.

##  Appendix C Connection to Related Works

This paper evaluates causal queries directly through a parameterized model in
the form of an SCM trained on available data from the true SCM
$\mathcal{M}^{*}$, which we call a proxy model-based approach. In this
section, we discuss the differences between the NCM presented in this paper
and the models used in other works that also implement a proxy approach.

First, we note that many of these works implicitly benefited from the
constraints of enforcing a causal diagram into the model structure without
proving any properties about these constraints. One contribution of our work
is to provide some explicit results about these constraints via Thms. 1 and 2.
Second, and more importantly, most of these works operated under different
assumptions and problem settings, to be discussed in individual sections. When
applied to our problem setting, these methods may not be as general as the
method discussed in this paper. Some of these losses of generality are
highlighted in the following propositions and illustrated through the further
examples provided in App. D.

######  Proposition 1.

Any model in which the Markovianity assumption (also known as causal
sufficiency) is enforced is not $\mathcal{L}_{3}$-$\mathcal{G}$ expressive
(i.e. an equivalent statement of Thm. 2 cannot be proven about such a model
class). $\blacksquare$

###### Proof.

See counterexample provided in Example 6 in Sec. D.3. ∎

######  Proposition 2.

There exist settings in which an $\mathcal{L}_{3}$ query that is not
identifiable from $\mathcal{L}_{1}$ data ($P(\mathbf{V})$) alone may be
identifiable given data sets from $\mathcal{L}_{2}$. $\blacksquare$

###### Proof.

Examples include ID settings 9, 10, 11, and 12 from Fig. 4, in which their
counterparts without the $P(\mathbf{V}_{x})$ datasets are not identifiable.
Another example is shown in Example 9 in Sec. D.3. ∎

######  Proposition 3.

There exist settings in which an $\mathcal{L}_{3}$ query is not identifiable
even in the Markovian case and given both $\mathcal{G}$ and all datasets from
$\mathcal{L}_{2}$. $\blacksquare$

###### Proof.

Example is shown in Example 8 from Sec. D.3. Further, this is an example of a
case where the query cannot even be bounded. ∎

We discuss more specifically the differences with each paper in the following
subsections. Whenever notation differs between the two works, we use our
notation to maintain consistency and make the comparison more transparent.

###  C.1 Kocaoglu et al. (2018)

Similar to our work, Kocaoglu et al. (2018) used an implicit approach to fit a
proxy SCM on the data via generative adversarial networks. The method was
specifically trained on observational ($\mathcal{L}_{1})$ data and used to
evaluate interventional ($\mathcal{L}_{2}$) queries in the Markovian setting.
Like for a $\mathcal{G}$-NCM, their approach used a proxy model
$\widehat{M}=\langle\widehat{\mathbf{U}},\mathbf{V},\widehat{\mathcal{F}},P(\widehat{\mathbf{U}})\rangle$,
where $\widehat{\mathbf{U}}$ consists of independent fixed noise variables and
each $\hat{f}_{V_{i}}\in\widehat{\mathcal{F}}$ is defined such that
$\hat{f}_{V_{i}}$ takes $\mathbf{Pa}_{V_{i}}$ and
$\widehat{\mathbf{U}}_{V_{i}}$ as input. Here, $\mathbf{Pa}_{V_{i}}$ is the
set of parents of $V_{i}$ in the graph $\mathcal{G}$, and
$\widehat{\mathbf{U}}_{V_{i}}$ is a set of noise variables specific to
variable $V_{i}$.

This work introduced two novel GAN architectures for $\mathcal{L}_{2}$
inference, CausalGAN and CausalBEGAN, focused on applications in image
generation. For this challenging task, the CausalGAN architecture was designed
in a unique fashion in which the model was separated into two parts: one for
generating the labels, and another for generating the image conditioned on
these labels. This was in contrast to having one generator learning the entire
joint distribution of labels and images. The approach for learning the image
generator is similar to a conditional GAN in that the objective of the
generator includes an additional term to maximize the prediction accuracy of a
labeling network. However, this often results in a problem they described as
_label-conditioned mode collapse_ , where the generator may have been
complacent in generating a small number of fixed samples that are easy to
label. Hence, in addition to this labeler, CausalGAN also introduced an anti-
labeler network. The anti-labeler also attempts to label the generated images,
but the objective of CausalGAN simultaneously maximizes the accuracy of the
labeler while minimizing the accuracy of the anti-labeler. This prevents
label-conditioned mode collapse because images that are too easily labeled
will fail to minimize the accuracy of the anti-labeler. More formally, for a
binary label $l$ with probability $\rho=P(l=1)$, the objective of the
CausalGAN generator is:

$\displaystyle\min_{G}\mathbb{E}_{x\sim
p_{g}(x)}\left[\log\left(\frac{1-D(x)}{D(x)}\right)\right]$
$\displaystyle-\rho\mathbb{E}_{x\sim
p_{g}^{1}(x)}[\log(D_{LR}(x)]-(1-\rho)\mathbb{E}_{x\sim
p_{g}^{0}(x)}[log(1-D_{LR}(x))]$ $\displaystyle+\rho\mathbb{E}_{x\sim
p_{g}^{1}(x)}[\log(D_{LG}(x)]+(1-\rho)\mathbb{E}_{x\sim
p_{g}^{0}(x)}[log(1-D_{LG}(x))],$ (17)

where $G$ is the generator, $D$ is the discriminator, $p_{g}(x)$ is the
distribution induced by the generator, $p_{g}^{i}(x)$ is the distribution
conditioned on $l=i$, $D_{LR}$ is the labeler, and $D_{LG}$ is the anti-
labeler. The loss is composed of the loss from the discriminator subtracted
with the loss from the labeler and adding the loss from the anti-labeler. The
paper proved that the generator successfully matched the true conditional
image distribution when $G$ was the optimum of this objective (Kocaoglu et
al., 2018, Thm. 2). CausalBEGAN was another architecture introduced where
optimization considered image quality margin and label margin. In both cases,
interventions were computed using the mutilation procedure on the generating
model $\widehat{M}$.

CausalGAN and CausalBEGAN were empirically shown to successfully learn the
data distribution in several image generation tasks. Experiments demonstrated
several settings in which intervening on a label resulted in a different image
compared to conditioning on the label.

The reference produced impressive results in generating images under causal
interventions using novel implicit GAN architectures. While the approach is
similar to our work in that a GAN proxy model is used, there are a number of
differences with our work.

  1. 1.

The paper worked in the space of inferring $\mathcal{L}_{2}$ quantities given
observational data from $\mathcal{L}_{1}$. Our work generalizes this setting
to inferring any counterfactual query from $\mathcal{L}_{3}$ given any
combination of datasets from $\mathcal{L}_{1}$ and $\mathcal{L}_{2}$.

  2. 2.

Unlike this work, we do not make the Markovianity assumption, meaning we allow
unobserved confounding between variables. Many real world settings often have
unobserved confounding, so the Markovianity assumption can be limiting. In a
situation where the true SCM has unobserved confounding, a model class which
enforces Markovianity will not be expressive enough to model the true SCM
(Prop. 1).

  3. 3.

By (Kocaoglu et al., 2018, Def. 1), the architecture of CausalGAN took the
form of a causal diagram $\mathcal{G}$, similar to how the
$\mathcal{G}$-constrained NCM is constructed in Def. 2. In our work, we
discuss the theoretical properties of the implied constraints via Thms. 1 and
2. While $\mathcal{G}^{(\mathcal{L}_{2})}$-consistency likely held for
CausalGAN due to the nature of its construction, this property was not
explicitly proven or discussed.

  4. 4.

The fact that all $\mathcal{L}_{2}$ quantities are identifiable from
$\mathcal{L}_{1}$ data was acknowledged (Kocaoglu et al., 2018, Prop. 1). The
identification question constitutes an additional challenge in our paper since
we consider layer 3 queries and allow for unobserved confounding. It’s
understood that without Markovianity, even inferences about $\mathcal{L}_{2}$
queries are not always identifiable from $\mathcal{L}_{1}$ data.

###  C.2 Pawlowski et al. (2020)

Pawlowski et al. (2020) implemented a deep proxy SCM approach in which each
mechanism $\hat{f}_{V_{i}}$ models a conditional distribution
$P(V_{i}\mid\mathbf{Pa}_{V_{i}})$ under the Markovianity assumption (lack of
unobserved confounders). Under this assumption, the joint observational
distribution $P(\mathbf{V})$ can be factorized as

$P(\mathbf{V})=\prod_{V_{i}\in\mathbf{V}}P(V_{i}\mid\mathbf{Pa}_{V_{i}}),$
(18)

so each mechanism represents one conditional distribution in this
factorization. They proposed three novel approaches to modeling this
conditional distribution.

  1. 1.

The first approach involved modeling $P(X\mid\mathbf{Pa}_{X})$ using
normalizing flows. $\widehat{\mathbf{U}}_{X}$ consists of a noise variable
$\epsilon$, and $\hat{f}_{X}$ is an invertible function such that
$x=\hat{f}_{X}(\epsilon;\mathbf{pa}_{X})$. Then the explicit likelihood could
be computed as
$p(x\mid\mathbf{pa}_{x})=p(\epsilon)\cdot|\det\nabla_{\epsilon}\hat{f}_{X}(\epsilon;\mathbf{pa}_{x})|^{-1}$,
where $\epsilon$ could be computed as $\hat{f}_{X}^{-1}(x;\mathbf{pa}_{x})$.

  2. 2.

The second approach extended the first approach to high dimensional settings
by splitting $\epsilon$ into two independent noise variables, $z$ and $u$, and
splitting $\hat{f}_{X}$ into two components: a non-invertible function $g$ and
an invertible function $h$. Then
$\hat{f}_{X}(\epsilon,\mathbf{pa}_{x})=h(u;g(z;\mathbf{pa}_{x}),\mathbf{pa}_{x})$,
where $g$ is expected to capture the high-level aspects of the data. This is
helpful in high dimensional settings as the size of the invertible component
$h$ can be reduced because working with the outputs of $g$ is easier than
directly working with the data. Since $p(x\mid\mathbf{pa}_{x})$ is no longer
tractable, the model is instead trained using variational inference by
maximizing the evidence lower bound (ELBO) with variational distribution
$Q(z\mid x,\mathbf{pa}_{x})$,

$\log p(x\mid\mathbf{pa}_{x})\geq\mathbb{E}_{Q}(z\mid x,\mathbf{pa}_{x})[\log
p(x\mid z,\mathbf{pa}_{x})]-D_{\operatorname{KL}}[Q(z\mid
x,\mathbf{pa}_{x})||P(z)].$ (19)

  3. 3.

Although not demonstrated in the paper, a third approach was presented
suggesting the use of a GAN-based approach to implicitly model the likelihood.

The method primarily used in the paper was the 2nd approach. They evaluated
counterfactual quantities using the 3-step procedure from Pearl (2000, Thm.
7.1.7):

  1. 1.

Abduction – The observations $\mathbf{v}$ of all variables $\mathbf{V}$ are
used to predict $P(\mathbf{U}\mid\mathbf{v})$, where
$P(\mathbf{U}\mid\mathbf{v})$ decomposes as $\prod_{V_{i}\in
V}P(\epsilon_{V_{i}}\mid v_{i},\mathbf{pa}_{V_{i}})$ due to the Markovianity
assumption. This is then computed as

$\displaystyle P(\epsilon_{V_{i}}\mid v_{i},\mathbf{pa}_{V_{i}})$
$\displaystyle=P(z_{V_{i}}\mid v_{i},\mathbf{pa}_{V_{i}})P(u_{V_{i}}\mid
z_{V_{i}},v_{i},\mathbf{pa}_{V_{i}})$ $\displaystyle\approx Q(z_{V_{i}}\mid
e_{V_{i}}(v_{i};\mathbf{pa}_{V_{i}}))\delta(u_{V_{i}})$ (20)

where $e(v_{i};\mathbf{pa}_{V_{i}})$ is a deterministic encoder mapping
$v_{i}$ and $\mathbf{pa}_{V_{i}}$ to $z_{V_{i}}$, and $\delta$ is a Dirac
delta distribution centered at
$h_{V_{i}}^{-1}(v_{i};g_{V_{i}}(z_{V_{i}};\mathbf{pa}_{V_{i}}),\mathbf{pa}_{V_{i}})$.

  2. 2.

Action – For an intervention over sets of variables $\mathbf{X}$, each
mechanism $\hat{f}_{X}$ for $X\in\mathbf{X}$ is replaced either with a fixed
value $x$ (atomic) or a surrogate mechanism
$\tilde{f}_{X}(\epsilon_{X}\mid\widetilde{\mathbf{pa}}_{X})$ (stochastic).

  3. 3.

Prediction – Under this distribution of $P(\mathbf{U}\mid\mathbf{v})$, and
with the modified mechanisms, the deep SCM is evaluated to produce a result
for the counterfactual quantity $P(\mathbf{V}_{\mathbf{x}}\mid\mathbf{v})$.

The approach was demonstrated empirically on a synthetic MNIST dataset and on
a real world brain imaging dataset.

The paper presented a powerful approach to modeling a deep SCM using
normalizing flows and variational inference, allowing for efficient and
precise training even in high-dimensional cases. Still, there are a number of
differences between this paper and our work.

  1. 1.

Unlike this work, we do not make the Markovianity assumption, meaning we allow
unobserved confounding between variables. Many real world settings often have
unobserved confounding, so the Markovianity assumption can be limiting. In a
situation where the true SCM has unobserved confounding, a model class which
enforces Markovianity will not be expressive enough to model the true SCM
(Prop. 1). Moreover, it is not a trivial task to extend the normalizing flow
approach in this paper to a case with unobserved confounders. In such a case,
Eq. 18 would no longer be valid, and it is not immediately clear how to train
invertible mechanisms that can share noise.

  2. 2.

Like our work, this paper assumed that the causal graph $\mathcal{G}$ is
provided. However, in our work, we discuss the theoretical properties of the
implied constraints via Thms. 1 and 2, while this paper implicitly benefited
from these constraints without discussing these properties.

  3. 3.

The identification problem is not considered in this work while it is a major
focus in our work. As shown in Prop. 3, there exist counterfactual quantities
that cannot be identified given the graph $\mathcal{G}$, and all available
data from $\mathcal{L}_{1}$ and $\mathcal{L}_{2}$, even under the Markovianity
assumption. If these quantities were to be estimated from a model trained on
such data, the result will likely be incorrect and/or misleading. For that
reason, we make sure to identify all queries before attempting to estimate
them.

  4. 4.

While this paper only considered input in the form of observational data from
$\mathcal{L}_{1}$, our work considers cases where additional experimental
datasets from $\mathcal{L}_{2}$ are provided (even cases where perhaps the
observational distribution $P(\mathbf{V})$ is not provided). As shown in Prop.
2, there do exist counterfactual queries that are not identifiable given
observational data but are identifiable given datasets from $\mathcal{L}_{2}$.

  5. 5.

This paper evaluated counterfactual quantities using the 3-step procedure of
abduction, action, and prediction. While this is a sound algorithm for
evaluating applicable quantities, it does not apply generally to all
counterfactual quantities. Specifically, it only applies in quantities of the
form $P(\mathbf{Y}_{\mathbf{x}}\mid\mathbf{z})$, in which only one
intervention is applied before the conditioning bar, and observations are
present after the conditioning bar. This excludes quantities such as
$P(y_{x},y^{\prime}_{x^{\prime}})$, which is a joint probability of two
counterfactual terms in different worlds (one in $X=x$ and another in
$X=x^{\prime}$), as well as quantities such as $P(y_{x}\mid
y^{\prime}_{x^{\prime}})$, in which a term in an intervened world
$(X=x^{\prime})$ is after the conditioning bar.

  6. 6.

Finally, although an implicit approach via GANs was discussed in (Pawlowski et
al., 2020), it was not explored. Our paper realizes an implementation of GAN
approach in the non-Markovian setting.

###  C.3 Witty et al. (2021)

Witty et al. (2021) introduced simulation-based identifiability (SBI), which
uses a Bayesian approach to solve the identification problem for
$\mathcal{L}_{2}$ and some $\mathcal{L}_{3}$ queries given $\mathcal{L}_{1}$
data. In their approach, the goal was to identify a particular query $Q$ given
a fixed observational distribution $P(\mathbf{V})$ and a causal diagram
$\mathcal{G}$. An SCM
$\mathcal{M}^{*}=\langle\mathbf{U}^{*},\mathbf{V},\mathcal{F}^{*},P(\mathbf{U}^{*})\rangle$
is sampled from a prior over SCMs (specifically their functions and exogenous
variables), $P(\mathcal{F},\mathbf{U})$. In each iteration, $n$ samples are
drawn from $P^{\mathcal{M}^{*}}(\mathbf{V})$, the observational dataset of
$\mathcal{M}^{*}$. Similar to our work, two parameterized SCMs,
$\widehat{M}_{1}$ and $\widehat{M}_{2}$ are optimized to simultaneously fit
$P^{\mathcal{M}^{*}}(\mathbf{V})$ while maximizing the difference in the
induced queries of the two SCMs, $\widehat{Q}_{1}$ and $\widehat{Q}_{2}$. This
is done by maximizing the objective

$L(\widehat{M}_{1},\widehat{M}_{2},\\{\mathbf{v}_{i}\\}_{i=1}^{n})=\sum_{i=1}^{n}\left(\log
P^{\widehat{M}_{1}}(\mathbf{v}_{i})+\log
P^{\widehat{M}_{2}}(\mathbf{v}_{i})\right)+\lambda|\widehat{Q}_{1}-\widehat{Q}_{2}|,$
(21)

where $\\{\mathbf{v}_{i}\\}_{i=1}^{n}$ are a collection of $n$ samples from
$P^{\mathcal{M}^{*}}(\mathbf{V})$. Maximizing this objective is equivalent to
jointly maximizing the likelihood of the data for both $\widehat{M}_{1}$ and
$\widehat{M}_{2}$ as well as maximizing the gap between $\widehat{Q}_{1}$ and
$\widehat{Q}_{2}$. Identifiability is concluded if, following a hypothesis
testing procedure, $|\widehat{Q}_{1}-\widehat{Q}_{2}|$ is under a certain
threshold. Non-identifiability is concluded otherwise. They proved that $Q$ is
identifiable w.r.t. $\mathcal{M}^{*}$ if and only if there does not exist any
parameterized SCM $\widehat{M}$ fitted to the data from
$P^{\mathcal{M}^{*}}(\mathbf{V})$ such that $P^{\widehat{M}}(\mathbf{V})$
converges to $P^{\mathcal{M}^{*}}(\mathbf{V})$ as $n\rightarrow\infty$,
$Q^{\widehat{M}}\neq Q^{\mathcal{M}^{*}}$, and $\widehat{M}$ is in the support
of the prior (Witty et al., 2021, Thm. 3.2).

In non-parametric settings, there exist sound and complete symbolic methods
such as do-calculus that could efficiently and correctly identify any
$\mathcal{L}_{2}$ query Pearl (2000). Hence, as we understand, the paper
primarily focused its attention on parametric settings, where the model class
of SCMs was further restricted to specific function classes (e.g. linear), and
symbolic approaches to identification under this constrained class were not as
well-studied.

For the identification task, SBI shares many similarities with our approach in
Alg. 1, notably that they also trained two proxy SCMs to minimize/maximize the
query while simultaneously fitting the data. Still, there are many differences
between the two works.

  1. 1.

The paper considered a general form of the query
$Q(\mathbf{V}_{\mathbf{X}=\mathbf{x}},\mathbf{V}_{\mathbf{X}=\mathbf{x}^{\prime}})$,
which was a function of the counterfactual terms of $\mathbf{V}$ in the two
worlds of $\mathbf{X}=\mathbf{x}$ and $\mathbf{X}=\mathbf{x}^{\prime}$. For
example, average treatment effect (ATE) given a binary treatment $X$ could be
computed as
$Q(\mathbf{V}_{\mathbf{X}=\mathbf{x}},\mathbf{V}_{\mathbf{X}=\mathbf{x}^{\prime}})=\mathbb{E}[Y_{X=1}]-\mathbb{E}[Y_{X=0}]$,
which was what was used in their experiments. While in theory, this may have
included certain counterfactual queries, these queries were not explored, so
the work largely focused on the space of $\mathcal{L}_{2}$ queries. Moreover,
they only considered data from $\mathcal{L}_{1}$. Our work generalizes this
setting to inferring any counterfactual query from $\mathcal{L}_{3}$ given any
combination of datasets from $\mathcal{L}_{1}$ and $\mathcal{L}_{2}$.

  2. 2.

Like our work, this paper assumed that the causal graph $\mathcal{G}$ is
provided. However, in our work, we discuss the theoretical properties of the
implied constraints via Thms. 1 and 2, while this paper implicitly benefited
from these constraints without discussing these properties.

  3. 3.

The paper was largely concerned with identifying quantities in parametric
cases (and functions of SCMs were fixed to certain parametric families), while
our work is largely focused on the non-parametric case. We make no claims
about expressivity or $\mathcal{G}$-consistency in parametric settings.

  4. 4.

Our work assumes that there exists a true SCM $\mathcal{M}^{*}$, which
generates a specific dataset from which we aim to identify and estimate our
query. In contrast, (Witty et al., 2021) specifically aimed to identify a
query without a true SCM in mind and instead matched the distribution with an
SCM sampled from a prior distribution over the space of SCMs. While estimation
could theoretically be performed after concluding identifiability with SBI,
this was not a focus of (Witty et al., 2021), and the estimation would not be
meaningful unless $\mathcal{M}^{*}$ corresponded to some true SCM in a real
application.

  5. 5.

In terms of the approach, (Witty et al., 2021) differs from our work in that
the two parametrized models, $\widehat{M}_{1}$ and $\widehat{M}_{2}$ are
optimized jointly, while in our case, the two models are optimized separately.
The goal in SBI is to maximize the gap between the two models in the query,
while Alg. 1 does this implicitly when maximizing or minimizing the query in
the two models separately. Further, our implementation in practice utilizes
GANs, which is not explored in this paper.

###  C.4 Xia et al. (2021)

Xia et al. (2021) introduced the NCM, which is the model class considered in
our paper. They proved that NCMs are expressive on all three layers (Xia et
al., 2021, Thm. 1) but then subsequently proved the Neural Causal Hierarchy
Theorem (Xia et al., 2021, Corol. 1) which states that despite this
expressivity, an NCM trained only on lower layer data would almost never match
the true SCM on higher layers. To proceed, they introduced a graphical
inductive bias in the form of a causal diagram $\mathcal{G}$ derived from the
true SCM. They defined a $\mathcal{G}$-constrained version of the NCM (Def. 2
in our paper).

They proved that the class of $\mathcal{G}$-constrained NCMs is both
$\mathcal{G}$-consistent (Xia et al., 2021, Thm. 2) and
$\mathcal{L}_{2}$-$\mathcal{G}$ expressive (Xia et al., 2021, Thm. 3).
$\mathcal{G}$-consistency ($\mathcal{G}^{(\mathcal{L}_{2})}$-consistency in
the context of our work) constrained the class to a smaller set of models
which may have potentially agreed on the causal query, while
$\mathcal{L}_{2}$-$\mathcal{G}$ expressiveness guaranteed that
$\mathcal{G}$-NCMs could still model any SCM within this constrained space.
Still, these results only considered up to layer 2 of the causal hierarchy,
and made no claims about layer 3. In contrast, Thms. 1 and 2 in our work are
much stronger claims, proving these results on layer 3. That is,
$\mathcal{G}$-NCMs have all of the constraints implied on the third layer
while still being able to express any true SCM on all three layers within this
constrained setting.

The paper then defined the neural identification problem specifically for
$\mathcal{L}_{2}$ queries. They proved a duality stating that neural
identification for $\mathcal{L}_{2}$ quantities is equivalent to graphical
identification (Xia et al., 2021, Thm. 4) and showed that the mutilation
procedure applied to a proxy $\mathcal{G}$-NCM can be used to estimate
identifiable queries (Xia et al., 2021, Corol. 2). This definition of neural
identification did not consider $\mathcal{L}_{3}$ quantities or
$\mathcal{L}_{2}$ input data. In our case, we define neural counterfactual
identification (Def. 4) generalizing the neural identification problem to
handle $\mathcal{L}_{3}$ queries given arbitrary datasets from
$\mathcal{L}_{1}$ and $\mathcal{L}_{2}$. We prove a stronger duality in this
setting (Thm. 3) and that counterfactual quantities can be evaluated directly
from the proxy $\mathcal{G}$-NCM, akin to the mutilation procedure for
$\mathcal{L}_{2}$ quantities (Corol. 1).

Following these identification results, (Xia et al., 2021) introduced an
algorithm (Xia et al., 2021, Alg. 1) which is both sound and complete for
identifying and estimating $\mathcal{L}_{2}$ queries (Xia et al., 2021, Corol.
4). The algorithm was implemented in practice through a maximum likelihood
approach (MLE-NCM) through (Xia et al., 2021, Alg. 2). This involved training
two $\mathcal{G}$-NCM parameterizations, where one was trained to maximize the
query, while the other was trained to minimize it, and both models were
trained to fit the observational data $P(\mathbf{V})$. The algorithm in our
counterfactual setting (Alg. 1) is similar on the surface, but there are many
additional challenges in solving this more general problem. First, while an
interventional query might be easy to compute and optimize, it is not as clear
how to compute a counterfactual quantity with gradients such that it can be
maximized or minimized. After all, counterfactual quantities often have
complicated forms such as nesting (e.g. $P(Y_{Z_{X}})$) or having the same
variable in multiple worlds (e.g. $P(Y_{X=0},Y_{X=1})$). Second, fitting a
single dataset is standard practice, which is what needs to be done when only
observational data is available. It is not as clear how to simultaneously fit
several datasets. Our work solves both of these problems in practice using the
GAN-NCM approach in Sec. 4, but also explains how one might solve this problem
using the MLE-NCM in Sec. B.6.1.

To demonstrate the ability of the MLE-NCM to identify and estimate
$\mathcal{L}_{2}$ queries in practice, the paper showed results from
identifying average treatment effect (ATE) from 8 typical graphs encountered
in the literature. They also provided estimation results for the identifiable
cases. While our experiments mainly focused on identifying $\mathcal{L}_{3}$
queries, we reproduced the experiments of (Xia et al., 2021) in Sec. B.1 using
the GAN-NCM to demonstrate that our approach is general and can be applied to
the $\mathcal{L}_{2}$ case as well. Further, although the reason why one might
prefer to use an NCM for causal inference is due to the ability of neural
models to scale inferences to high dimensions, this ability was not emphasized
in the experimental results of this paper. We explicitly showed how the GAN-
NCM can perform nearly as well even in 16-dimensional cases while showing that
the MLE-NCM scales poorly due to tractability issues (see Fig. 5 in Sec. 5).
One additional issue of the MLE-NCM is its bias towards returning ID results
even in non-ID cases. This is evident in Fig. 4 and can be potentially
problematic if used in practice since attempting to estimate a non-ID query
will likely produce an incorrect and misleading result. For more details, see
the discussion surrounding Fig. 14 in Sec. B.1.

To summarize, this paper introduced the NCM and proved several properties
about it, including results on expressiveness, constraints, and
identification. Our work also uses the NCM and proves many analogous results,
but there are still several differences:

  1. 1.

The paper primarily worked in the space of inferring $\mathcal{L}_{2}$
quantities given observational data from $\mathcal{L}_{1}$. Our work
generalizes this setting to inferring any counterfactual query from
$\mathcal{L}_{3}$ given any combination of datasets from $\mathcal{L}_{1}$ and
$\mathcal{L}_{2}$. This is reflected in the theoretical results, the
algorithms, and the practical implementation. We show that the
$\mathcal{G}$-NCM has all $\mathcal{L}_{3}$ constraints (Thm. 1) but can still
express any model on all three layers in this constrained setting (Thm. 2). We
show that the $\mathcal{G}$-NCM can be used to solve the identification and
estimation problem in the more general counterfactual setting (Thm. 3, Corol.
1, Alg. 1, Corol. 2). These results drastically generalize the capabilities of
the NCM to $\mathcal{L}_{3}$.

  2. 2.

While the MLE-NCM implemented in this paper is an option for counterfactual
identification and estimation (as shown in Sec. B.6.1), we develop a new type
of architecture for the NCM using a generative adversarial approach. Our
approach, detailed in Sec. 4, has many advantages over the MLE-NCM, notably
its design simplicity, its ability to scale to high dimensions, and its
robustness in solving the ID problem. We even show that the GAN-NCM can solve
all of the experimental settings originally presented in (Xia et al., 2021) in
Sec. B.1.

From all of these results, both theoretical and empirical, it is evident that
our approach, in fact, subsumes the work in (Xia et al., 2021).

###  C.5 Zecevic et al. (2021)

Zecevic et al. (2021) related the class of graph neural networks (GNNs) to the
class of Markovian SCMs. Given a DAG $\mathcal{G}$ with variables
$\mathbf{V}$, they defined GNN layers using the functions

$h_{V_{i}}=\phi\left(d_{V_{i}},\bigoplus_{V_{j}\in\mathcal{N}(V_{i})}\psi(d_{V_{i}},d_{V_{j}})\right),$
(22)

where $d_{V_{i}}$ is either data associated with $V_{i}$ or a feature computed
from a previous layer, $\mathcal{N}(V_{i})$ is the set of neighbors of $V_{i}$
in $\mathcal{G}$, $\bigoplus$ refers to an aggregating function such as sum,
and $h_{V_{i}}$ is a feature of $V_{i}$ for the current layer. They proved
(Zecevic et al., 2021, Thm. 1), which states that for any Markovian SCM
$\mathcal{M}=\langle\mathbf{U},\mathbf{V},\mathcal{F},P(\mathbf{U})\rangle$,
there exists a GNN construction such that $h_{V_{i}}=f_{V_{i}}$ for all
$f_{V_{i}}\in\mathcal{F}$, and such a construction (but without shared $\phi$
and $\psi$) is therefore a special case of the NCM introduced in (Xia et al.,
2021) ((Zecevic et al., 2021, Corol. 1)). Using a GNN interpretation of an
SCM, one could compute quantities from an intervention $\mathbf{X}=\mathbf{x}$
through the mutilation approach by evaluating Eq. 22 but with edges into
variables of $\mathbf{X}$ removed from $\mathcal{N}$.

Leveraging results from variational graph autoencoders (VGAEs), (Zecevic et
al., 2021) demonstrated an approach to learning the observational distribution
with a GNN architecture. They introduced interventional VGAEs (iVGAEs) which
defined inferences from VGAEs on layer 2 quantities via the mutilation
procedure on the GNN layers. They proved a wealth of results related to this
class of models following similar results proven in (Xia et al., 2021):
$\mathcal{L}_{2}$-expressivity (Zecevic et al., 2021, Thm. 3), a causal
hierarchy theorem result between layers 1 and 2 (Zecevic et al., 2021, Corol.
2), and equivalence between neural-identification from iVGAEs and classical
identification of $\mathcal{L}_{2}$ quantities (Zecevic et al., 2021, Thm. 4).

(Zecevic et al., 2021) introduced an interesting discussion on the
relationship between GNNs, SCMs, and NCMs. The work differs from our work in a
variety of ways.

  1. 1.

Like (Xia et al., 2021), (Zecevic et al., 2021) primarily worked in the space
of inferring $\mathcal{L}_{2}$ quantities given observational data from
$\mathcal{L}_{1}$. Our work generalizes this setting to inferring any
counterfactual query from $\mathcal{L}_{3}$ given any combination of datasets
from $\mathcal{L}_{1}$ and $\mathcal{L}_{2}$.

  2. 2.

Unlike this work, we do not make the Markovianity assumption, meaning we allow
unobserved confounding between variables. Many real world settings often have
unobserved confounding, so the Markovianity assumption can be limiting. In a
situation where the true SCM has unobserved confounding, a model class which
enforces Markovianity will not be expressive enough to model the true SCM
(Prop. 1). The results proven in this paper do not extend trivially to cases
with unobserved confounders. Traditional GNNs are typically defined on
adjacency matrices, which cannot describe bidirected edges in causal diagrams.

  3. 3.

Like our work, this paper assumed that the causal graph $\mathcal{G}$ was
provided. However, in our work, we discuss the theoretical properties of the
implied constraints via Thms. 1 and 2. There does not appear to be a result in
the paper demonstrating that $\mathcal{G}$-GNNs or iVGAEs are
$\mathcal{G}^{(L_{2})}$-consistent. It is not immediately obvious that
$\mathcal{G}$-constrained GNNs have all of the $\mathcal{L}_{2}$ constraints
of $\mathcal{G}$ given the unique architecture of GNNs.

  4. 4.

The identification problem is one of the main focuses of our work. While the
paper discussed the identification problem of $\mathcal{L}_{2}$ queries given
$\mathcal{L}_{1}$ data, they did not provide a method of solving the problem
using iVGAEs. Rather, they focused on the estimation problem under the
assumption of identifiability. In fact, all $\mathcal{L}_{2}$ quantities are
identifiable from $\mathcal{L}_{1}$ data under the Markovianity assumption
(Bareinboim et al., 2022, Corol. 2), so an approach for solving the
identification problem in this setting may not be meaningful.

  5. 5.

Finally, we note that while this paper explored some novel aspects of using
GNN and VGAE architectures for causal inference, our paper focuses on GANs, an
entirely different class of models.

###  C.6 Sanchez-Martin et al. (2021)

Sanchez-Martin et al. (2021) introduced Variational Causal Autoencoders
(VACA), a model for estimating interventional and counterfactual quantities
under the Markovianity assumption using a variational graph autoencoder
architecture (VGAE) and graph neural network (GNN) architecture like (Zecevic
et al., 2021). VACA is defined as a tuple $\langle
A,p(\mathbf{Z}),p_{\theta},q_{\phi}\rangle$, where

  * •

$A$ is an adjacency matrix representing the causal DAG $\mathcal{G}$;

  * •

$p(\mathbf{Z})=\prod_{V_{i}\in\mathbf{V}}p(Z_{V_{i}})$ is a prior distribution
over independent latent variables
$\mathbf{Z}=\\{Z_{V_{i}}:V_{i}\in\mathbf{V}\\}$;

  * •

The decoder $p_{\theta}(\mathbf{V}\mid\mathbf{Z},A)$ is a GNN parameterized by
$\theta$ mapping $\mathbf{Z}$ and $A$ to the likelihood of $\mathbf{V}$ given
these values;

  * •

The encoder $q_{\phi}(\mathbf{Z}\mid\mathbf{V},A)$ is a GNN parameterized by
$\phi$ mapping $\mathbf{V}$ and $A$ to the posterior approximation of
$\mathbf{Z}$ given these values.

The model was trained to maximize the likelihood of observational data
$P(\mathbf{V})$, which can be expressed as

$P(\mathbf{V})\approx\int
p_{\theta}(\mathbf{V}\mid\mathbf{Z},A)p(\mathbf{Z})d\mathbf{Z}.$ (23)

Here, $\mathbf{Z}$ is similar to the role of $\mathbf{U}$ in an SCM. Although
$\mathbf{Z}$ does not have to be equal to $\mathbf{U}$ from the true SCM, each
$Z_{V_{i}}$ provides a source of noise for variable $V_{i}$ that captures the
same information that $\mathbf{U}_{V_{i}}$ contributes to $V_{i}$ in the true
SCM. The difference between $\mathbf{Z}$ and $\mathbf{U}$ does not matter as
long as the VACA model can match the true observational distribution.

During evaluation, interventional ($\mathcal{L}_{2}$) quantities can be
estimated via

$p(\mathbf{V}_{\mathbf{X}=\mathbf{x}})\approx\int\int
p_{\theta}(\mathbf{V}\mid\tilde{\mathbf{Z}},\tilde{\mathbf{Z}}^{\mathcal{I}},A^{\mathcal{I}})p(\tilde{\mathbf{Z}})q_{\phi}(\tilde{\mathbf{Z}}^{\mathcal{I}}\mid
A^{\mathcal{I}},\mathbf{X})d\tilde{\mathbf{Z}}d\tilde{\mathbf{Z}}^{\mathcal{I}},$
(24)

where $A^{\mathcal{I}}$ is the adjacency matrix with edges into variables of
$\mathbf{X}$ removed, $\tilde{\mathbf{Z}}^{\mathcal{I}}$ is the set of
representations of $\mathbf{X}$, and $\tilde{\mathbf{Z}}$ is the subset of
$\mathbf{Z}$ corresponding to $\mathbf{V}\setminus\mathbf{X}$. In words,
intervened variables are assigned a representation from the encoder under a
graph with the edges into $\mathbf{X}$ mutilated, while the remaining
variables are assigned the same representation.

Further, the following counterfactual ($\mathcal{L}_{3}$) quantity can also be
estimated:

$\displaystyle
P(\mathbf{V}_{\mathbf{X}=\mathbf{x}}=\mathbf{v}^{CF}\mid\mathbf{v}^{F})$
$\displaystyle\approx\int\int
p_{\theta}(\mathbf{V}=\mathbf{v}^{CF}\mid\tilde{\mathbf{Z}}^{F},\tilde{\mathbf{Z}}^{\mathcal{I}},A^{\mathcal{I}})q_{\phi}(\tilde{\mathbf{Z}}^{\mathcal{I}}\mid\mathbf{x},A^{\mathcal{I}})q_{\phi}(\tilde{\mathbf{Z}}^{F}\mid\mathbf{v}^{F},A)d\tilde{\mathbf{Z}}^{\mathcal{I}}d\tilde{\mathbf{Z}}^{F},$
(25)

where $\mathbf{v}^{F}$ is the observed factual value of $\mathbf{V}$,
$\mathbf{v}^{CF}$ is the counterfactual value of $\mathbf{V}$ to be queried,
and $\tilde{\mathbf{Z}}^{F}$ is the set of representations associated with
$\mathbf{v}^{F}$ except for the variables in $\mathbf{X}$. In words, the
3-step procedure for evaluating counterfactuals (Pearl, 2000, Thm. 7.1.7) is
applied here. The abduction step is performed by computing the factual
representations $\tilde{\mathbf{Z}}^{F}$ through the encoder
$q_{\phi}(\tilde{\mathbf{Z}}^{F}\mid\mathbf{v}^{F},A)$ applied to the factual
observations $\mathbf{v}^{F}$. The action step is performed by assigning the
intervened representations $\tilde{\mathbf{Z}}^{\mathcal{I}}$ as discussed in
Eq. 24. Finally, the prediction step is performed by using the interventional
representations $\tilde{\mathbf{Z}}^{\mathcal{I}}$ from $\mathbf{X}$ and the
factual representations $\tilde{\mathbf{Z}}^{F}$ from
$\mathbf{V}\setminus\mathbf{X}$ in the encoder
$p_{\theta}(\mathbf{V}\mid\tilde{\mathbf{Z}}^{F},\tilde{\mathbf{Z}},A^{\mathcal{I}})$
with the mutilated graph $A^{\mathcal{I}}$ to produce a counterfactual outcome
$\mathbf{v}^{CF}$ for $\mathbf{V}$ in this setting.

VACA presented interesting ideas for evaluating interventional and
counterfactual quantities in the Markovian case using GNNs and VGAEs. Still,
they differ from our work in several ways.

  1. 1.

Unlike this work, we do not make the Markovianity assumption, meaning we allow
unobserved confounding between variables. Many real world settings often have
unobserved confounding, so the Markovianity assumption can be limiting. In a
situation where the true SCM has unobserved confounding, a model class which
enforces Markovianity will not be expressive enough to model the true SCM
(Prop. 1). VACA does not extend trivially to cases with unobserved
confounders. Traditional GNNs are typically defined on adjacency matrices,
which cannot describe bidirected edges in causal diagrams.

  2. 2.

Like our work, this paper assumed that the causal graph $\mathcal{G}$ was
provided. However, in our work, we discuss the theoretical properties of the
implied constraints via Thms. 1 and 2. This paper possibly benefited from
these constraints implicitly without discussing these properties, but it is
not actually clear if VACA is a $\mathcal{G}^{(\mathcal{L}_{3})}$-consistent
model.

  3. 3.

The identification problem was not considered in this work while it is a major
focus in our work. As shown in Prop. 3 above, there exist counterfactual
quantities that cannot be identified given the graph $\mathcal{G}$, and all
available data from $\mathcal{L}_{1}$ and $\mathcal{L}_{2}$, even under the
Markovianity assumption. If these quantities were to be estimated from a model
trained on such data, the result will likely be incorrect and misleading. For
that reason, we make sure to identify all queries before attempting to
estimate them.

  4. 4.

While this paper only considered input in the form of observational data from
$\mathcal{L}_{1}$, our work considers cases where additional experimental
datasets from $\mathcal{L}_{2}$ are provided (even cases where perhaps the
observational distribution $P(\mathbf{V})$ is not provided). As shown in Prop.
2, there do exist counterfactual queries that are not identifiable given
observational data but are identifiable given datasets from $\mathcal{L}_{2}$.

  5. 5.

This paper evaluated counterfactual quantities using the 3-step procedure of
abduction, action, and prediction. While this is a sound algorithm for
evaluating applicable quantities, it does not apply generally to all
counterfactual quantities. Specifically, it only applies in quantities of the
form $P(\mathbf{Y}_{\mathbf{x}}\mid\mathbf{z})$, in which only one
intervention is applied before the conditioning bar, and observations are
present after the conditioning bar. This excludes quantities such as
$P(y_{x},y^{\prime}_{x^{\prime}})$, which is a joint probability of two
counterfactual terms in different worlds (one in $X=x$ and another in
$X=x^{\prime}$), as well as quantities such as $P(y_{x}\mid
y^{\prime}_{x^{\prime}})$, in which a term in an intervened world
$(X=x^{\prime})$ is after the conditioning bar.

  6. 6.

Finally, we note that while this paper explored some novel aspects of using
GNN and VGAE architectures for causal inference, our paper focuses on GANs, a
different class of models.

##  Appendix D Further Discussions and Examples

In this section, we supplement the main points of the paper with further
discussions and examples. In Section D.1, we expand on the impossibility of
performing inferences across layers of the PCH in general and why assumptions
are necessary. In particular, we discuss the implications of the causal
diagram assumption and compare it with alternatives. In Section D.2, we expand
on the reasoning behind introducing the discussion of constraints separately
for each layer (as in Def. 3). In Section D.3, we ground the theoretical
results of Secs. 2 and 3 by walking through examples with actual constructions
of NCMs. In Section D.4, we discuss the Markovianity assumption, why it
simplifies causal inference, whether it’s plausible, and why the approach
proposed in this paper is more general for _not_ making this assumption.

###  D.1 Necessity of Causal Assumptions

As emphasized at the beginning of Sec. 2, the entire discussion on constraints
in this paper is due to the Neural Causal Hierarchy Theorem (N-CHT) (Xia et
al., 2021, Corol. 1), following the Causal Hierarchy Theorem from Bareinboim
et al. (2022) that states that distributions from lower layers of the
hierarchy almost surely underdetermine higher layers in a measure-theoretic
sense. In other words, it is not possible to make a statement about a higher
layer (e.g. interventions, counterfactuals) with only knowledge from lower
layers (e.g. observations). This result is well-understood in the causal
inference field, and Pearl eloquently states in the Book of Why that “without
the causal model, we could not go from rung [layer] one to rung [layer] two.
This is why deep-learning systems (as long as they only use rung-one data and
do not have a causal model) will never be able to answer questions about
interventions (…)” (Pearl & Mackenzie, 2018, p. 32).

To further ground this discussion, consider the following examples.

######  Example 1.

Consider a scenario where the board of education of a city is deciding on
policies for school curricula. Data is collected from schools across the city,
noting the textbooks ($X$) used in the schools and the GPA of students from
that school ($Y$). For simplicity, let $X$ and $Y$ be binary variables, where
$X=1$ represents the use of a new textbook, and $Y=1$ represents a high
average GPA. Suppose the true model $\mathcal{M}^{*}$ describing the
relationship between $X$ and $Y$ is specified as follows:

$\mathcal{M}^{*}=\begin{cases}\mathbf{U}&=\\{U_{XY},U_{Y1},U_{Y2},U_{Y3}\\}\\\
\mathbf{V}&=\\{X,Y\\}\\\ \mathcal{F}&=\begin{cases}f_{X}(u_{XY})&=u_{XY}\\\
f_{Y}(x,u_{XY},u_{Y1},u_{Y2},u_{Y3})&=\begin{cases}\neg x\wedge
u_{Y1}&\text{if }u_{XY}=0\\\ (x\oplus u_{Y2})\vee u_{Y3}&\text{if
}u_{XY}=1\end{cases}\end{cases}\\\ P(\mathbf{U})&=U_{XY},U_{Y1}\sim
Bern(0.2),U_{Y2}\sim Bern(0.5),U_{Y3}\sim Bern(0.6)\end{cases}$ (26)

Here, $\neg$, $\wedge$, $\vee$, and $\oplus$ indicate bitwise negation, AND,
OR, and XOR operators, respectively. The board of education would like to see
if they should implement a policy to switch all schools to the new textbook.
Formally, they will decide to do this if

$P(Y_{X=1}=1)>P(Y_{X=0}=1)$ (27)

that is, if the effect of incorporating the new textbook on high GPA is higher
than not incorporating it.

$X$ | $Y$ | $P(X,Y)$
---|---|---
0 | 0 | 0.64
0 | 1 | 0.16
1 | 0 | 0.04
1 | 1 | 0.16
Table 2: Observational data for Example 1.

However, the board does not have access to the true SCM $\mathcal{M}^{*}$.
Instead, they have only collected observational data, shown in terms of
probabilities in Table 2. One can verify that these are the true observational
probabilities as induced by $\mathcal{M}^{*}$. Given this data, one may be
tempted to train a model to match this observational distribution. For
example, suppose the board’s data science team learns the following model:

$\widehat{M}=\begin{cases}\widehat{\mathbf{U}}&=\\{\widehat{U}_{X},\widehat{U}_{Y}\\}\\\
\mathbf{V}&=\\{X,Y\\}\\\
\widehat{\mathcal{F}}&=\begin{cases}\hat{f}_{X}(\widehat{u}_{X})&=\mathbbm{1}\\{\widehat{u}_{X}<0.2\\}\\\
\hat{f}_{Y}(x,\widehat{u}_{Y})&=\begin{cases}\mathbbm{1}\\{\widehat{u}_{Y}<0.2\\}&\text{if
}x=0\\\ \mathbbm{1}\\{\widehat{u}_{Y}<0.8\\}&\text{if
}x=1\end{cases}\end{cases}\\\
P(\widehat{\mathbf{U}})&=\widehat{U}_{X},\widehat{U}_{Y}\sim
Unif(0,1)\end{cases}$

One can verify that $\widehat{M}$ induces the same $P(\mathbf{V})$ and
therefore matches $\mathcal{M}^{*}$ in $\mathcal{L}_{1}$. Computing the
interventional effect from this model leads to:

$\displaystyle P^{\widehat{M}}(Y_{X=0}=1)=P(\widehat{U}_{Y}<0.2)=0.2,$ (28)
$\displaystyle P^{\widehat{M}}(Y_{X=1}=1)=P(\widehat{U}_{Y}<0.8)=0.8$ (29)

Since $P^{\widehat{M}}(Y_{X=1}=1)>P^{\widehat{M}}(Y_{X=0}=1)$, this seems to
suggest that implementing the new textbook is beneficial to the students’
population. However, computing the same values from $\mathcal{M}^{*}$ reveals
that

$\displaystyle P^{\mathcal{M}^{*}}(Y_{X=0}=1)=0.32,$ (30) $\displaystyle
P^{\widehat{M}}(Y_{X=1}=1)=0.16.$ (31)

In other words, implementing the new textbook would actually harm students’
performance on average.

The implemented model $\widehat{M}$ assigns $P(y_{x})=P(y\mid x)$, failing to
account for confounding between $X$ and $Y$. Still, even if confounding is
considered, it is unclear to what extent $X$ and $Y$ are confounded. In fact,
even the causal direction is not known. For example, perhaps schools with
better-performing students opted to use the new textbook to challenge the
students, resulting in a flipped causal relationship. Without further
information about the generating model, one cannot infer anything about the
interventional distribution $P(y_{x})$, even if the trained model matches the
true model on the observational level.,

$\blacksquare$

######  Example 2.

Consider again the board of education situation from Example 1. Suppose, after
conducting further experimental studies, the board discovers that (1) $Y$ has
no causal effect on $X$ ($P(x_{y})=P(x)$), and (2) $P(Y_{X=0}=1)=0.32$ and
$P(Y_{X=1}=1)=0.16$ as discovered in the previous example. Now, all
$\mathcal{L}_{2}$ information from $\mathcal{M}^{*}$ is provided as well. The
board realizes that implementing the new textbook may not be the best idea and
is now wondering whether schools that have already implemented the textbook
and have poorly performing students would see improvement if the new textbook
was no longer used. Formally, they are wondering about the counterfactual
probability $P(Y_{X=0}=1\mid X=1,Y=0)$. They construct the following model:

$\widehat{M}=\begin{cases}\widehat{\mathbf{U}}&=\\{\widehat{U}_{XY}\\}\\\
\mathbf{V}&=\\{X,Y\\}\\\
\widehat{\mathcal{F}}&=\begin{cases}\hat{f}_{X}(\widehat{u}_{XY})&=\mathbbm{1}\\{\widehat{u}_{XY}<0.2\\}\\\
\hat{f}_{Y}(x,\widehat{u}_{XY})&=\begin{cases}1&\text{if
}\widehat{u}_{XY}<0.16\\\ 0&\text{if }0.16\leq\widehat{u}_{XY}<0.84\\\ \neg
x&\text{otherwise}\end{cases}\end{cases}\\\
P(\widehat{\mathbf{U}})&=\widehat{U}_{XY}\sim Unif(0,1)\end{cases}$ (32)

Note that $\widehat{M}$ induces the same $\mathcal{L}_{1}$ and
$\mathcal{L}_{2}$ distributions as $\mathcal{M}^{*}$. When querying
$\widehat{M}$ with the query of interest, one obtains

$P^{\widehat{M}}(Y_{X=0}=1\mid X=1,Y=0)=0,$ (33)

indicating that removing the book from the curriculum would have no changes on
student performance. However, in reality,

$P^{\mathcal{M}^{*}}(Y_{X=0}=1\mid X=1,Y=0)=1,$ (34)

meaning that removing the book would guarantee an improvement.

Once again, despite fully matching $\mathcal{L}_{1}$ and $\mathcal{L}_{2}$ of
the true model, $\widehat{M}$ is not able to provide any meaningful
information about an $\mathcal{L}_{3}$ query. $\blacksquare$

These examples are written in SCM notation, but equivalent NCM constructions
can be found with the same distributions (following from (Xia et al., 2021,
Thm. 1)). These examples illustrate that there can be many models which match
on a lower layer but fail to match on a higher layer. Hence, a model (such as
an NCM) trained on data from a lower layer cannot be expected to produce
correct inferences from higher layers (i.e. crossing between layers is
impossible). For more examples, see Examples 7-9 from Bareinboim et al.
(2022). Also, Example 2 from Xia et al. (2021) shows a case where four
different NCMs all match in the given observational data but have drastically
differing outputs for an $\mathcal{L}_{2}$ query, and none of them are
correct. Evidently, NCMs cannot circumvent the impossibility described by the
CHT despite their expressiveness, as expressivity does not guarantee that the
learned content will match that of the true generating model.

Given this impossibility, the natural course of action to proceed in problems
of cross-layer inference is to study under what assumptions inferences could
potentially be made. Obviously, the problem could be solved if the mechanisms
and noise of the true generating SCM are known, but this is almost never the
case in practice. In this work, we assume that we are given partial knowledge
of the true generating SCM in the form of a causal diagram, a qualitative
description of the underlying mechanisms and noise. This assumption is
commonly used in the causal inference literature starting from Pearl (1995),
where the field started in its modern form, and later shown in several classic
causal inference applications such as in the identification of causal effects
with do-calculus (Pearl, 2000, Ch. 4) or counterfactual identification and
their axioms (Pearl, 2000, Ch. 7). There are compelling cases in which a data
analyst may have a qualitative understanding of the causal relationships of
the environment to form a causal graph (Bottou et al., 2013; Pearl, 2000).
These relationships are explored in several previous counterfactual studies
using symbolic approaches such as Galles & Pearl (1998), Shpitser & Pearl
(2007), and Correa et al. (2021), to cite a few.

The reason why the causal diagram assumption can be utilized to solve
counterfactual inference problems is due to the constraints implied by it.
While many SCMs may agree on lower layer quantities but disagree on higher
layer quantities, the pool of plausible models shrinks dramatically after
falsifying all SCMs that do not match the given constraints. Specifically
these constraints are the $\mathcal{L}_{i}$ equality constraints discussed in
Appendix A.1, which led to the notion of
$\mathcal{G}^{(\mathcal{L}_{i})}$-consistency discussed in Sec. 2. Still, even
with these assumptions, inferring counterfactual quantities is a non-trivial
task. Despite having a constrained setting, it is still not guaranteed that a
counterfactual quantity is computable from observational and interventional
data, nor is it obvious how to compute such a quantity even if it is
computable. This results in the counterfactual identification and estimation
problems discussed in Secs. 3 and 4.

$X$$\mathbf{Z}$$Y$ Figure 16: Graph with backdoor set $\mathbf{Z}$.

While assumptions are necessary to perform counterfactual inference, other
forms of assumptions aside from the causal diagram may also lead to positive
results. For example, many works in the intersection of neural models and
causal inference have achieved success under the assumption of conditional
ignorability (see references cited in Footnote 1). This assumption states
that, under treatment $X$, outcome $Y$, and co-variates $\mathbf{Z}$, we have

$X\mathbin{\mathchoice{\hbox
to0.0pt{\hbox{\set@color$\displaystyle\perp$}\hss}\kern 3.46875pt{}\kern
3.46875pt\hbox{\set@color$\displaystyle\perp$}}{\hbox
to0.0pt{\hbox{\set@color$\textstyle\perp$}\hss}\kern 3.46875pt{}\kern
3.46875pt\hbox{\set@color$\textstyle\perp$}}{\hbox
to0.0pt{\hbox{\set@color$\scriptstyle\perp$}\hss}\kern 2.36812pt{}\kern
2.36812pt\hbox{\set@color$\scriptstyle\perp$}}{\hbox
to0.0pt{\hbox{\set@color$\scriptscriptstyle\perp$}\hss}\kern 1.63437pt{}\kern
1.63437pt\hbox{\set@color$\scriptscriptstyle\perp$}}}Y_{x}\mid\mathbf{Z},$
(35)

or in words, $X$ is independent of the counterfactual variable $Y_{x}$ when
conditioned on $\mathbf{Z}$. This assumption is also known as back-door
admissibility and is encoded in causal diagrams like Fig. 16. In this paper,
we take the graphical language to validate assumptions, which may appear in
different forms than Eq. 35 above.

Another possibility is to assume the availability of all $\mathcal{L}_{2}$
data. This data, combined with the faithfulness assumption, can be used to
learn the true causal diagram (Kocaoglu et al., 2017). If this is the approach
used, then the work in this paper can be applied on top of this learned causal
diagram. There are examples of counterfactual queries that are still not
identifiable even when given the graph and all $\mathcal{L}_{2}$ information
(see Example 8). Still, this relies on the user having an abundance of
experimental data, which is not the most common setting in practice. The work
in this paper is agnostic to the approach for acquiring the causal diagram, so
it can still be used whether the graph is obtained from experiments or from
other means such as through expert knowledge.

Finally, although it is out of the scope of this work, one may envision ways
to relax the causal diagram assumption. One example is to work with an
equivalence class (EC) of causal diagrams which can be learned from purely
observational ($\mathcal{L}_{1}$) data (Jaber et al., 2020). Performing causal
inferences using the EC comes with its own merits and faces a different set of
challenges. In terms of merits, the work relies on fewer assumptions as it is
entirely data-driven. In terms of challenges, for example, many interventional
and counterfactual quantities are not identifiable since the equivalence class
is strictly weaker than the causal diagram itself (for obvious reasons). Such
works are complementary to this work since the data scientist can choose for
themselves whether they would prefer to make more assumptions for stronger
results.

Ultimately, the goal of this work is not to promote one type of assumption
over another but to offer a toolbox to solve counterfactual inference problems
once the assumptions are made. The results of this paper are simply the first
steps for solving counterfactual inferences in a well-understood setting. A
promising direction for future work is to further relax and generalize these
approaches when the diagram is not fully known.

###  D.2 Discussion on $\mathcal{L}_{i}$ Constraints

######  Example 3.

$X$$Z$$Y$ (a) True graph.

$X$$Z$$Y$ (b) Alternate graph.

Figure 17: Graphs for illustrating equality constraints.

To illustrate this concept of graphical constraints on different layers with
an example, consider the graphs in Fig. 17. Suppose graph (a) is the causal
diagram of the true SCM $\mathcal{M}^{*}$. Graph (a) implies different
constraints across all three layers. For example, on layer 1, we know by
d-separation that $P^{*}(y\mid z,x)=P^{*}(y\mid z)$. One could easily train a
model $\mathcal{M}$ given data on $P^{*}(X,Y,Z)$ such that the model outputs
$P^{\mathcal{M}}(X,Y,Z)=P^{*}(X,Y,Z)$, and hence, the equality constraint
$P^{\mathcal{M}}(y\mid z,x)=P^{\mathcal{M}}(y\mid z)$ holds in $\mathcal{M}$
as well. In this case $\mathcal{M}$ would be considered
$\mathcal{G}^{(\mathcal{L}_{1})}$-consistent with $\mathcal{M}^{*}$.

However, this is not to say that $\mathcal{M}$ takes the constraints of
$\mathcal{G}$ on higher layers. For example, $\mathcal{M}$ could be another
SCM with causal diagram (b), which also implies the constraint $P(y\mid
z,x)=P(y\mid z)$. The similarities break down on higher layer constraints. For
example, on layer 2, graph (a) implies that $P^{*}(x_{z})=P^{*}(x)$ (i.e. $Z$
has no causal effect on $X$) and $P^{*}(z_{x})=P^{*}(z\mid x)$ (i.e. the
causal effect of $X$ on $Z$ is equal to the association). If graph (b) truly
reflects the mechanisms of $\mathcal{M}$, then neither of these two
constraints would hold for $\mathcal{M}$ since the causal direction is
reversed. $\mathcal{M}$ would _not_ be
$\mathcal{G}^{(\mathcal{L}_{2})}$-consistent with $\mathcal{M}^{*}$, greatly
reducing its inferential capabilities on layer 2 quantities induced by
$\mathcal{M}^{*}$.

Layer 1 and layer 2 constraints are formalized through the concepts of
bayesian networks and causal bayesian networks, respectively. Although not yet
formally studied, there are further constraints encoded by $\mathcal{G}$
corresponding to layer 3 distributions. For example, graph (a) implies that
$P^{*}(y_{z}\mid z^{\prime})=P^{*}(y_{z})$ and that
$P^{*}(y_{x^{\prime}},z_{x})=P^{*}(y_{z},z_{x})$. These constraints are vital
for performing inference on layer 3, which is why
$\mathcal{G}^{(\mathcal{L}_{3})}$-consistency is an important property.

$\blacksquare$

###  D.3 Properties of the $\mathcal{G}$-NCM

Sec. 2 emphasizes that the $\mathcal{G}$-NCM has the properties of
$\mathcal{G}^{(\mathcal{L}_{3})}$-consistency and
$\mathcal{L}_{3}$-$\mathcal{G}$ expressiveness. In this section, we expand on
why these two properties are nontrivial.

The proof that the $\mathcal{G}$-NCM is
$\mathcal{G}^{(\mathcal{L}_{3})}$-consistent arises naturally from its
structure and the fact that it is a type of SCM. We first note that any other
type of structure may struggle to achieve
$\mathcal{G}^{(\mathcal{L}_{3})}$-consistency. To illustrate this point,
consider the following example.

######  Example 4.

$X$$Z$$Y$ Figure 18: Graph for following examples.

When representing a collection of distributions, perhaps the most simple model
that comes to mind is an expression wrapper.

######  Definition 12 (Expression Wrapper).

Given a collection of distributions $\mathbb{P}$, an _expression wrapper_
$M^{\mathbb{P}}$ defined over $\mathbb{P}$ is a collection of probability
models $M_{\mathbb{P}}=\\{M_{P}:P\in\mathbb{P}\\}$, where each $M_{P}$ is an
$\mathcal{L}_{1}$-expressive model such as a normalizing flow that can express
any distribution $P$. $\blacksquare$

It is obvious by construction that regardless of $\mathbb{P}$, there exists an
expression wrapper defined over $\mathbb{P}$ that can express all of the
distributions in $\mathbb{P}$. If we are interested in an application modeled
by true SCM $\mathcal{M}^{*}$, and
$\mathbb{P}=\mathcal{L}_{3}(\mathcal{M}^{*})$, that implies that there exists
an expression wrapper that is $\mathcal{L}_{3}$-consistent with
$\mathcal{M}^{*}$. The trouble arises when we aim to enforce
$G^{(\mathcal{L}_{3})}$-consistency. Consider the simple 3 variable case where
$\mathcal{G}$ is shown in Fig. 18. Even on layer 2, $\mathcal{G}$ encodes
several constraints, including but not limited to:

  1. 1.

$P(x)=P(x_{z})=P(x_{y})=P(x_{zy})$

  2. 2.

$P(z)=P(z_{y})$

  3. 3.

$P(y_{xz})=P(y_{z})=P(y\mid z)=P(y\mid x,z)$

  4. 4.

$P(z_{xy})=P(z_{x})=P(z\mid x)$

  5. 5.

$P(y_{x})=P(y\mid x)$

  6. 6.

$P(\\{x,z\\}_{y})=P(x,z)$

  7. 7.

$P(\\{y,z\\}_{x})=P(y,z\mid x)$

There are even more constraints on Layer 3, such as
$P(\\{x,z\\}_{y})=P(x,z)=P(x,z_{y})=P(x_{y},z)=P(x_{y},z_{y^{\prime}})=P(x_{y^{\prime}},z_{y})=P(x_{y^{\prime}},z_{y^{\prime}})=P(x_{y^{\prime}},z_{y^{\prime\prime}})=P(x_{z},z)=P(x_{z},z_{y})=P(x_{z},z_{y^{\prime}})=P(x_{z^{\prime}},z)=P(x_{z^{\prime}},z_{y})=P(x_{z^{\prime}},z_{y^{\prime}})=P(x_{z,y},z)=P(x_{z,y},z_{y})=P(x_{z,y},z_{y^{\prime}})=P(x_{z^{\prime},y},z)=P(x_{z^{\prime},y},z_{y})=P(x_{z^{\prime},y},z_{y^{\prime}})=P(x_{z,y^{\prime}},z)=P(x_{z,y^{\prime}},z_{y})=P(x_{z,y^{\prime}},z_{y^{\prime}})=P(x_{z,y^{\prime}},z_{y^{\prime\prime}})=P(x_{z^{\prime},y^{\prime}},z)=P(x_{z^{\prime},y^{\prime}},z_{y})=P(x_{z^{\prime},y^{\prime}},z_{y^{\prime}})=P(x_{z^{\prime},y^{\prime}},z_{y^{\prime\prime}})=P(x,z_{x})=P(x,z_{x,y})=P(x,z_{x,y^{\prime}})=P(x_{y},z_{x})=P(x_{y},z_{x,y})=P(x_{y},z_{x,y^{\prime}})=P(x_{y^{\prime}},z_{x})=P(x_{y^{\prime}},z_{x,y})=P(x_{y^{\prime}},z_{x,y^{\prime}})=P(x_{y^{\prime}},z_{x,y^{\prime\prime}})=P(x_{z},z_{x})=P(x_{z},z_{x,y})=P(x_{z},z_{x,y^{\prime}})=P(x_{z^{\prime}},z_{x})=P(x_{z^{\prime}},z_{x,y})=P(x_{z^{\prime}},z_{x,y^{\prime}})=P(x_{z,y},z_{x})=P(x_{z,y},z_{x,y})=P(x_{z,y},z_{x,y^{\prime}})=P(x_{z^{\prime},y},z_{x})=P(x_{z^{\prime},y},z_{x,y})=P(x_{z^{\prime},y},z_{x,y^{\prime}})=P(x_{z,y^{\prime}},z_{x})=P(x_{z,y^{\prime}},z_{x,y})=P(x_{z,y^{\prime}},z_{x,y^{\prime}})=P(x_{z,y^{\prime}},z_{x,y^{\prime\prime}})=P(x_{z^{\prime},y^{\prime}},z_{x})=P(x_{z^{\prime},y^{\prime}},z_{x,y})=P(x_{z^{\prime},y^{\prime}},z_{x,y^{\prime}})=P(x_{z^{\prime},y^{\prime}},z_{x,y^{\prime\prime}})$,
among several others. Furthermore, each of these constraints must hold for
every value of $x,y,z,x^{\prime},y^{\prime},z^{\prime}$, etc., in their
respective domains.

The sheer number of constraints even in a 3-variable case is so large that it
is difficult to list them all, let alone train a model to actually enforce all
of these constraints. The causal diagram parsimoniously encodes all of these
constraints implicitly, and therefore, the NCM enforces these constraints by
construction when it is used as an inductive bias. However, these constraints
must be enforced manually when using a model such as the expression wrapper.
Not only must every distribution $P\in\mathbb{P}$ be trained, which is already
a difficult task due to the size of $\mathbb{P}$, all of these exponentially
many constraints must be simultaneously enforced if
$\mathcal{G}^{(\mathcal{L}_{3})}$-consistency is desired. This is simply not
feasible. $\blacksquare$

There are also models that achieve
$\mathcal{G}^{(\mathcal{L}_{2})}$-consistency but not
$\mathcal{G}^{(\mathcal{L}_{3})}$-consistency. Consider the following example.

######  Example 5.

Recall the definition of Causal Bayesian networks.

######  Definition 13 (Causal Bayesian Network (CBN) (Bareinboim et al., 2022,
Def. 16)).

Given observed variables $\mathbf{V}$, let $\mathbf{P}^{*}$ be the collection
of all interventional distributions $P(\mathbf{V}\mid do(\mathbf{x}))$,
$\mathbf{X}\subseteq\mathbf{V}$, $\mathbf{x}\in\mathcal{D}_{\mathbf{X}}$. A
causal diagram $\mathcal{G}$ is a Causal Bayesian Network for $\mathbf{P}^{*}$
if for every intervention $do(\mathbf{X}=\mathbf{x})$ and every topological
ordering $<$ of $\mathcal{G}_{\overline{X}}$ through its directed edges,

  1. (i)

$P(\mathbf{V}\mid do(\mathbf{X}=\mathbf{x}))$ is semi-Markov relative to
$\mathcal{G}_{\overline{\mathbf{X}}}$.

  2. (ii)

For every $V_{i}\in\mathbf{V}\setminus\mathbf{X}$,
$\mathbf{W}\subseteq\mathbf{V}\setminus(\mathbf{Pa}_{i}^{\mathbf{X}+}\cup\mathbf{X}\cup\\{V_{i}\\})$:

$P(v_{i}\mid
do(\mathbf{x}),\mathbf{pa}_{i}^{\mathbf{x}+},do(\mathbf{w}))=P(v_{i}\mid
do(\mathbf{x}),\mathbf{pa}_{i}^{\mathbf{x}+})$

,

  3. (iii)

For every $V_{i}\in\mathbf{V}\setminus\mathbf{X}$, let
$\mathbf{Pa}_{i}^{\mathbf{X}+}$ be partitioned into two sets of confounded and
unconfounded parents, $\mathbf{Pa}_{i}^{c}$ and $\mathbf{Pa}_{i}^{u}$ in
$\mathcal{G}_{\overline{\mathbf{X}}}$. Then

$\displaystyle P(v_{i}\mid
do(\mathbf{x}),\mathbf{pa}_{i}^{c},do(\mathbf{pa}_{i}^{u}))$
$\displaystyle=P(v_{i}\mid
do(\mathbf{x}),\mathbf{pa}_{i}^{c},\mathbf{pa}_{i}^{u})$

Here,
$\mathbf{Pa}_{V_{i}}^{\mathbf{x}+}=\mathbf{Pa}(\\{V\in\mathbf{C}_{\overline{\mathbf{X}}}(V_{i}):V\leq
V_{i}\\})$, with $\mathbf{C}_{\overline{\mathbf{X}}}$ referring to the
corresponding C-component in $\mathcal{G}_{\overline{\mathbf{X}}}$ and $\leq$
referring to the topological ordering. $\blacksquare$

To clarify notation used in this definition, $P(\mathbf{y}\mid
do(\mathbf{x}),\mathbf{z})=P(\mathbf{y}_{\mathbf{x}}\mid\mathbf{z}_{\mathbf{x}})$
for any $\mathbf{X},\mathbf{Y},\mathbf{Z}\subseteq\mathbf{V}$, and
$\mathcal{G}_{\overline{\mathbf{X}}}$ refers to the graph $\mathcal{G}$ with
the edges into $\mathbf{X}$ removed.

In words, a causal diagram $\mathcal{G}$ is a CBN for a collection of
distributions $\mathbf{P}^{*}$ if $\mathbf{P}^{*}$ satisfies the
$\mathcal{L}_{2}$ equality constraints induced by $\mathcal{G}$. The CBN
explicitly states these constraints, and more details can be found in
Bareinboim et al. (2022).

Let $\mathcal{G}$ be some causal diagram induced by true SCM
$\mathcal{M}^{*}$. Let
$M_{\operatorname{\textsf{CBN}}}=\langle\mathcal{G},\mathbf{P}^{*}\rangle$ be
some CBN defined over distributions $\mathbf{P}^{*}$. Then, by definition,
$M_{\operatorname{\textsf{CBN}}}$ is
$\mathcal{G}^{(\mathcal{L}_{2})}$-consistent with $\mathcal{M}^{*}$. However,
$M_{\operatorname{\textsf{CBN}}}$ would not be
$\mathcal{G}^{(\mathcal{L}_{3})}$-consistent with $\mathcal{M}^{*}$ since the
constraints in the definition do not include any distributions outside of
layer 2. In fact, $\mathbf{P}^{*}$ does not even have to be defined with layer
3 distributions.

As an example, suppose $\mathcal{G}$ is the graph in Fig. 18, and that
$\mathcal{L}_{1}(\mathcal{M}^{*})=\mathcal{L}_{1}(M_{\operatorname{\textsf{CBN}}})$,
that is, $\mathbf{P}^{*}$ matches the layer 1 distribution of
$\mathcal{M}^{*}$. Then, since $M_{\operatorname{\textsf{CBN}}}$ is
$\mathcal{G}^{(\mathcal{L}_{2})}$-consistent, one can derive $P(y_{x})=P(y\mid
x)$ and conclude that
$P^{M_{\operatorname{\textsf{CBN}}}}(y_{x})=P^{\mathcal{M}^{*}}(y_{x})$,
successfully identifying a layer 2 quantity. However, if a model were
$\mathcal{G}^{(\mathcal{L}_{3})}$-consistent with $\mathcal{M}^{*}$, one could
derive quantities such as
$P(y_{x},z_{x})=P(y_{xz},z_{x})=P(y_{z},z_{x})=P(y_{z})P(z_{x})=P(y\mid
z)P(z\mid x)$, also computable from observational quantities. This is not
guaranteed to be true with the CBN, since counterfactual quantities like
$P(y_{z},z_{x})$ are not guaranteed to be defined in $\mathbf{P}^{*}$, and
even if they were, they could take any values since no part of the definition
of CBNs constrains them. $\blacksquare$

There also exist models that may be
$\mathcal{G}^{(\mathcal{L}_{3})}$-consistent but are not
$\mathcal{L}_{3}$-$\mathcal{G}$ expressive. Any model class in which
Markovianity is enforced is an example of such a case, as shown in the
following.

######  Example 6.

$X$$Y$ (a) True graph.

$X$$Y$ (b) Markovian subgraph.

Figure 19: Graphs for Example 6.

Given the causal diagram $\mathcal{G}$ from the true SCM $\mathcal{M}^{*}$,
consider $\widetilde{\Omega}(\mathcal{G})$, the set of all Markovian SCMs that
induce causal diagrams with a subset of the edges of $\mathcal{G}$. By
construction, $\widetilde{\Omega}(\mathcal{G})$ is
$\mathcal{G}^{(\mathcal{L}_{3})}$-consistent with $\mathcal{M}^{*}$ since no
constraint can be removed by removing edges, so any SCM in
$\widetilde{\Omega}(\mathcal{G})$ has at least as many $\mathcal{L}_{3}$
constraints as those induced by $\mathcal{G}$.

Suppose $\mathcal{G}$ is the bow graph as shown in Fig. 19, and the true SCM
$\mathcal{M}^{*}$ is defined as

$\mathcal{M}^{*}=\begin{cases}\mathbf{U}&=\\{U_{XY},U_{Y}\\},\text{ both
binary}\\\ \mathbf{V}&=\\{X,Y\\},\text{ both binary}\\\
\mathcal{F}&=\begin{cases}f_{X}(u_{XY})&=u_{XY}\\\
f_{Y}(x,u_{XY},u_{Y})&=u_{XY}\oplus u_{Y}\end{cases}\\\ P(\mathbf{U})&\text{
defined such that
}U_{XY}\sim\operatorname{Bernoulli}(0.5),U_{Y}\sim\operatorname{Bernoulli}(0.1)\end{cases}$

Evidently, $X$ and $Y$ are highly correlated, as $P(Y=1\mid X=1)=0.9$, but
$P(Y=1)=0.5$. However, this correlation is not causal, as the input $X$
remains unused by $f_{Y}$, implying that $P(Y_{X=1}=1)=P(Y=1)=0.5$. It is
clear that no model in $\widetilde{\Omega}(\mathcal{G})$ will be able to
represent $\mathcal{M}^{*}$ on all three layers. For example, one may
construct $\mathcal{M}^{\prime}\in\widetilde{\Omega}(\mathcal{G})$ as follows

$\mathcal{M}^{\prime}=\begin{cases}\mathbf{U}^{\prime}&=\\{U_{X},U_{Y}\\},\text{
both binary}\\\ \mathbf{V}&=\\{X,Y\\},\text{ both binary}\\\
\mathcal{F}^{\prime}&=\begin{cases}f^{\prime}_{X}(u_{X})&=u_{X}\\\
f^{\prime}_{Y}(x,u_{Y})&=x\oplus u_{Y}\end{cases}\\\
P(\mathbf{U}^{\prime})&\text{ defined such that
}U_{X}\sim\operatorname{Bernoulli}(0.5),U_{Y}\sim\operatorname{Bernoulli}(0.1)\end{cases}$

This model induces the graph $\mathcal{G}^{\prime}$ as shown in Fig. 19(b), so
it is a valid model from $\widetilde{\Omega}(\mathcal{G})$. Moreover, it can
even reproduce $\mathcal{L}_{1}$ quantities from $\mathcal{M}^{*}$, as it is
clear that $P^{\mathcal{M}^{\prime}}(Y=1\mid X=1)=0.9$, and
$P^{\mathcal{M}^{\prime}}(Y=1)=0.5$. This makes $\mathcal{M}^{\prime}$
$\mathcal{L}_{1}$-consistent with $\mathcal{M}^{*}$. However, on layer 2, the
similarities break down. While in the true model, $X$ has no effect on $Y$, in
this case, the Markovian nature of $\mathcal{M}^{\prime}$ means that
$P^{\mathcal{M}^{\prime}}(Y_{X=1}=1)=P^{\mathcal{M}^{\prime}}(Y=1\mid
X=1)=0.9$, resulting in a different value.

Clearly $\mathcal{M}^{\prime}$ is not $\mathcal{L}_{2}$-consistent with
$\mathcal{M}^{*}$, let alone $\mathcal{L}_{3}$-consistent. However, no other
model from $\widetilde{\Omega}(\mathcal{G})$ would be able to achieve
$\mathcal{L}_{2}$-consistency either. In a Markovian model, the correlation
between $X$ and $Y$ must be modeled through a causal effect from either $X$ to
$Y$ or from $Y$ to $X$, and both of these options would violate
$\mathcal{L}_{2}$-consistency with $\mathcal{M}^{*}$. $\blacksquare$

Other examples of model classes that are
$\mathcal{G}^{(\mathcal{L}_{3})}$-consistent but not
$\mathcal{L}_{3}$-$\mathcal{G}$ expressive include parametric SCMs in which
the parameterization is not fully expressive, such as linear SCMs.

As shown in Thm. 2, the class of $\mathcal{G}$-NCMs is
$\mathcal{L}_{3}$-$\mathcal{G}$ expressive. It may not be immediately clear
why this expressivity holds on all three layers instead of only on layers 1 or
2. The intuition is that with the expressivity of neural networks, a
$\mathcal{G}$-NCM can directly model the functions of the true SCM
$\mathcal{M}^{*}$ instead of indirectly matching its distributions. By
matching the functions and utilizing the probability integral transform on
uniform random variables to match $P(\mathbf{U})$, it is guaranteed that one
can construct a $\mathcal{G}$-NCM to directly mimic the behavior of
$\mathcal{M}^{*}$ on the functional level, implying equivalence on all three
layers. Consider the following example.

######  Example 7.

$X$$Y$ Figure 20: Graph for Example 7.

Suppose we are studying the effects of a drug ($X$) on the recovery of a
disease ($Y$). For simplicity, $X$ and $Y$ are both binary, $X=1$ represents
the drug being taken, and $Y=1$ represents a full recovery. There are two
hypothesized models. In $\mathcal{M}_{1}$, taking the drug helps people
without a certain gene ($U_{Y}=1$) but actively harms people with the gene
($U_{Y}=0$). In $\mathcal{M}_{2}$, the drug has no effect. The two SCMs are
clarified as follows:

$\mathcal{M}_{1}=\begin{cases}\mathbf{U}_{1}&=\\{U_{X},U_{Y}\\},\text{ both
binary}\\\ \mathbf{V}&=\\{X,Y\\},\text{ both binary}\\\
\mathcal{F}_{1}&=\begin{cases}f^{1}_{X}(u_{X})&=u_{X}\\\
f^{1}_{Y}(x,u_{Y})&=x\oplus u_{Y}\end{cases}\\\ P(\mathbf{U}_{1})&\text{
defined such that }U_{X},U_{Y}\sim\operatorname{Bernoulli}(0.5)\end{cases}$
$\mathcal{M}_{2}=\begin{cases}\mathbf{U}_{2}&=\\{U_{X},U_{Y}\\},\text{ both
binary}\\\ \mathbf{V}&=\\{X,Y\\},\text{ both binary}\\\
\mathcal{F}_{2}&=\begin{cases}f^{2}_{X}(u_{X})&=u_{X}\\\
f^{2}_{Y}(x,u_{Y})&=u_{Y}\end{cases}\\\ P(\mathbf{U}_{2})&\text{ defined such
that }U_{X},U_{Y}\sim\operatorname{Bernoulli}(0.5)\end{cases}$

It turns out that these two hypothesized models match in both
$\mathcal{L}_{1}$ and $\mathcal{L}_{2}$. For example,
$P^{\mathcal{M}_{1}}(Y_{X=1}=1)=P^{\mathcal{M}_{2}}(Y_{X=1}=1)=0.5$. The
difference lies in their induced distributions for $\mathcal{L}_{3}$. For
example, $P^{\mathcal{M}_{1}}(Y_{X=1}=1\mid X=0,Y=1)=0$. In words, in
$\mathcal{M}_{1}$, if it is observed that someone did not take the drug and
recovered from the disease, then it is guaranteed that they would fail to
recover if they had taken the drug. On the other hand,
$P^{\mathcal{M}_{2}}(Y_{X=1}=1\mid X=0,Y=1)=P^{\mathcal{M}_{2}}(Y=1\mid
Y=1)=1$ since $X$ simply has no effect on the disease.

The reasoning for this difference at the functional level is that both
$f^{1}_{Y}$ and $f^{2}_{Y}$ have two different behaviors with respect to $X$
depending on the setting of $U_{Y}$. When $U_{Y}=0$, $f^{1}_{Y}(x)=x$, while
$f^{2}_{Y}(x)=0$. When $U_{Y}=1$, $f^{1}_{Y}(x)=1-x$, while $f^{2}_{Y}(x)=1$.
In both cases, the functions are different between $\mathcal{M}_{1}$ and
$\mathcal{M}_{2}$. Notably, the output of $f^{1}_{Y}$ depends on $X$ while the
output of $f^{2}_{Y}$ does not.

We can reproduce this difference in behavior when constructing
$\mathcal{G}$-NCMs to represent these SCMs.

Recall the binary step activation function,

$\sigma(z)=\begin{cases}0&z<0\\\ 1&z\geq 0.\end{cases}$ (36)

Using this activation function, we can create various useful neural network
functions like the following:

$\displaystyle\hat{f}_{OR}(z_{1},\dots,z_{n})$
$\displaystyle=\sigma\left(\sum_{i=1}^{n}z_{i}-1\right)$ bitwise OR
$\displaystyle\hat{f}_{AND}(z_{1},\dots,z_{n})$
$\displaystyle=\sigma\left(\sum_{i=1}^{n}z_{i}-n\right)$ bitwise AND
$\displaystyle\hat{f}_{NOT}(z)$ $\displaystyle=\sigma(-z)$ bitwise NOT
$\displaystyle\hat{f}_{\geq p}(z)$ $\displaystyle=\sigma(z-p)$
$\displaystyle\text{1 if }z\geq p$ $\displaystyle\hat{f}_{<p}(z)$
$\displaystyle=\hat{f}_{NOT}(\hat{f}_{\geq p}(z))$ $\displaystyle\text{1 if
}z<p$

With these functions, consider the following $\mathcal{G}$-NCMs,
$\widehat{M}_{1}$ and $\widehat{M}_{2}$ that replicate the behavior of
$\mathcal{M}_{1}$ and $\mathcal{M}_{2}$ respectively.

$\widehat{M}_{1}=\begin{cases}\widehat{\mathbf{U}}_{1}&=\\{\widehat{U}_{X},\widehat{U}_{Y}\\}\\\
\mathbf{V}&=\\{X,Y\\},\text{ both binary}\\\
\widehat{\mathcal{F}}_{1}&=\begin{cases}\hat{f}^{1}_{X}(\widehat{u}_{X})&=\hat{f}_{<0.5}(\widehat{u}_{X})\\\
\hat{f}^{1}_{Y}(x,\widehat{u}_{Y})&=\hat{f}_{OR}\begin{cases}\hat{f}_{AND}(x,\hat{f}_{NOT}(\hat{f}_{<0.5}(\widehat{u}_{Y})))\\\
\hat{f}_{AND}(\hat{f}_{NOT}(x),\hat{f}_{<0.5}(\widehat{u}_{Y}))\end{cases}\end{cases}\\\
P(\widehat{\mathbf{U}}_{1})&\text{ defined such that
}\widehat{U}_{X},\widehat{U}_{Y}\sim\operatorname{Unif}(0,1)\end{cases}$
$\widehat{M}_{2}=\begin{cases}\widehat{\mathbf{U}}_{2}&=\\{\widehat{U}_{X},\widehat{U}_{Y}\\}\\\
\mathbf{V}&=\\{X,Y\\},\text{ both binary}\\\
\mathcal{F}_{2}&=\begin{cases}\hat{f}^{2}_{X}(\widehat{u}_{X})&=\hat{f}_{<0.5}(\widehat{u}_{X})\\\
\hat{f}^{2}_{Y}(x,\widehat{u}_{Y})&=\hat{f}_{<0.5}(\widehat{u}_{Y})\end{cases}\\\
P(\widehat{\mathbf{U}}_{2})&\text{ defined such that
}\widehat{U}_{X},\widehat{U}_{Y}\sim\operatorname{Unif}(0,1)\end{cases}$

Note that $\widehat{M}_{1}$ and $\widehat{M}_{2}$ match $\mathcal{M}_{1}$ and
$\mathcal{M}_{2}$, respectively, on all three layers. For example,
$P^{\widehat{M}_{1}}(Y_{X=1}=1\mid X=0,Y=1)=0$ and
$P^{\widehat{M}_{2}}(Y_{X=1}=1\mid X=0,Y=1)=1$. This is because the
construction of $\widehat{M}_{1}$ and $\widehat{M}_{2}$ were designed to
directly match the functions of $\mathcal{M}_{1}$ and $\mathcal{M}_{2}$.
$\blacksquare$

As a disclaimer, note that $\mathcal{L}_{3}$-$\mathcal{G}$ expressiveness
means the existence of a model that matches on all three layers, but it does
not necessarily mean that such a model can be found. In the case of Example 7,
we could create $\mathcal{G}$-NCMs to match the hypothesized SCMs because the
full specifications of the SCMs were provided. In practice, this is rarely the
case, and it may not be possible to determine the correct NCM from training on
data alone, which is why the identification problem is important to determine
when inferences can be made from a trained NCM.

Determining whether a counterfactual quantity is identifiable from the given
data and constraints is an important problem even in the Markovian case.
Consider the following simple example.

######  Example 8.

$X$$Y$ Figure 21: Graph for Example 8.

Consider a simple study on the effect of a new drug ($X$) on survival rate of
a sick patient ($Y$), modeled by the following SCM $\mathcal{M}$:

$\mathcal{M}^{*}=\begin{cases}\mathbf{U}&=\\{U_{X},U_{Y1},U_{Y2}\\},\text{ all
binary}\\\ \mathbf{V}&=\\{X,Y\\},\text{ both binary}\\\
\mathcal{F}&=\begin{cases}f_{X}(u_{X})&=u_{X}\\\
f_{Y}(x,u_{Y1},u_{Y2})&=\begin{cases}x\oplus u_{Y1}&u_{Y2}=0\\\
u_{Y1}&u_{Y2}=1\end{cases}\end{cases}\\\ P(\mathbf{U})&\text{ defined such
that }U_{X},U_{Y1},U_{Y2}\sim\operatorname{Bernoulli}(0.5)\end{cases}$

This SCM induces the graph $\mathcal{G}$ in Fig. 21. Evidently from observing
$\mathcal{M}^{*}$, this drug has some effect on patients with this disease.
However, $\mathcal{M}^{*}$ is unobserved, and we only have data. Suppose we
have both observational data and experimental data from testing on $X$.
Experiments on $Y$ are not needed since $Y$ has no causal effect on $X$, so
this essentially means we have data on all of layer 1 and layer 2. This data
is shown in terms of probabilities in Table. 3.

$X$ | $Y$ | $P(X,Y)$
---|---|---
0 | 0 | 0.25
0 | 1 | 0.25
1 | 0 | 0.25
1 | 1 | 0.25

$x$ | $Y_{X=x}$ | $P(Y_{X=x})$
---|---|---
0 | 0 | 0.5
0 | 1 | 0.5
1 | 0 | 0.5
1 | 1 | 0.5

Table 3: Data for Example 8

The data seems to suggest that the drug has no effect on the disease since
$P(y_{x})=P(y\mid x)=P(y)$ for all values of $y$ and $x$. However, just to be
sure, we are interested in evaluating the query $P(Y_{X=1}=1\mid X=0,Y=0)$,
the probability that someone who did not take the drug and did not recover
_would have_ recovered had they taken the drug. This quantity is called the
probability of sufficiency (Pearl, 2000, Def. 9.2.2).

Running the procedure in Alg. 1 on $\mathcal{G}$ and
$\mathbb{Z}=\\{P(\mathbf{V}),P(\mathbf{V}_{x})\\}$, we may achieve the
following two $\mathcal{G}$-NCMs (with the notation from Example 7)

$\widehat{M}_{1}=\begin{cases}\widehat{\mathbf{U}}&=\\{\widehat{U}_{X},\widehat{U}_{Y}\\}\\\
\mathbf{V}&=\\{X,Y\\}\\\
\widehat{\mathcal{F}}_{1}&=\begin{cases}\hat{f}_{X}^{(1)}(\widehat{u}_{X})&=\hat{f}_{<0.5}(\widehat{u}_{X})\\\
\hat{f}_{Y}^{(1)}(x,\widehat{u}_{Y})&=\hat{f}_{<0.5}(\widehat{u}_{Y})\end{cases}\\\
P(\widehat{\mathbf{U}})&\text{ defined such that
}\widehat{U}_{X},\widehat{U}_{Y}\sim\operatorname{Unif}(0,1)\end{cases}$

$\widehat{M}_{2}=\begin{cases}\widehat{\mathbf{U}}&=\\{\widehat{U}_{X},\widehat{U}_{Y}\\}\\\
\mathbf{V}&=\\{X,Y\\}\\\
\widehat{\mathcal{F}}_{2}&=\begin{cases}\hat{f}_{X}^{(2)}(\widehat{u}_{X})&=\hat{f}_{<0.5}(\widehat{u}_{X})\\\
\hat{f}_{Y}^{(2)}(x,\widehat{u}_{Y})&=\hat{f}_{OR}(\hat{f}_{AND}(x,\hat{f}_{<0.5}(\widehat{u}_{Y})),\hat{f}_{AND}(\hat{f}_{NOT}(x),\hat{f}_{\geq
0.5}(\widehat{u}_{Y})))\end{cases}\\\ P(\widehat{\mathbf{U}})&\text{ defined
such that
}\widehat{U}_{X},\widehat{U}_{Y}\sim\operatorname{Unif}(0,1)\end{cases}$

Note that $\hat{f}_{Y}^{(2)}$ simply resembles the XOR function. Evaluating
$P^{\widehat{M}_{1}}(X,Y)$, $P^{\widehat{M}_{2}}(X,Y)$,
$P^{\widehat{M}_{1}}(Y_{x})$, and $P^{\widehat{M}_{2}}(Y_{x})$, one can verify
that they match the results in Table 3. However, for our query, observe that

$\displaystyle P^{\widehat{M}_{1}}(Y_{X=1}=1\mid X=0,Y=0)$
$\displaystyle=\frac{P^{\widehat{M}_{1}}(Y_{X=1}=1,X=0,Y=0)}{P^{\widehat{M}_{1}}(X=0,Y=0)}$
$\displaystyle=\frac{P(\widehat{U}_{Y}<0.5,\widehat{U}_{X}\geq
0.5,\widehat{U}_{Y}\geq 0.5)}{P(\widehat{U}_{X}\geq 0.5,\widehat{U}_{Y}\geq
0.5)}$ $\displaystyle=\frac{0}{0.25}=0,$

and

$\displaystyle P^{\widehat{M}_{2}}(Y_{X=1}=1\mid X=0,Y=0)$
$\displaystyle=\frac{P^{\widehat{M}_{2}}(Y_{X=1}=1,X=0,Y=0)}{P^{\widehat{M}_{2}}(X=0,Y=0)}$
$\displaystyle=\frac{P(\widehat{U}_{Y}<0.5,\widehat{U}_{X}\geq
0.5,\widehat{U}_{Y}<0.5)}{P(\widehat{U}_{X}\geq 0.5,\widehat{U}_{Y}<0.5)}$
$\displaystyle=\frac{0.25}{0.25}=1.$

In words, according to $\widehat{M}_{1}$, someone who did not take the drug or
recover is guaranteed _not_ to recover even if they had taken the drug. On the
other hand, according to $\widehat{M}_{2}$, the opposite is true, and that
person would be guaranteed to recover if they had taken the drug. Remarkably,
despite the fact that we have the graph $\mathcal{G}$ and _all_ data from
$\mathcal{L}_{1}$ and $\mathcal{L}_{2}$, we cannot narrow down the value of
our query at all, as it can take any possible value between 0 and 1.

Moreover, when we evaluate the true result of the query from
$\mathcal{M}^{*}$, we get

$\displaystyle P^{\mathcal{M}^{*}}(Y_{X=1}=1\mid X=0,Y=0)$
$\displaystyle=\frac{P^{\mathcal{M}^{*}}(Y_{X=1}=1,X=0,Y=0)}{P^{\mathcal{M}^{*}}(X=0,Y=0)}$
$\displaystyle=\frac{P(U_{X}=0,U_{Y1}=0,U_{Y2}=0)}{P(U_{X}=0,U_{Y1}=0)}$
$\displaystyle=\frac{0.125}{0.25}=\frac{1}{2}.$

In other words, neither $\widehat{M}_{1}$ nor $\widehat{M}_{2}$ captured the
true value of the query. This query is not identifiable from $\mathbb{Z}$ and
$\mathcal{G}$. Arbitrarily training a model to fit $\mathbb{Z}$ and
$\mathcal{G}$ is almost guaranteed to fail to capture the true query in this
case.

$\blacksquare$

Although Example 8 shows a compelling case where one would not want to
estimate the query due to non-identifiability, the NCM serves as a robust
proxy model in cases where the query is identifiable.

######  Example 9.

$R$$J$$E$ Figure 22: Graph for Example 9.

Suppose we are conducting a study on whether race ($R$) affects an applicant’s
chances of getting a certain job ($J$). We also collect data on the
applicant’s education background ($E$). For simplicity, suppose all three are
binary variables, $R=1$ implies that the applicant is from a protected class
(e.g. black, hispanic), $E=1$ implies that the applicant graduated from a
prestigious university, and $J=1$ implies that the applicant got the job.
Experts suggest that race may affect application results both directly and
indirectly through education background, but race and education may share
confounding factors such as location, as shown in Fig. 22. The full
relationship between variables is modeled through SCM $\mathcal{M}^{*}$:

$\mathcal{M}^{*}=\begin{cases}\mathbf{U}&=\\{U_{RE},U_{E},U_{J}\\},\text{ all
binary}\\\ \mathbf{V}&=\\{R,E,J\\}\\\
\mathcal{F}&=\begin{cases}f_{R}(u_{RE})&=u_{RE}\\\
f_{E}(r,u_{RE},u_{E})&=(\neg u_{RE}\wedge\neg r)\oplus u_{E}\\\
f_{J}(r,e,u_{J})&=(\neg r\vee e)\oplus u_{J}\end{cases}\\\
P(\mathbf{U})&\text{ defined such that
}U_{RE}\sim\operatorname{Bernoulli}(0.5),U_{E},U_{J}\sim\operatorname{Bernoulli}(0.25)\end{cases}$

Relevant quantities and their corresponding probability for each setting of
$\mathbf{U}$ is shown in Table 4.

$P(\mathbf{U})$ | $U_{RE}$ | $U_{E}$ | $U_{J}$ | $R$ | $E$ | $J$ | $E_{R=0}$ | $E_{R=1}$ | $J_{R=0}$ | $J_{R=1}$ | $J_{R=1,E_{R=0}}$
---|---|---|---|---|---|---|---|---|---|---|---
$p_{0}=9/32$ | 0 | 0 | 0 | 0 | 1 | 1 | 1 | 0 | 1 | 0 | 1
$p_{1}=3/32$ | 0 | 0 | 1 | 0 | 1 | 0 | 1 | 0 | 0 | 1 | 0
$p_{2}=3/32$ | 0 | 1 | 0 | 0 | 0 | 1 | 0 | 1 | 1 | 1 | 0
$p_{3}=1/32$ | 0 | 1 | 1 | 0 | 0 | 0 | 0 | 1 | 0 | 0 | 1
$p_{4}=9/32$ | 1 | 0 | 0 | 1 | 0 | 0 | 0 | 0 | 1 | 0 | 0
$p_{5}=3/32$ | 1 | 0 | 1 | 1 | 0 | 1 | 0 | 0 | 0 | 1 | 1
$p_{6}=3/32$ | 1 | 1 | 0 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1
$p_{7}=1/32$ | 1 | 1 | 1 | 1 | 1 | 0 | 1 | 1 | 0 | 0 | 0
Table 4: Relevant quantities calculated for each setting of $\mathbf{U}$ in
$\mathcal{M}^{*}$ for Example 9.

We are given the causal diagram $\mathcal{G}$ from Fig. 22 as well as the
datasets
$\mathbb{Z}=\\{P(\mathbf{V}),P(\mathbf{V}_{R=0}),P(\mathbf{V}_{R=1})\\}$. The
observational data from $P(\mathbf{V})$ was collected by simply asking random
participants to list their race, education level, and application results. The
experimental data from $P(\mathbf{V}_{R=0})$ and $P(\mathbf{V}_{R=1})$ is
obtained similarly, except participants are asked to write their race as $R=0$
or $R=1$ respectively in their school and job applications regardless of their
actual race. This data, computed from Table 4, is shown in Table 5.

$R$ | $E$ | $J$ | $P(\mathbf{V})$
---|---|---|---
0 | 0 | 0 | $p_{3}=1/32$
0 | 0 | 1 | $p_{2}=3/32$
0 | 1 | 0 | $p_{1}=3/32$
0 | 1 | 1 | $p_{0}=9/32$
1 | 0 | 0 | $p_{4}=9/32$
1 | 0 | 1 | $p_{5}=3/32$
1 | 1 | 0 | $p_{7}=1/32$
1 | 1 | 1 | $p_{6}=3/32$

$r$ | $E_{R=r}$ | $J_{R=r}$ | $P(\mathbf{V}_{R=r})$
---|---|---|---
0 | 0 | 0 | $p_{3}+p_{5}=4/32$
0 | 0 | 1 | $p_{2}+p_{4}=12/32$
0 | 1 | 0 | $p_{1}+p_{7}=4/32$
0 | 1 | 1 | $p_{0}+p_{6}=12/32$
1 | 0 | 0 | $p_{0}+p_{4}=18/32$
1 | 0 | 1 | $p_{1}+p_{5}=6/32$
1 | 1 | 0 | $p_{3}+p_{7}=2/32$
1 | 1 | 1 | $p_{2}+p_{6}=6/32$

Table 5: Data for Example 9

When looking at the average treatment effect of race on job success from the
data, we see that

$P(J_{R=1}=1)-P(J_{R=0}=1)=\frac{12}{32}-\frac{24}{32}=-\frac{12}{32}=-0.375$

which is a rather large disparity between the success rates of protected
classes compared to others. To avoid accusations of discrimination, the
company claims that this disparity is not a result of the company’s hiring
practices. The company says that they only look at education levels, and the
reason why protected classes are less likely to get the job is because the
fault lies with the education system for discriminating against protected
classes in their admissions process. To verify this claim, we may be
interested in the _natural direct effect_ query, formally defined in this
context as

$Q=P(J_{R=1,E_{R=0}}=1)-P(J_{R=0}=1).$ (37)

In words, the first term is the probability that someone will get the job if
we change their race to a protected class, but change their education level as
if their race was not from a protected class. The second term is simply their
natural chance of getting the job if their race was fixed to a non-protected
class. In other words, this query calculates the direct effect of race on the
hiring process, irrespective of education. We do not have data on the first
term of $Q$, but it turns out that $Q$ is identifiable from $\mathcal{G}$ and
$\mathbb{Z}$. Running Alg. 1 would return two NCMs that induce the same result
for $Q$. For example, the first term can be derived symbolically as

$\displaystyle P(J_{R=1,E_{R=0}}=1)$
$\displaystyle=P(J_{R=1,E=0}=1,E_{R=0}=0)+P(J_{R=1,E=1}=1,E_{R=0}=1)$
Unnesting
$\displaystyle=P(J_{R=1,E=0}=1)P(E_{R=0}=0)+P(J_{R=1,E=1}=1)P(E_{R=0}=1)$
C-factor decomp. $\displaystyle=P(J=1\mid R=1,E=0)P(E_{R=0}=0)$ Rule 2
$\displaystyle+P(J=1\mid R=1,E=1)P(E_{R=0}=1)$
$\displaystyle=\left(\frac{3}{12}\right)\left(\frac{16}{32}\right)+\left(\frac{3}{4}\right)\left(\frac{16}{32}\right)=\frac{1}{8}+\frac{3}{8}=\frac{1}{2}$
From data

From the data, we also have $P(J_{R=0}=1)=\frac{24}{32}=\frac{3}{4}$. Hence,
we can compute

$\displaystyle
Q=P(J_{R=1,E_{R=0}}=1)-P(J_{R=0}=1)=\frac{1}{2}-\frac{3}{4}=-\frac{1}{4}=-0.25,$

revealing that while the direct effect is slightly lower than the average
treatment effect, the company still shares some responsibility in the
discrimination of protected classes in their hiring practices. Indeed, one
could construct a $\mathcal{G}$-NCM that matches $\mathcal{M}^{*}$ on the
given data sets, for example:

$\widehat{M}=\begin{cases}\widehat{\mathbf{U}}&=\\{\widehat{U}_{RE},\widehat{U}_{J}\\}\\\
\mathbf{V}&=\\{R,E,J\\}\\\
\widehat{\mathcal{F}}&=\begin{cases}\hat{f}_{R}(\widehat{u}_{RE})&=\hat{f}_{\geq
0.5}(\widehat{u}_{RE})\\\
\hat{f}_{E}(r,\widehat{u}_{RE})&=\hat{f}_{OR}\begin{cases}\hat{f}_{AND}\left(\hat{f}_{<0.375}(\widehat{u}_{RE}),\hat{f}_{NOT}(r)\right)\\\
\hat{f}_{AND}\left(\hat{f}_{<0.125}(\widehat{u}_{RE}),r\right)\\\
\hat{f}_{\geq 0.875}(\widehat{u}_{RE})\end{cases}\\\
\hat{f}_{J}(r,e,\widehat{u}_{J})&=\hat{f}_{OR}\begin{cases}\hat{f}_{AND}\left(\hat{f}_{NOT}(r),\hat{f}_{<0.75}(\widehat{u}_{J})\right)\\\
\hat{f}_{AND}\left(r,\hat{f}_{NOT}(e),\hat{f}_{<0.25}(\widehat{u}_{J})\right)\\\
\hat{f}_{AND}\left(r,e,\hat{f}_{<0.75}(\widehat{u}_{J})\right)\\\
\end{cases}\end{cases}\\\ P(\widehat{\mathbf{U}})&\text{ defined such that
}\widehat{U}_{RE},\widehat{U}_{J}\sim\operatorname{Unif}(0,1)\end{cases}$

One can verify that the quantities induced by $\widehat{M}$ match those in
Table 5. Interestingly, the functions of $\widehat{M}$ are quite different
from those of $\mathcal{M}^{*}$, and in fact, they will likely not match in
non-identifiable quantities. However, if we evaluate $\widehat{M}$ on our
query of interest $Q$, we see that

  1. 1.

$E_{R=0}=1$ when $\widehat{U}_{RE}<0.375$ or $\widehat{U}_{RE}\geq 0.875$, 0
otherwise;

  2. 2.

When $E_{R=0}=0$, $J_{R=1,E_{R=0}}=1$ when $\widehat{U}_{J}<0.25$, 0
otherwise;

  3. 3.

When $E_{R=0}=1$, $J_{R=1,E_{R=0}}=1$ when $\widehat{U}_{J}<0.75$, 0
otherwise.

Hence, we can compute

$\displaystyle P^{\widehat{M}}(J_{R=1,E_{R=0}}=1)$
$\displaystyle=P(0.375\leq\widehat{U}_{RE}<0.875)P(\widehat{U}_{J}<0.25)+P(\widehat{U}_{RE}<0.375\vee\widehat{U}_{RE}\geq
0.875)P(\widehat{U}_{J}<0.75)$ $\displaystyle=(0.5)(0.25)+(0.5)(0.75)=0.5,$

which matches the symbolic result. Additionally, we see that
$P^{\widehat{M}}(J_{R=0}=1)=P(\widehat{U}_{J}<0.75)=0.75$, so we see that the
result induced by $\widehat{M}$ for $Q$ is indeed $0.5-0.75=-0.25$.

Computing the true value of the query from $\mathcal{M}^{*}$ using values from
Table 4, we see that

$\displaystyle
P^{\mathcal{M}^{*}}(J_{R=1,E_{R=0}}=1)-P^{\mathcal{M}^{*}}(J_{R=0}=1)$
$\displaystyle=p_{0}+p_{3}+p_{5}+p_{6}-p_{0}-p_{2}-p_{4}-p_{6}=p_{3}+p_{5}-p_{2}-p_{4}$
$\displaystyle=\frac{1}{32}+\frac{3}{32}-\frac{3}{32}-\frac{9}{32}=-\frac{8}{32}=-0.25.$

Remarkably, by simply fitting the $\mathcal{G}$-NCM on the data and querying
it as if it were $\mathcal{M}^{*}$ itself, we achieve the correct value of $Q$
matching both the true value from $\mathcal{M}^{*}$ and the result from the
symbolic derivation. As long as the training succeeds in matching the data, it
does not matter what functions are learned, and the evaluation requires no
knowledge of the symbolic approach, demonstrating the strength of using the
NCM for evaluating identifiable queries. $\blacksquare$

###  D.4 Discussion on Markovianity

$X$$\mathbf{Z}$$Y$ (a) Example of a graph from a Markovian SCM. Note the lack
of bidirected edges. $Z$ is an observed confounder between $X$ and $Y$.

$X$$\mathbf{Z}$$Y$ (b) If $Z$ is unobserved, then Markovianity no longer holds
because $Z$ becomes an unobserved confounder between $X$ and $Y$.

$X$$Y$ (c) In non-Markovian settings, unobserved confounding between two
observed variables is denoted as a bidirected edge.

Figure 23: Graphs for Markovianity discussion.

One of the major contributions emphasized in this work is the generality of
the approach. In particular, Alg. 1 solves the counterfactual identification
and estimation problems for any counterfactual query with any given input in
any causal graph setting. Notably, we do not make the _Markovianity_
assumption about the generating SCM. Markovianity restricts the set of
unobserved variables $\mathbf{U}$, requiring them to be independent and having
each $U\in\mathbf{U}$ as an input to at most one function in $\mathcal{F}$.
This translates to the lack of bidirected edges in any causal diagram
generated by a Markovian SCM (Fig. 23(a), for example). In other words, since
no unobserved variable can be a common cause of two observed variables,
Markovianity encodes the fact that there is no unobserved confounding.

Markovianity is frequently assumed in many works (including several related
works in Appendix C) because it drastically simplifies the problem. When
attempting to infer $\mathcal{L}_{2}$ queries from $\mathcal{L}_{1}$ data and
the causal diagram $\mathcal{G}$, it turns out that any interventional query
$P(\mathbf{y}_{\mathbf{x}})$ can trivially be identified as

$P(\mathbf{y}_{\mathbf{x}})=\sum_{\mathbf{z}}P(\mathbf{y}\mid\mathbf{z},\mathbf{x})P(\mathbf{z}),$
(38)

where $\mathbf{Z}$ is the set of non-descendants of $\mathbf{X}$; for details,
see (Bareinboim et al., 2022, Corol. 2).

Several observations follow from this observation. First, this implies that
works that attempt to infer $\mathcal{L}_{2}$ queries from $\mathcal{L}_{1}$
data under the Markovianity assumption do not even need to consider the
identification problem, as they are all identifiable given the graph. Kocaoglu
et al. (2018) is an example of a work that can directly estimate the query
without having to be concerned about the identification problem due to the
Markovianity assumption.

Second, when attempting to estimate the $\mathcal{L}_{2}$ query, the estimand
in Eq. 38 is the same for every query, so only one form of the query needs to
be considered when designing an estimator. For example, although not
explicitly assuming Markovianity, the works in Footnote 1 assume conditional
ignorability, which implies the lack of unobserved confounding between the
treatment $X$ and the outcome $Y$. This allows the causal effect
$P(\mathbf{y}_{\mathbf{x}})$ to be identified as stated in Eq. 38, so these
works can focus directly on estimating the specific equation in Eq. 38. On the
counterfactual level, the identification problem may still need to be
considered as there are non-identifiable $\mathcal{L}_{3}$ queries even in the
Markovian case. Still, the problem is drastically simplified as confounding is
no longer an issue. In fact, works like Pawlowski et al. (2020) and Sanchez-
Martin et al. (2021) rely on this simplification for counterfactual estimation
in the Markovian case and do not consider the identification problem despite
potential the lack of identifiability of various queries.

Third, another notable simplification introduced by Markovianity is that
distributions from the SCM can often be factorized into a convenient form. For
example, the observational distribution $P(\mathbf{v})$ can be factorized as
the product:

$P(\mathbf{v})=\prod_{V\in\mathbf{V}}P(v\mid\mathbf{pa}_{V}),$

where $\mathbf{Pa}_{V}$ refer to the set of variables with an edge into $V$ in
$\mathcal{G}$, and $v$ and $\mathbf{pa}_{V}$ are consistent with $\mathbf{v}$.
Similarly, an interventional distribution $P(\mathbf{v}_{\mathbf{x}})$ can be
factorized as

$P(\mathbf{v}_{\mathbf{x}})=\prod_{V\in\mathbf{V}\setminus\mathbf{X}}P(v\mid\mathbf{pa}_{V}),$

where $v$ and $\mathbf{pa}_{V}$ are consistent with $\mathbf{v}$ and
$\mathbf{x}$. Counterfactual distributions follow a similar pattern. For
example, when considering the graph $\mathcal{G}$ in Fig. 23(a),
$P(y_{x},x^{\prime}_{z})$ factorizes as

$P(y_{x},x^{\prime}_{z})=\sum_{z^{\prime}}P(y\mid
x,z^{\prime})P(x^{\prime}\mid z)P(z^{\prime}),$

where $x\neq x^{\prime}$. The pattern is that in every case, each term in the
estimand can be written as the probability of a variable given its parents.
This allows for models to encode the mechanism for each variable $V$ as simply
the conditional distribution $P(V\mid\mathbf{Pa}_{V})$, abstracting away the
functions $\mathcal{F}$ and exogenous noise $P(\mathbf{U})$ from the SCM.
Architectures that benefit from this include DeepSCM from Pawlowski et al.
(2020), which utilizes normalizing flows to model each conditional probability
$P(V\mid\mathbf{Pa}_{V})$, and VACA from Sanchez-Martin et al. (2021), which
fits a variational graph autoencoder to the factorized distributions. In
contrast, in non-Markovian cases, these factorizations typically do not hold.
For example, if the true SCM induces the graph in Fig. 23(c), then
$P(y_{x})\neq P(y\mid x)$ despite $X$ being the only observed parent of $Y$.
Works that rely on Markovian factorizations often do not have a
straightforward generalization to non-Markovian cases.

Fourth, the works cited in this section that leverage the Markovianity
assumption have still provided important contributions to the understanding
and practical approaches to causal and counterfactual inference, and they are
widely praised for their successes. However, Markovianity is a very
restricting assumption and is unrealistic in many practical settings. For
example, consider a simple study on the effects of exercise $X$ on cholesterol
$Y$, with age $Z$ as a confounder. The relationship between variables can be
modeled as shown in Fig. 23(a). If, in the study, data on age is simply not
collected, then $Z$ becomes an unobserved confounder (see Fig. 23(b)), and
Markovianity is immediately violated. The corresponding graph would look like
Fig. 23(c) with $Z$ omitted. In many cases, collecting the data may not even
be an option. For example, if $X$ is smoking, $Y$ is cancer, and $Z$ is the
expressiveness of a gene, perhaps collecting data on gene expressivity is too
difficult or even impossible. In general, one cannot expect to be able to
collect data on every relevant variable in the study, so unobserved
confounding is to be expected in realistic settings.

The solution presented in this paper is therefore important to understand the
more complex, realistic, and general setting of potentially having unobserved
confounding. Since the lack of the Markovianity assumption greatly complicates
the problem of counterfactual inference, many of the theoretical results in
this work (e.g. Thms. 2 and 3) required more involved proofs. The result is
that using Alg. 1, $\mathcal{G}$-NCMs can solve any instance of counterfactual
identification and estimation even under unobserved confounding. This is
demonstrated in the experiments, which consider several settings in which
Markovianity is violated (e.g. the second, third, and fourth rows in Fig. 4,
where the unobserved confounding is shown visually in the causal diagrams as
red dashed bidirected arrows).

##  Appendix E Frequently Asked Questions

  1. Q1.

Why do you assume that the causal diagram exists instead of learning it? Is it
reasonable to expect that the causal diagram is available?

Answer. For a detailed discussion on this topic, see Sec. D.1. In summary, the
assumption of the causal diagram is made out of necessity. The Neural Causal
Hierarchy Theorem (N-CHT) (Xia et al., 2021, Corol. 1) states that
distributions from lower layers almost always underdetermine higher layers in
a measure-theoretic sense (Bareinboim et al., 2022, Thm. 1). This means that
despite its universal expressivity, an NCM trained only on lower layer data
(e.g., observational) cannot be reliably used for inferences on higher layers
(e.g., interventional or counterfactual). Therefore, assumptions are necessary
to make inferences across the PCH’s layers. The causal diagram is a flexible
assumption that is used throughout the literature and encodes a qualitative
description of the generating model, which is often much easier to obtain than
the actual mechanisms of the generating SCM (Pearl, 2000; Spirtes et al.,
2000; Peters et al., 2017). Given the N-CHT, any attempts to learn the causal
diagram will be incorrect or not meaningful unless other assumptions are made
in this process (e.g., faithfulness). While other types of assumptions can be
made, all such assumptions have their own benefits and drawbacks. The goal of
this paper is not to decide which set of assumptions is the best, but rather
to provide the toolkit for data analysts to perform the inferences once the
assumptions have already been made.

  2. Q2.

How could the graphical assumption be relaxed?

Answer. Relaxing the graphical assumption is out of the scope of this work
since the goal of this paper is first to establish a principled approach of
using neural models to perform counterfactual inferences under general well-
understood assumptions. Still, finding a solution to the problem when only
partial information about the graph is available is an interesting direction
for future work. One such possible approach, discussed at the end of Sec. D.1,
is to utilize equivalence classes of causal diagrams, learnable from data, to
perform inferences (Jaber et al., 2018; 2019; 2022, forthcoming). However,
this approach has its own drawbacks, notably, that equivalence classes are not
as descriptive as the true causal diagram, and therefore fewer queries will be
identifiable. Still, this is one approach that can be considered in cases
where assuming the availability of the full causal diagram is too stringent of
a condition.

  3. Q3.

Why do you assume non-Markovianity? Is this more general than solutions that
assume Markovianity?

Answer. That is an interesting question. We do not assume “non-Markovianity”,
but rather we do not make any assumption on Markovianity at all. Typically,
Markovianity is the assumption and not the other way around. Markovianity
requires that there exists no unobserved confounding, but works that do not
assume Markovianity do not enforce such a requirement and are therefore more
general because they work in both Markovian and non-Markovian cases. As
mentioned earlier in the manuscript, it is hard to ascertain Markovianity in
many practical settings, so not making the assumption is more reasonable in
many natural settings. Unlike many related works like those in Appendix C,
this paper considers cases with unobserved confounding, which is represented
graphically through dashed bidirected edges. This approach is more general in
the sense that Alg. 1 can be used even in the presence of unobserved
confounders, which is one of the major contributions of this paper.
Markovianity is often assumed to simplify the problem setting, but it is often
unrealistic to expect that data on all potential confounding variables are
collected in a study. For a detailed discussion on this topic, see Appendix
D.4.

  4. Q4.

Why do you use the $\mathcal{G}$-NCM definition in Def. 2 with feedforward
neural networks and uniform noise? Could the results in this paper also hold
for other function classes and noise?

Answer. The architectural choices made in Def. 2 are standard choices used to
ground the discussion to a more concrete model. The same architecture is also
used in the experiments. Still, the specific choice of architecture is not
important to the theoretical results as long as key properties persist (e.g.
$\mathcal{G}^{(\mathcal{L}_{3})}$-consistency, $\mathcal{L}_{3}$-$\mathcal{G}$
expressiveness), which will often be the case, since the arguments to prove
the theoretical results in this work are mostly independent of the
architecture. In particular, Sec. D.1 of Xia et al. (2021), which introduces
the $\mathcal{G}$-NCM, discusses generalizations of NCMs in terms of functions
and noise which still guarantee the relevant theoretical results, so long as
the functions are sufficiently expressive and the noise variables have well-
defined probability density functions.

  5. Q5.

Given that neural networks are universal approximators, it seems unsurprising
that $\mathcal{G}$-NCMs can perform counterfactual inference. What makes
$\mathcal{G}$-NCMs so special compared to other model classes?

Answer. The structure of this paper follows the delicate balance between
expressiveness and constraint. As emphasized by the Causal Hierarchy Theorem
(Bareinboim et al., 2022, Thm. 1) and (Xia et al., 2021, Corol. 1) ,
expressivity alone is not enough to infer higher layer quantities from lower
layer data. Without assumptions, any model, including NCMs, trained only on
$\mathcal{L}_{1}$ or $\mathcal{L}_{2}$ data will not contain any meaningful
information about $\mathcal{L}_{3}$ quantities. Hence, assumptions in the form
of constraints are discussed in Sec. 2. One important question is how to
properly incorporate such constraints into a model. Constraining a model too
much may hinder its expressiveness, and it may no longer be able to express
the true model. Therefore, a proper model class for counterfactual inference
must be able to incorporate the constraints without losing expressivity within
the constrained class.

The $\mathcal{G}$-NCM model class has both the critical properties of
expressiveness and constraint (Thms. 1 and 2) that are not enjoyed by many
model classes. Examples of these properties are shown in Appendix D.3. Example
4 shows an example of a model class that, although expressive on all three
layers, fails to properly incorporate the constraints of the graph,
$\mathcal{G}$. Examples 5 and 6 show examples of model classes that
incorporate the constraints of $\mathcal{G}$ but fail to achieve
expressiveness within the constrained class. In contrast to these examples,
the $\mathcal{G}$-NCM model class has a natural and intuitive approach to
incorporating the constraints of $\mathcal{G}$ without damaging its
expressiveness. These properties empower the $\mathcal{G}$-NCM the capability
of solving the problems of counterfactual identification and estimation (see
Thm. 3, Corols. 1 and 2, and Alg. 1). Examples 8 and 9 concretely demonstrate
this ability.

When compared to similar works, like the ones in Appendix C, these properties
are not as clearly established. While many of the works do incorporate
graphical constraints, these constraints are often taken for granted as they
are not formally discussed. In contrast, the ideas of graphical constraints
are developed formally in this work in Sec. 2 and Appendix A.1. Many models in
the related works assume Markovianity, which fails in terms of expressiveness
for the same reason as Example 6. Some of the works only focus on the
interventional ($\mathcal{L}_{2}$) level rather than the counterfactual
($\mathcal{L}_{3}$) level.

  6. Q6.

Why would one prefer to identify counterfactual quantities using NCMs? Would
it not be easier first to identify the query using symbolic approaches, then
use the NCM for estimation?

Answer. For a detailed discussion on the difference between symbolic and
optimization approaches, we encourage the reader to read Appendix C.4 of Xia
et al. (2021). In summary, the approach in this paper utilizes
$\mathcal{G}$-NCMs to serve as a proxy model for the true generating SCM. In
contrast, existing symbolic approaches like Correa et al. (2021) work directly
with the graph and data, which is on a level of abstraction without the
generating SCM. These works typically derive an exact answer for the target
query by directly using the constraints of the graph. Indeed, one could use a
hybrid approach of using a symbolic algorithm to identify a counterfactual
query followed by training a $\mathcal{G}$-NCM to obtain an estimate. This is
certainly an option, and the goal of this paper is to provide another solution
route, not to decide for the data analyst which approach is better.

The proxy model approach presented in this paper is insightful to the problem
of counterfactual inference because it shows that $\mathcal{G}$-NCMs are
capable of jointly solving both counterfactual identification and estimation,
which typically appear together in practice. Solving both problems together
indicates that the $\mathcal{G}$-NCM is capable of both reasoning whether the
query has any counterfactual meaning through identification and providing the
practical benefits of statistical modeling through estimation. There are many
reasons why one may prefer a proxy-based approach. First, the generative
properties of the learned proxy model allow it to serve as an infinite source
of samples for the user. Second, the trained model could provide quick
estimation results for multiple queries without the need for retraining,
provided the queries are identifiable. Third, the idea of minimizing and
maximizing the query in Alg. 1 is central to identification problems. Hence, a
similar procedure can be used in similar identification problem settings
without an existing symbolic solution. Fourth, this idea can potentially be
used to bound non-identifiable queries, which is an exciting avenue for future
work.

  7. Q7.

Why would one want to learn an entire SCM to solve the estimation problem for
a specific query instead of just directly estimating the query?

Answer. Indeed, if the user simply wants to estimate one query and is
knowledgeable about the estimand, it may be more efficient, statistically or
computationally, to focus efforts on just estimating that query. However,
there are many reasons why the proxy model approach used in this work may be
beneficial. First, by learning an entire SCM and estimating queries using Eq.
4, one can apply the same procedure to any potential query, reducing the
complexity of the solution. This relieves the user of the need to derive an
estimand. Even if the estimand is known, developing a specific approach to
estimating this estimand can be difficult depending on the form of the
estimand. Second, if the user wants to estimate multiple queries, they could
potentially reuse the same $\mathcal{G}$-NCM without retraining, provided that
the queries are identifiable. Third, as a generative model, the
$\mathcal{G}$-NCM can be used as an infinite source of samples should the user
desire this in addition to having an estimate for the query.

  8. Q8.

Why is the GAN architecture used for the experiments?

Answer. We found that using a GAN approach was the most straightforward way to
implement the $\mathcal{G}$-NCM in practice. When using an implicit generative
model approach, one can construct a $\mathcal{G}$-NCM exactly as described in
Def. 2 and train it as specified in Alg. 3. Using a discriminator network as
the divergence function ($\mathbb{D}_{P}$) provides a flexible and practical
approach to generating a loss for optimization, as understood in GAN
literature. Alternative approaches to modeling the divergence functions all
have practical drawbacks. For example, KL divergence is typically intractable
to compute directly. Appendix B.6.1 explains why an explicit generative model
approach like the MLE-NCM (used in Xia et al. (2021)) is significantly more
complicated. We also include experiments comparing against the MLE-NCM, which
demonstrated some key flaws of the MLE approach, such as poor runtime
scalability (see Fig. 5) or having a bias to output “ID” even in non-ID
settings (see Fig. 4). For further results, see Appendix B.1. Still, the key
theoretical results in this work do not depend on the specific architecture
used in practice, and the choice of using GANs is unimportant in this
discussion. Users can use any architecture they see fit while still applying
the results in this work.

  9. Q9.

Given the training instability of GANs, is the GAN-NCM approach guaranteed to
work?

Answer. Indeed, the theoretical results described in Sec. 3 assume that
infinite data is provided and optimization is perfect, which is unrealistic in
practice. Hence, a hypothesis testing procedure is required to decide
identifiability, described in Appendix B.4. This means that the results from
solving the identification and estimation problems will have some degree of
error in practice. While deriving bounds on the error could be an interesting
direction for future work, the error is shown to be quite manageable
empirically using GANs, as shown in the experiments in Sec. 5. In the
estimation experiments, the error is shown to decrease as samples increase,
demonstrating a level of statistical consistency.

  10. Q10.

Does the work apply to cases with continuous variables?

Answer. The theoretical results in this work assume that the variable domains
are finite, and no claims are made about continuous variables. If extending to
the continuous case, on the theoretical side, counterfactual distributions and
queries would have to be properly defined and new machinery based on measure
theory evoked. It is non-trivial to prove some of the theoretical results of
this paper in continuous settings. While the exogenous variables and
structural functions of the true model can take any form, this paper relies on
the discreteness of the endogenous variables to ensure that there exists a
$\mathcal{G}$-NCM that can exactly match the true model (see also Zhang et al.
(2022)). It is currently an open problem to generalize these results to the
continuous case and in generality.

On the experimental side, sampling or maximizing/minimizing the query from a
$\mathcal{G}$-NCM may have complications depending on the definition of a
counterfactual query/distribution in a continuous setting. Intuitively,
however, if the query can be sampled via Alg. 2, then similar techniques
applied in our paper can be extended without much change (although currently
without theoretical guarantees). For example, if the query is $E[Y_{X=x}|X\geq
0]$ with both continuous $X$ and $Y$, one could (1) sample several instances
of $X$, (2) filter out negative samples, (3) sample $Y$ using the values of
$\mathbf{U}$ from the remaining instances, and (4) compute the average out of
these values of $Y$. Then, this Monte Carlo estimate can be used in Alg. 3 for
maximizing/minimizing the query and estimating the final result if
identifiable. If the query is not so easily compatible with Alg. 2, then the
approach may need to be adjusted to accommodate continuous definitions.

  11. Q11.

Why are there no empirical studies comparing other neural estimation methods
(such as the ones in Appendix C)?

Answer. Despite their notable contributions, as discussed earlier, related
works like the one in Appendix C are not as general as the approach described
in this work. For example, some works assume Markovianity, some only solve
queries from $\mathcal{L}_{2}$ or specific queries from $\mathcal{L}_{3}$, and
some do not consider identifiability (see Appendix C for a detailed discussion
on each of these works). The works in Footnote 1 assume conditional
ignorability or backdoor admissibility, which only hold in a few graphs, such
as graph (a) from Fig. 11 (it does not hold in any other graph in Fig. 11).
The approach presented in this paper is the first work to use neural models to
identify any $\mathcal{L}_{3}$ query given any combination of
$\mathcal{L}_{1}$ and $\mathcal{L}_{2}$ data and any graph setting (even non-
Markovian cases). It will also perform the estimation for identifiable cases,
but it will return “FAIL” instead of estimating non-identifiable cases. Hence,
comparisons with these specific related works are not appropriate since they
cannot handle such general settings.

  12. Q12.

How does the proposed approach scale to higher dimensional variables or with a
greater number of variables?

Answer. The training of the GAN-NCM scales tractably with higher
dimensionality or variables, as shown in Fig. 5. While having higher
dimensions or more variables may require a larger architecture, and therefore
more training time, typically this also comes with a higher data requirement.
Hence, the architecture can be scaled with the complexity of the problem, and
although complex problems will require more training time, the user may be
willing to spend more computational resources on a more complex problem. We
hope this work can be a baseline and act as a catalyst to the study of
different neural architectures with different computational and statistical
properties.