Deep End-to-end Causal Inference

Introduction

Causal-aware decision making is pivotal in many fields such as economics and healthcare . For example, in healthcare, caregivers may wish to understand the effectiveness of different treatments given only historical data. They aspire to estimate treatment effects from observational data, with incomplete or no knowledge of the causal relationships between variables. This is the end-to-end causal inference problem, displayed in Figure 1, where we discover the causal graph and estimate treatment effects together using weaker causal assumptions and observational data.

It is well known that any causal conclusion drawn from observational data requires assumptions that are not testable in the observational environment . Existing methods for estimating causal quantities from data, which we refer to as causal inference methods, commonly assume complete a priori knowledge of the causal graph. This is rarely available in real-world applications, especially when many variables are involved. On the other hand, existing causal graph discovery methods, i.e. those that seek to infer the causal graph from observational data, require assumptions about statistical properties of the data, which often require less human input . These methods often return a large set of plausible graphs, as shown in Figure 1. This incompatibility of assumptions and inputs/outputs makes the task of answering causal queries in an end-to-end manner non-trivial.

We tackle the problem of end-to-end causal inference (ECI) in a non-linear additive noise structural equation model (SEM) with no latent confounders. Our framework aims to allow practitioners to estimate causal quantities using only observational data as input. Our contributions are:

A deep learning-based end-to-end causal inference framework named DECI, which performs both causal discovery and inference. DECI is an autoregressive-flow based non-linear additive noise SEM capable of learning complex nonlinear relationships between variables and non-Gaussian exogenous noise distributions. DECI uses variational inference to learn a posterior distribution over causal graphs. Additionally, we show how the functions learnt by DECI can later be used for simulation-based estimation of (C)ATE. DECI is trained once on observational data; different causal quantities can then be efficiently extracted from the fitted structural equation model.

Theoretical analysis of DECI. We show that, under correct model specification, DECI asymptotically recovers the true causal graph and data generating process. Furthermore, we show that DECI generalizes a number of causal discovery methods, such as Notears , Grandag , and others , providing a unified view of functional causal discovery methods.

Extending DECI for applicability to real data. To make DECI applicable to real-data, we implement support for mixed type (continuous and categorical) variables and missing value imputation.

Insights into ECI performance with more than 1000 experiments. We systematically evaluate DECI, along with a range of combinations of existing discovery and inference algorithms. DECI performs very competitively with baselines from both the causal discovery and inference domains.

Related Work and Preliminaries

Related Work. Our work relates to both causal discovery and causal inference research. Approaches for causal discovery from observational data can be classified into three groups: constraint-based, score-based, and functional causal models . Recently, Zheng et al. framed the directed acyclic graph (DAG) structure learning problem as a continuous optimisation task. Extensions employ nonlinear function approximators, like neural networks, to model the relationships among connected variables. Our work combines this class of approaches with standard causal assumptions to obtain our main theorem about causal graph learning. We extend functional methods to handle mixed data types and missing values. Outside of functional causal discovery, functional relationships between variables (see Figure 1(3)) are typically not learned by discovery algorithms . Thus, distinct models, with potentially incompatible assumptions or inputs, must be relied upon for causal inference. However, when a DAG cannot be fully identified given the available data, constraint and score-based methods often return partially directed acyclic graphs (PAGs) or completed partially directed acyclic graphs (CPDAGs) . Instead of returning a summary graph representing a set, DECI returns a distribution over DAGs in such situation.

Causal inference methods assume that either the graph structure is provided or relevant structural assumptions are provided without the graph . Causal inference can be decomposed into two steps: identification and estimation. Identification focuses on converting the causal estimand (e.g. $P(Y|\text{do}(X=x),W)$) into an estimand that can be estimated using the observed data distribution (e.g. $P(Y|X,W)$). Common examples of identification methods include the back-door and front-door criteria , and instrumental variables . Causal estimation computes the identified estimand using statistical methods, such as simple conditioning, inverse propensity weighting , or matching . Machine learning-based estimators for CATE have also been proposed . Recent efforts to weaken structural assumption requirements allow for PAGs and CPDAGs. Our work takes steps in this direction, allowing inference with distributions over graphs.

Structural Equation Models (SEM). Let ${\mathbf{x}}=(x_{1},\dots,x_{D})$ be a collection of random variables. SEMs model causal relationships between the individual variables $x_{i}$. Given a DAG $G$ on nodes $\{1,\dots,D\}$, ${\mathbf{x}}$ can be described by $x_{i}=F_{i}\left({\mathbf{x}}_{\text{pa}(i;G)},z_{i}\right)$, where $z_{i}$ is an exogenous noise variable that is independent of all other variables in the model, $\text{pa}(i;G)$ is the set of parents of node $i$ in $G$, and $F_{i}$ specifies how variable $x_{i}$ depends on its parents and the noise $z_{i}$. In this paper, we focus on additive noise SEMs, also referred to as additive noise models (ANM), i.e.

Average Treatment Effects. The ATE and CATE quantities allow us to estimate the impact of our actions (treatments) . Assume that ${\mathbf{x}}_{T}$ (with $T\subset\{1,\dots,D\}$) are the treatment variables; the interventional distribution is denoted $p({\mathbf{x}}\mid\text{do}({\mathbf{x}}_{T}={\mathbf{a}}))$. The ATE and CATE on targets ${\mathbf{x}}_{Y}$ for treatment ${\mathbf{x}}_{T}{=}{\mathbf{a}}$ given a reference ${\mathbf{x}}_{T}{=}{\mathbf{b}}$, and conditional on ${\mathbf{x}}_{C}{=}{\mathbf{c}}$ for CATE, are given by

We consider the common scenario where the conditioning variables are not caused by the treatment.

DECI: Deep End-to-end Causal Inference

We introduce DECI, an end-to-end deep learning-based causal inference framework. DECI learns a distribution over causal graphs from observational data and (subsequently) estimates causal quantities. Section 3.1, describes our autoregressive flow based ANM SEM. Section 3.2 lays out the conditions under which DECI will recover the true causal graph given enough observational data (Theorem 1). Section 3.3 shows how the generative model learnt by DECI can be used to simulate samples from intervened distributions, allowing for treatment effect estimation. Section 3.4 extend’s DECI’s real-world applicability by adding support for non-Gaussian exogenous noise, mixed type data (continuous and discrete), and imputation for partially observed data.

DECI takes a Bayesian approach to causal discovery . We model the causal graph $G$ jointly with the observations ${\mathbf{x}}^{1},\ldots,{\mathbf{x}}^{N}$ as

We aim to fit $\theta$, the parameters of our non-linear ANM, using observational data. Once this model is fit, the posterior $p_{\theta}(G|{\mathbf{x}}^{1},\ldots,{\mathbf{x}}^{N})$ characterizes our beliefs about the causal structure.

Prior over Graphs. The graph prior $p(G)$ should characterize the graph as a DAG. We implement this by leveraging the continuous DAG penalty from Zheng et al. ,

which is non-negative and zero only if $G$ is a DAG. We then implement the prior as

where we weight the DAG penalty by $\alpha$ and $\rho$. These are gradually increased during training following an augmented Lagrangian scheme, ensuring only DAGs remain at convergence. We introduce prior knowledge about graph sparseness by penalising the norm $\lambda_{s}\|G\|_{F}$, with $\lambda_{s}$ a scalar.

Likelihood of Structural Equation Model. Following Khemakhem et al. , we factorise the observational likelihood $p_{\theta}({\mathbf{x}}^{n}|G)$ in an autoregressive manner. Rearranging the ANM assumption eq. 1, we have ${\mathbf{z}}=g_{G}({\mathbf{x}};\theta)={\mathbf{x}}-f_{G}({\mathbf{x}};\theta)$. The components of ${\mathbf{z}}$ are independent. If we have a distribution $p_{z_{i}}$ for component $z_{i}$, then we can write the observational likelihood as

Exogenous Noise Model $p_{{\mathbf{z}}}$. We consider two possible models for the distribution of ${\mathbf{z}}$. 1) A simple Gaussian $p_{z_{i}}(\cdot)=\mathcal{N}\left(\cdot|0,\sigma_{i}^{2}\right)$, where per-variable variances $\sigma_{i}^{2}$ are learnt. 2) A flow

We choose the learnable bijections $\kappa_{i}$ to be a rational quadratic splines , parametrised independently across dimensions. We do not couple across dimensions since our SEM requires independent noise variables. Spline flows are significantly more flexible than the Gaussian distributions employed in previous work .

Optimization and Inference Details. The model described presents two challenges. First, the true posterior over $G$ is intractable. Second, maximum likelihood cannot be used to fit the model parameters, due to the presence of the latent variable $G$. We simultaneously overcome both of these challenges using variational inference . We define a variational distribution $q_{\phi}(G)$ to approximate the intractable posterior $p_{\theta}(G|{\mathbf{x}}^{1},\ldots,{\mathbf{x}}^{N})$, and use it to build the ELBO, given by

where $H(q_{\phi})$ represents the entropy of the distribution $q_{\phi}$ and $p_{\theta}({\mathbf{x}}^{n}|G)$ takes the form of eq. 7 (derivation in Section B.3). We choose $q_{\phi}(G)$ to be the product of independent Bernoulli distributions for each potential directed edge in $G$. We parametrize edge existence and edge orientation separately, using the ENCO parametrization . The SEM parameters $\theta$ and variational parameters $\phi$ are trained by maximizing the ELBO. The Gumbel-softmax trick is used to stochastically estimate the gradients with respect to $\phi$. Section B.1 details the full optimisation procedure.

Unified View of functional Causal Discovery. We note that, like DECI, many functional causal discovery methods can be seen from a probabilistic perspective as fitting an autoregressive flow (with a hard acyclicity constraint) for different choices for the exogenous noise distribution $p_{\mathbf{z}}$ and transformation function $f$. We expand on the details of this perspective and formalising it in Appendix C. DECI employs NNs for $f$ and flexible, potentially non-Gaussian, distributions for $p_{\mathbf{z}}$, making it the most flexible member of this family.

2 Theoretical Considerations for DECI

We now show that maximizing the ELBO from eq. 10 recovers both the ground truth data generating process $p(\bm{x};G^{0})$ and true causal graph $G^{0}$ in the infinite data limit. This is formalized in Theorem 1. The assumptions required by the theorem, which are common in causal discovery research, can be informally summarized as (formal assumptions in Appendix A):

Minimality and Structural Identifiability, satisfied by a continuous non-linear ANM ,

Correct Specification, there exists $\theta^{*}$ such that $p_{\theta^{*}}({\mathbf{x}};G^{0})$ matches the data-generating process,

Causal Sufficiency, there are no latent confounders,

Under assumptions 1-5 (Appendix A), the solution $(\theta^{\prime},q^{\prime}_{\phi}({G}))$ from maximizing the ELBO (eq. 10) satisfies $q^{\prime}_{\phi}({G})=\delta({G}={G}^{\prime})$ where ${G}^{\prime}$ is a unique graph. In particular, ${G}^{\prime}={G}^{0}$ and $p_{\theta^{\prime}}(\mathbf{x};{G}^{\prime})=p(\mathbf{x};{G}^{0})$.

The proof is in Appendix A. It consists of two key steps: (i) the maximum likelihood estimate (MLE) of $(\theta,G)$ recovers $p(\bm{x};G^{0})$, (ii) solutions from maximizing the ELBO approach the MLE in the large data limit. Specifically, we show that DECI induces the same joint likelihood as the ground truth and the posterior $q_{\phi}({G})$ is a delta function $\delta(G=G^{0})$ concentrated on the true graph $G^{0}$.

3 Estimating Causal Quantities

We now show how the generative model learnt by DECI can be used to evaluate expectations under interventional distributions, and thus estimate ATE and CATE. As explained above, DECI returns $q_{\phi}(G)$, an approximation of the posterior over graphs given observational data. Then, interventional distributions and treatment effects can be obtained by marginalizing over graphs as

This can be seen as a probabilistic relaxation of traditional causal quantity estimators. When have observed enough data to be certain about the causal graph, i.e. $q_{\phi}(G)=\delta(G-G_{i})$, our procedure matches traditional causal inference. We go on to discuss how DECI estimates (C)ATE.

Estimating ATE. After training, we can use the model learnt by DECI to simulate new samples ${\mathbf{x}}$ from $p_{\theta}({\mathbf{x}}|G)$. We sample a graph $G\sim q_{\phi}(G)$ and a set of exogenous noise variables ${\mathbf{z}}\sim p_{{\mathbf{z}}}$. We then input this noise into the learnt DECI structural equation model to simulate ${\mathbf{x}}$, by applying eq. 1 and eq. 8 on ${\mathbf{z}}$ in the topological order defined by $G$. However, ATE estimation requires samples from the interventional distribution $p\left({\mathbf{x}}_{\setminus T}|\text{do}({\mathbf{x}}_{T}={\mathbf{b}}),G\right)$. These can be obtained by noting that

where $G_{\text{do}({\mathbf{x}}_{T})}$ is the “mutilated” graph obtained by removing incoming edges to ${\mathbf{x}}_{T}$. Thus, samples from this distribution can be obtained by following the sampling procedure explained above, but fixing the values ${\mathbf{x}}_{T}={\mathbf{b}}$ and using $G_{\text{do}({\mathbf{x}}_{T})}$ instead of $G$. Finally, we use these samples to obtain a Monte Carlo estimate of the expectations required for ATE computation eq. 2. Figure 2 illustrates that these samples are from a mixture distribution when the posterior has not collapsed to one graph.

Estimating CATE. We focus on CATE estimation for which the treatment ${\mathbf{x}}_{T}$ is not the cause of the conditioning set ${\mathbf{x}}_{C}$, i.e. there is no directed path from $T$ to $C$ in $G$. Under this assumption, we can estimate CATE by sampling from the interventional distribution $p\left({\mathbf{x}}_{\setminus T}|\text{do}({\mathbf{x}}_{T}),G\right)$ and then estimating the conditional distribution of ${\mathbf{x}}_{Y}$ given ${\mathbf{x}}_{C}$. To make this precise, we let $Y=X\setminus{(T\cup C)}$ denote all variables that we do not intervene or condition on. Conditional densities

making $g_{G}$ approximate the conditional mean of ${\mathbf{x}}_{Y}$. We choose $g_{G}$ to be a basis-function linear model with random Fourier basis functions . As illustrated in Figure 3, we train two separate surrogate models, one for our intervention ${\mathbf{x}}_{T}={\mathbf{a}}$ and one for the reference ${\mathbf{x}}_{T}={\mathbf{b}}$. We estimate CATE as the difference between their outputs evaluated at ${\mathbf{x}}_{C}={\mathbf{c}}$. This process is repeated for multiple posterior graphs samples $G\sim q_{\phi}$, allowing us to marginalise the posterior graphs

General ECI Framework. The probabilistic treatment of the DAG, and the re-use of functional causal discovery generative models for simulation-based causal inference are principles that can be applied beyond DECI. Constraint-based and score-based discovery methods often output a set of DAGs compatible with the data, i.e. a PAG or CPDAG. It is natural to interpret these equivalence classes as uniform distributions over members of sets of graphs. We can then use eq. 11 to estimate causal quantities by marginalizing over these distributions. The quantities inside the expectations over graphs can be estimated using any existing causal inference method, such as linear regression , Double ML , etc. Our experiments explore combinations of discovery methods that return graph equivalence classes with standard causal inference methods. We take expectations over causal graphs since these return the quantity that minimises the posterior expected squared error in our (C)ATE estimates while noting that the best statistic will be application dependent.

4 DECI for Real-world Heterogeneous Data

We extend DECI to handle mixed-type (continuous and discrete) data and data with missing values, which often arise in real-world applications.

Handling Mixed-type Data. For discrete-valued variables, we remove the additive noise structure and directly parameterise parent-conditional class probabilities

where $P_{i}\left({\mathbf{x}}_{\text{pa}(i;G)};\theta\right)$ is a normalised probability mass vector over the number of classes of $x_{i}$, obtained by applying the softmax operator to $f_{i}({\mathbf{x}}_{\text{pa}(i;G)})$. This means that for discrete variables, the output of $f_{i}$ is a vector of length equal to the number of classes for variable $i$. This approach gives a valid likelihood for $p_{\theta}({\mathbf{x}}^{n}|G)$ which we use to train DECI. However, since the full generative model is no longer an ANM, we cannot guarantee that Theorem 1 applies in this setting.

Handling Missing Data. We propose an extension of DECI to partially observed data.We assume that values are missing (completely) at random, the most common setting . We use ${\mathbf{x}}^{n}_{o}$ to denote the observed components of ${\mathbf{x}}^{n}$, ${\mathbf{x}}^{n}_{u}$ to denote the unobserved components, and their joint density in the observational environment is $p_{\theta}({\mathbf{x}}^{n}_{o},{\mathbf{x}}^{n}_{u})$. We approximate the posterior $p(G,{\mathbf{x}}^{n}_{u}|{\mathbf{x}}^{n}_{o})$ with the variational distribution,

which yields the following learning objective

We parameterize the Gaussian imputation distribution $q_{\psi_{n}}({\mathbf{x}}^{n}_{u}|{\mathbf{x}}^{n}_{o})$ using an amortization network , whose input is ${\mathbf{x}}^{n}_{o}$, and output the mean and variance of the imputation distribution $q_{\psi}({\mathbf{x}}^{n}_{u}|{\mathbf{x}}^{n}_{o})$.

Experiments

We evaluate DECI on both causal discovery and causal inference tasks. A full list of results and details of the experimental set-up are in Appendices B and E. Our code is in the supplement.

Datasets. We consider synthetic, pseudo-real, and real data. For the synthetic data, we follow Lachapelle et al. and Zheng et al. by sampling a DAG from two different random graph models, Erdős-Rényi (ER) and scale-free (SF), and simulating each ANM $x_{i}=f_{i}({\mathbf{x}}_{\text{pa}(i;G)})+z_{i}$, where $f_{i}$ is a nonlinear function (randomly sampled spline). We consider two noise distributions for $z_{i}$, a standard Gaussian and a more complex one obtained by transforming samples from a standard Gaussian with an MLP with random weights. We consider number of nodes $d\in\{16,64\}$ with number of edges $e\in\{d,4d\}$. The resulting datasets are identified as ER$(d,e)$ and SF$(d,e)$. All datasets have $n{=}5000$ training samples.

For the pseudo-real data we consider the SynTReN generator , which creates synthetic transcriptional regulatory networks and produces simulated gene expression data that mimics experimental data. We use the datasets generated by ($d{=}20$), and take $n{=}400$ for training. Finally, for the real dataset, we use the protein measurements in human cells from Sachs et al. . We use a training set with $n{=}800$ observational samples and $d{=}11$.

Baselines. We run DECI using two models for exogenous noise: a Gaussian with learnable variance (identified as DECI-G) and a spline flow (DECI-S). We compare against PC , (linear) Notears , the nonlinear variants Notears-MLP and Notears-Sob , Grandag , and ICALiNGAM . When a CPDAG is the output, e.g., from PC, we treat all possible DAGs under the CPDAG as having the same probability. All baselines are implemented with the gcastle package .

Causality Metrics. We report F1 scores for adjacency, orientation and causal accuracy . For DECI, we report the expected values of these metrics estimated over the graph posterior.

Figure 4 shows the results for the data generated with non-Gaussian noise. We observe that DECI achieves the best results across all metrics. Additionally, using the flexible spline model for the exogenous noise (DECI-S) yields better results than the Gaussian model (DECI-G). This is expected, as the noise used to generate the data is non-Gaussian. For Gaussian noise (see Figure 7), both DECI-S and DECI-G perform similarly. Moreover, when data are partially observed (PO), the strong performance of DECI remains, showing that DECI can handle missing data efficiently.

2 End-to-end Causal Inference

We evaluate the end-to-end pipeline, taking in observational data and returning (C)ATE estimates.

Datasets. We generate ground-truth treatment effects to compare against for the ER and SF synthetic graphs that were described in Section 4.1 by applying random interventions on these synthetic SEMs, ensuring at most 3 edges between the intervention and effect variables. For more detailed analysis, we hand-craft a suite of synthetic SEMs, which we name CSuite. CSuite datasets elucidate particular features of the model, such as identifiability of the causal graph, correct specification of the SEM, exogenous noise distributions, and size of the optimal adjustment set. We draw conditional samples from CSuite SEMs with HMC, allowing us to evaluate CATE. Finally, we include two semi-synthetic causal inference benchmark datasets for ATE evaluation: Twins (twin birth datasets in the US) and IHDP (Infant Health and Development Program data) . See Appendix D for all experimental details.

Baselines. To thoroughly evaluate end-to-end inference, we consider different ways of combining discovery and inference algorithms. For DECI, we can use a trained model to immediately estimate (C)ATE. We also consider using the learned DECI graph posterior in combination with existing methods for causal inference on a known graph: DoWhy-Linear and DoWhy-Nonlinear which implement linear adjustment and Double Machine Learning (DML) methods for backdoor adjustment respectively. We also pair other discovery methods with DECI and DoWhy treatment effect estimation, namely the PC algorithm as a baseline and the ground truth graph (when available) as a check. We evaluate end-to-end causal inference on all valid combinations that arise from combining discovery methods in {DECI-Gaussian (DGa), DECI-Spline (DSp), PC, and True graph (T)} with causal inference methods in {DECI-Gaussian (DGa), DECI-Spline (DSp), DoWhy-Linear (L), DoWhy-Nonlinear (N)}.

Metrics. We report RMSE between (C)ATE estimates and the ground truth.

Table 1 provides a high-level summary of our results. For each dataset, we estimated the ATE using each combination of methods, computed the RMSE and took the median over random seeds. We then ranked methods for each dataset (with 1 being the best) and aggregated over the 27 datasets. We find that DECI Spline has the overall best (lowest) rank. RMSE scores are in Appendix E.

In Table 2, we present detailed results for six CSuite datasets. Lin. Exp is a two node linear SEM with exponential noise, only DECI Spline can recover the true graph, ATE estimation quality is similar for different estimator once the true graph is found. Nonlin. Gauss is a two node non-linear SEM with Gaussian noise, only DECI can fit the highly non-linear functional relationship, with equal performance between DECI-Gaussian and -Spline. Large backdoor is a larger nonlinear SEM with non-Gaussian noise in which adjusting for all confounders is valid, but of high variance. For DECI-Spline, which performs well on discovery, the ATE estimation is best using DECI, as DoWhy takes the maximal adjustment set thereby increasing estimator variance. Weak arrows is a similar SEM to Large backdoor, except that a maximal adjustment set is now necessary. Here, DECI-Spline is best for discovery, but is somewhat less accurate for ATE estimation given the right graph. Nonlin. Simpson is an adversarially constructed dataset where 1) the true graph is theoretically identifiable (it is a non-linear ANM), but difficult to discover in practice, 2) ATE estimation is very poor given the wrong graph (Simpson’s paradox). All methods perform equally badly. Symprod Simpson is a similar but slightly easier dataset, for which DECI-Spline with DML does well.

We performed similar analysis for ATE estimation on ER and SF datasets, on additional CSuite datasets that contain discrete variables or are not theoretically identifiable, and CATE estimation on a subset of CSuite. See Appendix E.

On the semi-synthetic benchmark datasets, Twins and IHDP, we evaluated both ATE and CATE estimation as shown in Figure 5. For ATE estimation, DECI-Spline is fractionally better than baselines on Twins and significantly better for IHDP. On IHDP, it appears that only DECI-Spline was successful at causal discovery, and given the right graph, DECI-Spline is the best method for computing ATE. For CATE estimation, a similar pattern.

Summary Across all experiments we see that DECI enables end-to-end causal inference with competitive performance on both synthetic and more realistic data. DECI particularly performs well compared to other methods when its ability to handle nonlinear functional relationship and non-Gaussian noise distributions comes into play in causal discovery or causal inference. Other ECI method combinations can achieve strong performance, but have weak performance if either step’s assumptions are violated. We find DECI-Spline particularly attractive given its high degree of flexibility—it generally performs on par with or better than other methods.

Discussion, scope, and limitations

Causal inference requires causal assumptions on the relationships between variables of interest. The field of causal discovery aims to learn about these relationships from observational data, given some non-causal assumptions on the data generating process. Motivated by a real-world application where our knowledge of causal relationships is incomplete, DECI combines ideas from causal discovery and inference to go directly from observations to causal predictions. This formulation requires us to adopt assumptions, namely, that the data is generated with a non-linear ANM and that there are no unobserved confounders. Empirically, we find DECI to perform well when these assumptions are satisfied, validating the viability of an end-to-end approach. However, the non-linear ANM assumptions made by DECI are impossible to check in most real-world scenarios. Thus, combining the output of discovery methods with incomplete causal assumptions is an attractive avenue to make end-to-end methods more robust in the future. Interestingly, even in our experiments where DECI’s assumptions are violated (missing data, discrete type observations, etc), we do not find its performance to degrade severely. This encouraging result motivates us to extend our theoretical analysis to the mixed type and missing data settings in future work.

Acknowledgments and Disclosure of Funding

We would like to thank Vasilis Syrgkanis for insightful discussions regarding causal inference methods and EconML usage; we thank Yordan Zaykov for engineering support; we thank Biwei Huang and Ruibo Tu for feedback that improved this manuscript; we thank Maria Defante, Karen Fassio, Steve Thomas and Dan Truax for insightful discussions on real-world needs which inspired the whole project.

References

Appendix A Theoretical Considerations for DECI

DECI can be categorized as a functional score-based causal discovery approach, which aims to find the model parameters $\theta$ and mean-field posterior $q_{\phi}(G)$ by maximizing the ELBO (eq. 10). A key statistical property of DECI is whether it is capable of recovering the ground truth data generating distribution and true graph $G^{0}$ when DECI is correctly specified and with infinite data. In the following, we will show that DECI is indeed capable of this under standard assumptions. The main idea is to first show that the maximum likelihood estimate (MLE) recovers the ground truth due to the correctly specified model. Then, we prove that optimal solutions from maximizing the ELBO are closely related to the MLE under mild assumptions.

For a distribution $p_{\theta}({\mathbf{x}};{G})$ generated by DECI with graph ${G}$ and parameter $\theta$, we assume the minimality condition holds . Namely, the distribution $p_{\theta}({\mathbf{x}};{G})$ does not satisfy the local Markov condition with respect to any sub-graph of ${G}$.

To satisfy this assumption in practice, one can leverage Proposition 17 from Peters et al. , stating that the minimality condition can be satisfied if the the model is a continuous additive noise model (ANM) and its structural assignments are not a constant with respect to any of its arguments. In practice, one can always add an edge pruning step to remove spurious edges .

We assume that the DECI model satisfies the structural identifiability. Namely, for a distribution $p_{\theta}({\mathbf{x}};{G})$, the graph ${G}$ is said to be structural identifiable from $p_{\theta}({\mathbf{x}};{G})$ if there exists no other distribution $p_{\theta^{\prime}}({\mathbf{x}};{G}^{\prime})$ such that ${G}\neq{G}^{\prime}$ and $p_{\theta}({\mathbf{x}};G)=p_{\theta^{\prime}}({\mathbf{x}};G^{\prime})$.

For general SEM, this assumption does not hold. In fact, one can always search for functions resulting in the independence of cause and mechanisms in both directions . However, by correctly restricting the function family and the form of structural assignments, one can obtain structural identifiability . From the formulation of DECI, it is a special case of the non-linear ANM . Given the non-linear ANM assumption, together with the minimality condition and some additional mild assumptions, Theorem 20 from Peters et al. proves that our DECI model is structural identifiable.

We assume the DECI model is correctly specified. Namely, there exists a parameter $\theta^{*}$ such that $p_{\theta^{*}}({\mathbf{x}};G^{0})=p({\mathbf{x}};G^{0})$.

In practice, this assumption is hard to check in general. However, we can leverage the universal approximation capacity of neural networks , meaning that they can approximate continuous functions arbitrarily well. This flexibility gives us a higher chance that this assumption indeed holds.

We assume DECI and the ground truth are causally sufficient. Namely, there are no latent confounders in the model.

We assume for all parameters $\theta$ and possible graphs $G$, the following holds:

A.2 MLE Recovers Ground Truth

The likelihood has often been used as the score function for causal discovery. For example, Carefl adopts the likelihood ratio test in the bivariate case, which is equivalent to selecting the causal directions with the maximized likelihood. However, they did not explicitly show that the resulting model recovers the ground truth for the multivariate case. In addition, Zhang et al. proved that maximizing likelihood for bivariate causal discovery is equivalent to minimizing the dependence between the cause and the noise variable. With the correctly specified, structural identifiable model, the resulting noise and cause are independent through maximizing the likelihood, indicating the graph is indeed causal. However, it is non-trivial to generalize this to the multivariate case that we treat in DECI. In the following, we will show that under a correctly specified model and with maximum likelihood training with infinite data, DECI can recover the unique ground truth graph $G^{*}=G^{0}$ and the true data generating distribution $p_{\theta^{*}}({\mathbf{x}};G^{*})=p({\mathbf{x}};G^{0})$, where $(\theta^{*},G^{*})$ are MLE solutions.

Assuming assumptions 1–5 hold, we denote $(\theta^{*},G^{*})$ as the MLE solution with infinite training data. Then, we have

The key idea is to show that with arbitrary $(\theta,G)$, we have the following:

where the inequality is due to $\log t\leq t-1$. With assumption 3–4, we know there are no latent confounders and the model is correctly specified. Then, the above equality holds when $(\theta^{*},G^{*})$ induces the same join likelihood $p({\mathbf{x}};G^{0})$. Since the model is structural identifiable, we must have $G^{*}=G^{0}$. ∎

A.3 DECI Recovers the Ground Truth

Assume a variational distribution $q_{\phi}(G)$ over a space of graphs $\mathcal{G}_{\phi}$, where each graph $G\in\mathcal{G}_{\phi}$ has a non-zero associated weight $w_{\phi}(G)$. With the soft prior $p(G)$ defined as eq. 6 and bounded $\lambda,\rho,\alpha$, we have

First, we write down the definition of KL divergence

where $Z$ is the normalizing constant for the soft prior. From the definition and assumptions, it is trivial to know that $\log w_{\phi}(G)$, $\lambda\|G\|_{F}^{2}$ are bounded for all $G\in\mathcal{G}_{\phi}$. In the following, we show that $h(G)$ and $\log Z$ are also bounded.

From the definition of the DAG penalty, we have $h(G)=tr(\exp(G\odot G))-D$. The matrix exponential is defined as

Thus, it is obvious that for any $k$, the number of closed loops with length $k$ must be finite. Hence, it is trivial that $h(G)<\infty$. Therefore, with bounded $\lambda,\rho,\alpha$, the un-normalized soft prior

Thus, the normalizing constant $Z$ must be finite since there are only finite number of possible graphs.

Therefore, these must exists a constant $M_{\phi,G}$ such that $\log w_{\phi}(G)+\lambda\|G\|_{F}^{2}+\rho h(G)^{2}+\alpha h(G)+\log Z<M_{\phi,G}$. Hence, we have

where the third inequality is obtained by using Cauchy-Schwarz inequality. ∎

Now, we can prove that DECI can recover the ground truth. Recalling Theorem 1,

Assuming assumptions 1–5 are satisfied, the solution $(\theta^{\prime},q^{\prime}_{\phi}(G))$ from maximizing ELBO (eq. 10) in the infinite data limit satisfies $q^{\prime}_{\phi}(G)=\delta(G=G^{\prime})$ where $G^{\prime}$ is a unique graph. In particular, we have $G^{\prime}=G^{0}$ and $p_{\theta^{\prime}}({\mathbf{x}};G^{\prime})=p({\mathbf{x}};G^{0})$.

In terms of optimization, it is equivalent to re-write the ELBO (eq. 10) as

Now, under the infinite data limit and the definition of $q_{\phi}$, we have

where the second and third equalities are from Lemma 1 and the law of large numbers, respectively. Let $(\theta^{*},G^{*})$ be the solutions from MLE (Proposition 1). Then, since $\sum_{G\in\mathcal{G}_{\phi}}w_{\phi}(G)=1$, $w_{\phi}(G)>0$, we have

with the equality holding when every graph $G\in\mathcal{G}_{\phi}$ and associated parameter $\theta_{G}$ satisfies

From proposition 1, under correctly specified model, we have

Thus, for a $G^{\prime}\in\mathcal{G}_{\phi}$ and associated parameter $\theta^{\prime}$, the condition in eq. 17 becomes

which implies $p_{\theta^{\prime}}({\mathbf{x}};G^{\prime})=p({\mathbf{x}};G^{0})$. Since DECI is structural identifiable, this means $G^{\prime}=G^{0}$ and it is unique. Thus, the graph space $\mathcal{G}_{\phi}$ only contains one graph $G^{\prime}$, and $q^{\prime}_{\phi}(G)=\delta(G=G^{\prime})$. ∎

One should note that we do not explicitly restrict the noise distribution, indicating it still holds with the spline noise (Section 3.4). However, the above theorem implicitly assumes that DECI is a special case of ANM for structural indentifiability and that the data has no missing values. Thus, it is not applicable for DECI with the mixed-type and missing value extensions. We leave a more general theoretical guarantee to future work.

Appendix B Additional Details for DECI

As mentioned in the main text, we gradually increase the values of $\rho$ and $\alpha$ as optimization proceeds, so that non-DAGs are heavily penalized. Inspired by Notears, we do this with a method that resembles the updates used by the augmented Lagrangian procedure for optimization . The optimization process interleaves two steps: (i) Optimize the objective for fixed values of $\rho$ and $\alpha$ for a certain number of steps; and (ii) Update the values of the penalty parameters $\rho$ and $\alpha$. The whole optimization process involves running the sequence (i)–(ii) until convergence, or until the maximum allowed number of optimization steps is reached.

Step (i). Optimizing the objective for some fixed values of $\rho$ and $\alpha$ using Adam . We optimize the objective for a maximum of $6000$ steps or until convergence, whichever happens first (we stop early if the loss does not improve for $1500$ optimization steps. If so, we move to step (ii)). We use Adam, initialized with a step-size of $0.01$. During training, we reduce the step-size by a factor of $10$ if the training loss does not improve for $500$ steps. We do this a maximum of two times. If we reach the condition a third time, we do not decrease the step-size and assume optimization has converged, and move to step (ii).

B.2 Other Hyperparameters.

We use $\lambda_{s}=5$ in our prior over graphs eq. 6. For ELBO MC gradients we use the Gumbel softmax method with a hard forward pass and a soft backward pass with temperature of $0.25$.

For the non-Gaussian noise model in eq. 9, the bijection $\kappa$ is an 8 bin rational quadratic spline with learnt parameters.

In section 3.3, for ATE estimation we compute expectations by drawing 1000 graphs from DECI’s graph posterior $q_{\phi}$ and for each graph we draw 2 samples of ${\mathbf{x}}_{Y}$ for a total of 2000 samples. For CATE estimation, we need to train a separate surrogate predictor per graph samples. We draw 10 different graph samples and 10000 $({\mathbf{x}}_{C},{\mathbf{x}}_{Y})$ pair samples for each graph. We use these to train the surrogate models.

Our surrogate predictor is a basis function linear model with 3000 random Fourier features drawn such that the model approximates a Gaussian process with a radial basis function kernel of lengthscale equal to $1$ .

B.3 ELBO Derivation

The goal of maximum likelihood involves maximizing the likelihood of the observed variables. For DECI (with fully observed datasets) this corresponds to the log-marginal likelihood

Marginalising $G$ in the equation above is intractable, even for moderately low dimensions, since the number of terms in the sum grows exponentially with the size of $G$ (which grows quadratically with the data dimensionality $D$).

Variational inference proposes to use a distribution $q_{\phi}(G)$ to build the ELBO, a lower bound of the objective from eq. 18, as follows:

B.4 Intervened Density Estimation with DECI

Apart from (C)ATE estimation, DECI may also be used to evaluate densities under intervened distributions. For a given graph, the density of some observation vector ${\mathbf{a}}$ is computed by evaluating the base distribution density after inverting the SEM

noting that the transformation Jacobian is the identity. We then marginalise the graphs using Monte Carlo:

In the rest of this section we derive methods that allow using DECI to estimate causal quantities.

Under $G_{\text{do}({\mathbf{x}}_{T})}$, $i\in T$ correspond to parent nodes and we have the following factorisation: $p({\mathbf{x}}|G_{\text{do}({\mathbf{x}}_{T})})=p({\mathbf{x}}_{\setminus{T}}|G_{\text{do}({\mathbf{x}}_{T})})\prod_{i\in T}p({\mathbf{x}}_{i})$. We can then evaluate the interventional density of an observation ${\mathbf{x}}_{\setminus T}={\mathbf{a}}$ with DECI as

which amounts to evaluating the density of the exogenous noise correspondint to non-intervened variables. We can then marginalise the graph using Monte Carlo as in eq. 26.

B.5 Relationship with Khemakhem et al. [30]

Khemakhem et al. introduced Carefl, a method that uses autoregressive flows to learn causal-aware models, using the variables’ causal ordering to define the autoregressive transformations. The method’s main benefit is its ability to model complex nonlinear relationships between variables. However, Carefl alone is insufficient for causal discovery, as it requires the causal graph structure as an input. The authors propose a two-step approach. First, run a traditional constraint-based method (e.g., PC) to find the graph’s skeleton and orient as many edges as possible, and second, fit several flow models to determine the orientation of the remaining edges. The drawbacks of this approach include the dependence on an external causal discovery methods (which will inherently limit Carefl’s performance to that of the method used), and the cost of fitting multiple flow models to orient the edges that are left unoriented after the first step. Our method extends Khemakhem et al. to learn the causal graph among multiple variables and perform end-to-end causal inference.

B.6 Discussion on Causal Discovery Methods

When performing causal discovery, DECI returns a posterior over graphs. Most other causal discovery methods return either a single graph or an equivalence class of graphs. However, we can re-cast these methods in the probabilistic framework used by DECI by noting that a posterior over graphs takes the form

In this equation, the likelihood measures the degree of compatibility of a certain DAG architecture with the observed data. For score-based discovery methods we take the score to be $\log p({\mathbf{X}}|G)$. For functional discovery methods we use the exogenous variable log-density. Constraint-based methods can also be cast in this light by assuming a uniform distribution over all graphs in their outputted equivalence class $\mathcal{G}$: $\log p({\mathbf{X}}|G)=-\log|\mathcal{G}|,\forall\,G\in\mathcal{G}$. To what degree these methods succeed at constraining the space of possible graphs will depend on how well their respective assumptions are met and the amount of data available .

Appendix C Unified View of Causal Discovery Methods

This section introduces a simple analysis showing that, similarly to DECI, most causal discovery methods based on continuous optimization can be framed from a probabilistic perspective as fitting a flow. The benefits of this unified perspective are twofold. First, it allows a simple comparison between methods, shedding light on the different assumptions used by each one, their benefits and drawbacks. Second, it simplifies the development of new tools to improve these methods, since any improvements to one of them can be easily mapped to the others by framing them in this unified framework (e.g. our extensions to handle missing values and flexible noise distributions can be easily integrated with Notears).

where $h(\cdot)$ is the algebraic characterization of DAGs from eq. 5.

The acyclicity constraint is enforced by constraining the optimization domain to $\Theta=\{\theta:h(W(\theta))=0\}$. Then, the maximum likelihood objective can be written as

where the first equality we use the change of variable formula, valid because the transformation ${\mathbf{z}}=g({\mathbf{x}};\theta)={\mathbf{x}}-f({\mathbf{x}};\theta)$ is invertible for any $\theta\in\Theta$, and the second equality uses that the function $f({\mathbf{x}}^{n};\theta)$ has Jacobian-determinant equal to $1$, due to the constraint $\Theta=\{\theta:h(W(\theta))=0\}$. ∎

Lemma 2 is the main building block in the formulation of continuous optimization-based causal discovery methods from a probabilistic perspective as fitting flow models. This is simply because the objective used by each of these methods can be exactly recovered from eq. 29 with specific choices for $f({\mathbf{x}};\theta)$ and $p_{\mathbf{z}}$.

uses a standard isotropic Gaussian for $p_{\mathbf{z}}$ and a linear transformation for $f({\mathbf{x}},\theta)$. (This is similar to DECI-Gaussian, although DECI permits fully nonlinear functions.)

uses a standard isotropic Gaussian for $p_{\mathbf{z}}$ and $D$ independent multi-layer perceptrons, one for each component of $f({\mathbf{x}},\theta)$.

uses a standard isotropic Gaussian for $p_{\mathbf{z}}$ and a weighted linear combination of nonlinear basis functions.

uses a standard isotropic Gaussian for $p_{\mathbf{z}}$ and a GNN for $f({\mathbf{x}},\theta)$.

uses a factorized Gaussian with mean zero and learnable scales for $p_{\mathbf{z}}$ and $D$ multi layer perceptrons, one for each component of $f({\mathbf{x}},\theta)$.

This is a linear method whose original formulation was already in a probabilistic perspective, using a linear transformation for $f({\mathbf{x}};\theta)$.

In summary, recently proposed causal discovery methods based on continuous optimization can be formulated from a probabilistic perspective as fitting a flow with different constraints, transformations, and base distributions. This unified formulation sheds light on the assumptions done by each method (e.g. a Gaussian noise assumption, either implicitly as in Notears or explicitly as in Grandag) and, more importantly, simplifies the development of new tools to improve them. For instance, the ideas proposed to deal with partially-observed datasets and non-Gaussian noise are readily applicable to any of the causal discovery methods mentioned in this section, addressing some of their limitations .

Appendix D Datasets Details

Our two benchmark datasets are constructed following similar procedures described in Louizos et al. .

IHDP . This dataset contains measurements of both infants (birth weight, head circumference, etc.) and their mother (smoked cigarettes, drank alcohol, took drugs, etc) during real-life data collected in a randomized experiment. The main task is to estimate the effect of home visits by specialists on future cognitive test scores of infants. The outcomes of treatments are simulated artificially as in ; hence the outcomes of both treatments (home visits or not) on each subject are known. Note that for each subject, our models are only exposed to only one of the treatments; the outcomes of the other potential/counterfactual outcomes are hidden from the mode, and are only used for the purpose of ATE/CATE evaluation. To make the task more challenging, additional confoundings are manually introduced by removing a subset (non-white mothers) of the treated children population. In this way we can construct the IHDP dataset of 747 individuals with 6 continuous covariates and 19 binary covariates. We use 10 replicates of different simulations based on setting B (log-linear response surfaces) of , which can downloaded from https://github.com/AMLab-Amsterdam/CEVAE. We use a 70%/30% train-test split ratio. Before training our models, all continuous covariates are normalized.

Ground Truth ATE and CATE Estimation for TWINS and IHDP. In both benchmark datasets, since the held-out hypothetical outcomes of counterfactual treatments are already known, the the ground truth ATE can be naively estimated by averaging the difference between the factual and counterfactual outcomes across the entire dataset. The CATE estimation is a bit tricky, since both datasets contains covariates collected from real-world experiments, in which the underlying ground truth causal graph structure is unknown. As a result, exact CATE estimation is generally impossible for continuous conditioning sets. Therefore, when evaluating the CATE estimation performance on TWINS and IHDP, we focus only on discrete variables (binary and categorical) as conditioning set. This allows unbiased estimation of ground truth CATE by simply averaging the treatment effects on subgroups of subjects in the dataset, that have the corresponding discrete value in the conditioning set. We consider only single conditioning variable at a time, and estimate the corresponding CATE for evaluation.

We develop Causal Suite (CSuite), a number of small to medium (2–12 nodes) synthetic datasets generated from hand-crafted Bayesian networks with the intention of testing different capabilities of causal discovery and inference methods. All continuous-only datasets take the form of additive noise models.

Each dataset comes with a training set of 2000 samples, and between 1 and 2 intervention test sets. Each intervention test set has a treatment variable, treatment value, reference treatment value and effect variable. We estimate the ground truth ATE by drawing 2000 samples from the treated and reference intervened distributions. For the datasets used to evaluate CATE, we generate samples from conditional intervened distributions by using Hamiltonian Monte Carlo. We employ a burn-in of 10k steps and a thinning factor of 5 to generate 2000 conditional samples, which we then use to compute our ground truth CATE estimate. We note that because all ground truth causal quantities are estimated from samples, there is a lower bound on the expected error that can be obtained by our methods. When methods obtain an error equal or lower we say that they have solved the task.

A two node graph (Figure 6(a)) with a linear relationship and Gaussian noise. We have $X_{1}\sim N(0,1)$ and $X_{2}=\tfrac{1}{2}X_{1}+\tfrac{\sqrt{3}}{2}Z_{2}$ where $Z_{2}\sim N(0,1)$ is independent of $X_{1}$. The observational distribution is symmetrical in $X_{1}\leftrightarrow X_{2}$. The graph is not identifiable. The best achievable performance on this dataset is obtained when there is a uniform distribution over edge direction.

linexp

A two node graph (Figure 6(a)) with a linear functional relationship, but with exponentially distributed additive noise. We have $X_{1}\sim N(0,1)$ and $X_{2}=\tfrac{1}{2}X_{1}+\tfrac{\sqrt{3}}{2}(Z_{2}-1)$ where $Z_{2}\sim\text{Exp}(1)$ is independent of $X_{1}$. By using non-Gaussian noise, the graph becomes identifiable. However, the inference problem will be more challenging for methods sensitive to outliers, such as those that assume Gaussian noise.

nonlingauss

A two node graph (Figure 6(a)) with a nonlinear relationship and Gaussian additive noise. We have $X_{1}\sim N(0,1)$ and $X_{2}=\sqrt{6}\exp(-X_{1}^{2})+\alpha Z_{2}$ where $Z_{2}\sim N(0,1)$ is independent of $X_{1}$ and $\alpha^{2}={1-6\left(\frac{1}{\sqrt{5}}-\frac{1}{3}\right)}$. Note $\text{Var}(X_{2})=1$ and $\text{Cov}(X_{1},X_{2})=0$. By having a linear correlation of zero between $X_{1}$ and $X_{2}$, this dataset creates a potential failure mode for causal inference methods that assume linearity.

nonlin_simpson

symprod_simpson

Another Simpson’s paradox using the graph Figure 6(c). This dataset is similar to nonlin_simpson with 2 key differences: 1) the effect variable is the result of a product between the confounding variable and the treatment variable. This makes drawing causal inferences require non-linear function estimation. Additionally, the ATE is close to 0. The conditioning variable for the CATE task is a descendant of the confounding variable. This dataset probes for methods’ capacity to reduce their uncertainty about a confounding variables based on values of its child variables.

large_backdoor

A nine node graph, as shown in Figure 6(d). This dataset is constructed so that there are many possible choices of backdoor adjustment set. While both minimal and maximal adjustment sets can result in a correct solution, the a minimal adjustment set results in a much lower-dimensional adjustment problem and thus will result in lower variance solutions. The conditioning node for the CATE task is a child of the root variable. Thus the CATE task probes for methods’ capacity to infer the value of an observed confounder from one of its children. All variables are continuous, with nonlinear structural equations and non-Gaussian additive noise.

weak_arrows

A nine node graph, as shown in Figure 6(e). Unlike the previous dataset, when the true graph is known, a large adjustment set must be used. The causal discovery challenge revolves around finding all arrows, which are scaled to be relatively weak, but which have significant predictive power for $X_{9}$ in aggregate. This dataset tests methods’ capacity to identify the full adjustment set and adjust for a large number of variables simultaneously.

cat_to_cts

A two node (Figure 6(a)) graph with categorical $X_{1}$ and continuous $X_{2}$ with an additive noise model. We have $X_{1}\sim\text{Cat}\left(\tfrac{1}{4},\tfrac{1}{4},\tfrac{1}{2}\right)$ takes values in $\{0,1,2\}$ and $X_{2}=X_{1}+\tfrac{8}{5}(s(Z_{2})-1)$ where $s(x)=\log(\exp(x)+1)$ is the softplus function, and $Z_{2}\sim N(0,1)$ is independent of $X_{1}$.

cts_to_cat

A two node (Figure 6(a)) graph with continuous $X_{1}$ and categorical $X_{2}$. We take $X_{1}\sim U(-\sqrt{3},\sqrt{3})$ and $X_{2}$ categorical on $\{0,1,2\}$ with the following conditional probabilities

In this problem, we treat $X_{2}$ as the treatment and $X_{1}$ as the target, giving a theoretical ATE of zero.

mixed_simpson

Similar to the nonlin_simpson dataset, using the graph of Figure 6(b), but with $X_{3}$ categorical on three categories, and $X_{1}$ binary.

large_backdoor_binary_t

Similar to the large_backdoor dataset, using the graph of Figure 6(d), but with $X_{8}$ binary.

weak_arrows_binary_t

Similar to the weak_arrows dataset, using the graph of Figure 6(e), but with $X_{8}$ binary.

mixed_confounding

A large, mixed type dataset with 12 variables, as shown in Figure 6(f). In this dataset, $X_{1},X_{5}$ are binary, $X_{3},X_{6},X_{8}$ are categorical on three categories, and other variables are continuous. We utilise nonlinear structural equations and non-Gaussian additive noise.

Appendix E Additional Results

Figure 4 in the main text shows causal discovery results for the case where synthetic data was generated using non-Gaussian noise. In that case it was observed that using DECI together with a flexible noise model performed better than DECI with a Gaussian noise model. Figure 7 shows results for synthetic data generated using Gaussian noise. As expected, in this case using a Gaussian noise model is beneficial, although DECI with a spline noise mode still performs strongly.

E.2 Summary of CSuite Results

Comprehensive results on CSuite ATE and CATE performance are shown in figures 8 and 9. We first provide a summary of results here and then go into per-dataset analysis in the following subsection.

We find DECI to perform consistently well in our 2 node datasets. It learns a uniform posterior over graphs in the non-identifiable setting, it fits non-linear functions well and it is robust to heavy tailed noise when employing the spline noise model. We find linear and non-linear DML inference to also perform acceptably, with the exception of the heavy tailed noise case, where the methods overfit to outliers and thus estimate ATE poorly.

On the larger (4 and 12 node) datasets, when the true graph is available, DECI provides ATE estimates competitive with the well-established non-linear DML method. Notably, DECI outperforms backdoor adjustment methods when the number of possible adjustment sets is large. Choosing the optimal adjustment set is an np-hard problem and the most common approach is to simply choose the largest one. This leads to doWhy suffering from variance. DECI’s simulation-based approach avoids having to choose an adjustment set. On the other hand, for densely connected graphs where the strength of the connection between nodes is low, DECI struggles to capture the funcitonal relationships in the data and DML is most competitive. For CATE estimation DECI provides superior performance in all datasets and is able to completely solve all tasks but one.

When the graph is learnt from the data, the non-linear nature of our (4 and 12 node) datasets together with their heavy tailed noise make the discovery problem very challenging. We find that the PC algorithm provides very poor results or fails to find any causal DAGs compatible with the data when working with these datasets. We find both DECI to provide more acceptable performance with the DECI-spline variant producing more reliable results. In this learnt graph setting, causal inference performance deteriorates sharply as a consequence of imperfect causal discovery. However, our findings in terms of relative performance among inference methods stay the same.

E.3 Discussion of Continuous CSuite Results

lingauss: When the true graph is available, all our causal inference methods are able to solve this problem. However, when the graph needs to be identified from the data, causal discovery accuracy is around 50%. DECI discovery converges to a posterior with half of its mass on the right distribution resulting in DECI inference methods showing the lowest error.

linexp: The non-Gaussian noise causes difficulties for DECI-Gaussian, which identifies the wrong orientation in a majority of cases. As a result, inference algorithm yield poor results. Surprisingly, the PC algorithm is also unable to identify the causal graph, leading to overall poor inference performance. With the spline noise model, DECI successfully identifies the causal graph, allowing for all inference algorithms to solve the problem.

nonlingauss: The non-linear relationship between variables leads all DECI discovery runs to successfully recover the edge direction for this dataset while PC consistently identifies the wrong edge direction. As expected, linear ATE estimation performs poorly on this task. However, we find DoWhy non-linear to not fare much better, likely this is because DML still assumes a linear relationship between treatment and target. DECI solves the task successfully.

nonlin_simpson: Even when the true graph is available, none of our inference methods are able to recover the true ATE on this more difficult task. We observe non-linear methods (DECI and DoWhy-nonlinear) to perform similarly to each other and more strongly than the simple linear adjustment. For the CATE task, the true value is close to 0. This is correctly identified by both DECI-Gaussian and DECI-Spline. Interestingly, we find both linear and non-linear DoWhy variants to overestimate the causal effect when using the backdoor criterion. We attribute this to DECI solving a lower dimensional problem when estimating CATE. While DECI simply regresses the conditioning variable onto the effect variable. The backdoor adjustment employed by DoWhy requires regression from the joint space of conditioning variables and confounders onto the effect variables. The latter procedure involves estimating the relative strength of confounders and conditioning variables, which is a more challenging task.

This dataset provides a challenging causal discovery task. DECI identifies the correct edges with probability 0.9. It capacity to recover the edge orientation is slightly worse 0.65. This imperfect causal discovery leads to poor inference for all methods. (potentially because they get Simpson’s paradox the wrong way around).

symprod_simpson: Even with access to the true graph, no inference method is able to solve this problem. However, we find non-linear methods to clearly outperform linear adjustment for both CATE and ATE estimation. Among non-linear methods, performance is similar for ATE estimation, with DECI-Gaussian performing slightly better than DoWhy-nonlinear and DECI-Spline slightly worse. However, when estimating CATE, DECI inference present an error twice as low as nonlinear DoWhy. Again, we attribute this to the backdoor adjustment employed by DoWhy being a more challenging inference task than the 1d regression on simulated data employed by DECI.

In terms of causal discovery, results are similar to nonlin-simpson with PC failing completely and DECI obtaining an adjacency score of 0.92 and orientation of 0.7. The imperfect graph knowledge hurts causal inference. Again we see the non-linear backdoor adjustment to perform similarly to DECI for ATE estimation while DECI shows decisively stronger performance when estimating CATE. As expected, the linear adjustment method fares poorly in this strongly non-linear setting.

weak_arrows: When the true causal DAG is available we find that both DoWhy methods solve this ATE problem while both DECI methods predict slightly suboptimal ATE values.

In terms of causal discovery, DECI clearly outperforms PC with the spline noise model again proving more reliable and leading to better ATE estimates. Although no methods are able to solve the task, we find that non-linear DoWhy with the DECI-spline graphs performs best. We hypothesize that the amortised function structure employed by DECI suffers in very densely connected graphs with weak edges, like is the case here.

large_backdoor: With access to the true graph, DECI methods outperform both Dowhy variants for both ATE and CATE estimation. DECI-spline performs best and is able to solve both problems. When faced with many confounders, adjustment procedures suffer from large variance. As a result, despite the non-linearity of the functional relationships at play, the simpler linear backdoor adjustment outperforms the non-linear DML approach. On the other hand, DECI’s simulation based approach is not disadvantaged in this setting.

Following the trend of the previous datasets, PC performs poorly in terms of causal discovery, biasing downstream inference methods which perform poorly in terms of ATE and CATE estimation. DECI discovery is more reliable, an effect most noticeable when using the spline noise models. With the DECI-Spline posterior over graphs, both DECI-spline and linear DoWhy are able to solve the ATE problem and DECI-spline is the only method capable of solving the CATE task. For both tasks and noise models DECI outperforms non-linear DoWhy, again showing its invariance to the size of potential adjustment set.

E.4 Synthetic Graph Experiments

We test the performance of DECI on ATE estimation with random graphs as described in section 4.1. For each graph, we randomly generate interventional data for up to five random interventions. We chose the effect variable as the last variable in the causal order that has not yet been used for data generation. For each effect variable we chose the intervention by randomly traversing the graph up to three edges away from the effect variable.

Table 3 shows the performance of the ATE estimation of DECI and all baselines on the synthetic graph data. We only show results for methods that have a runtime of less than one day. Figure 10 shows the runtimes for the different methods. DECI has consistently the lowest runtime and scales best to larger graphs. While the runtime of DECI stays approximately constant for various graphs, the runtime of the ATE estimation baselines increases with more complex graphs. In general, the methods using the true graph outperform the methods that also perform causal discovery. Further, no method strongly outperforms all other methods with DECI being a strong competitor to the already established DML methods. Lastly, we can see that DECI is capable of performing causal discovery, data imputation and ATE estimation in an end-to-end fashion without degrading performance.

E.5 Learning in Non-identifiable Settings with the Help of Graph Priors

We investigate the utility of prior knowledge over causal graphs for causal discovery and end2end inference in non-identifiable and difficult to identify settings. Specifically, we generate 2 datasets composed of 2000 training examples each. The first is composed of only linear relationships between variables and Gaussian additive noise, making the causal graph non-identifiable. The second dataset also uses linear functions but has a mix of exponential and Tanh-Gaussian noise. Although identifiable, discovery in this latter setting is challenging.

We introduce prior knowledge about graph sparseness through the weighted adjacency matrix $W_{0}\in^{D\times D}$, with zero entries encouraging sparser graphs. The resulting informed DECI prior is

with the scalar $\lambda_{s}$ regulating the strength of the prior beliefs encoded in $W_{0}$.

We compare DECI inference with access to the true graph to end2end DECI inference. In the latter case we consider a PC prior, which has as its mean the CP-DAG provided by PC. We consider different prior strengths, i.e. the value of the entries of $W_{0}$, between 0 and 1. We also experiment with introducing the true-graph as a prior of this form, yielding what we refer to as the “informed prior”.

In the non-identifiable case, we find both DECI (prior strength 0) and PC discovery to provide incorrect graphs. Interestingly, providing the PC CPDAG as a prior for DECI can yield large gains in terms of causal discovery due to a variance reduction effect. These gains do not translate to better ATE estimation, where performance is not improved over the uninformative prior. Providing knowledge of the true graph does help causal inference, with a more confident prior yielding better results.

In the difficult identifiable case, the PC prior does provide gains to DECI. We find the optimal prior strength to be 0.5: a balanced combination of PC and DECI discovery is most reliable, while using exclusively one of the two algorithms yields worse results. In this identifiable setting informed DECI discovery is able to obtain perfect ATE estimation performance with a prior strength as low as 0.2.