Assumption violations in causal discovery and the robustness of score matching
Authors: Francesco Montagna
Atalanti A. Mastakouri
Elias Eulig
Nicoletta Noceti
Lorenzo Rosasco
Dominik Janzing
Bryon Aragam
Francesco Locatello
Venue: NeurIPS 2023 poster (2023)
Keywords: Causal discovery; empirical study; robust inference; benchmark
TL;DR: The challenge of causal structure learning lies in the modeling assumptions of the data, which are often untestable. We investigate the performance of recent causal discovery methods when assumptions on the data may be violated.

Abstract: When domain knowledge is limited and experimentation is restricted by ethical, financial, or time constraints, practitioners turn to observational causal discovery methods to recover the causal structure, exploiting the statistical properties of their data. Because causal discovery without further assumptions is an ill-posed problem, each algorithm comes with its own set of usually untestable assumptions, some of which are hard to meet in real datasets. Motivated by these considerations, this paper extensively benchmarks the empirical performance of recent causal discovery methods on observational _iid_ data generated under different background conditions, allowing for violations of the critical assumptions required by each selected approach. 
Our experimental findings show that score matching-based methods demonstrate surprising performance in the false positive and false negative rate of the inferred graph in these challenging scenarios, and we provide theoretical insights into their performance. This work is also the first effort to benchmark the stability of causal discovery algorithms with respect to the values of their hyperparameters. Finally, we hope this paper will set a new standard for the evaluation of causal discovery methods and can serve as an accessible entry point for practitioners interested in the field, highlighting the empirical implications of different algorithm choices.

Assumption violations in causal discovery and the robustness of score matching
Francesco Montagna MaL.a, Università di Genova
Atalanti Mastakouri AWS
Elias Fulig
Nicoletta Noceti
German Cancer Research Center (DKFZ)
MaLGa, Università di Genova
Lorenzo Rosasco
Dominik Janzing AWS
Bryon Aragam University of Chicago
MaLGa, Università di Genova MIT, CBMM Istituto Italiano di Tecnologia
Francesco Locatello
Institute of Science and Technology Austria (ISTA)
Abstract
When domain knowledge is limited and experimentation is restricted by ethical, financial, or time constraints, practitioners turn to observational causal discovery methods to recover the causal structure, exploiting the statistical properties of their data. Because causal discovery without further assumptions an ill-posed problem, each algorithm comes with its own set of usually untestable assumptions, some of which are hard meet in real datasets. Motivated by these considerations, this paper extensively benchmarks the empirical performance of recent causal discovery methods on observational iid data generated under different background conditions, allowing for violations of the critical assumptions required by each selected ap- proach. Our experimental findings show that score matching- saasching-based methods demon- strate surprising performance in the false positive and false negative rate of the inferred graph in these challenging scenarios, and we provide theoretical insights into their performance. This work is also the first effort to benchmark the stability of causal discovery algorithms with respect to the values of their hyperparameters. Fi- nally, hope this paper will ery methods and can serve as an accessible entry point for practitioners interested in the field, highlighting the empirical implications of different algorithm choices.
standard for the evaluation of causal discov-
1 Introduction
The ability to infer causal relationships from observational data, instead of simple statistical asso- ciations, is crucial to answer interventional and counterfactual queries without direct manipulation of a system [T] 2 3 4 The challenge of drawing causal conclusions from pure observations lies the data, which are often impossible verify. Methods based
in the modeling assumptions
conditional independence testing (e.g. PC, FCI and their variations 4 5 ( require faithfulness of
the distribution 1 2 7 to the causal graph, which formalizes the intuition that causal relations
manifest themselves in the form of statistical dependencies among the variables. The assumption of causal sufficiency (i.e. the absence of unobserved confounders (8) is common requirement for
37th Conference on Neural Information Processing Systems (NcurIPS 2023).

---

12,, which allows interpreting associations in the data as causal relationships. These strong conditions are arguably necessary but nevertheless hard or impossible to verify, and posit an entry barrier to the unobscured application of causal analysis in general settings. In addition to that, structure identifiability results define limitations on the parts of the causal graph that can be inferred from pure observations [ TO [3]. Traditional causal discovery methods (e.g. PC, FCI. GES [4] D
causal discovery 4 TO
are limited to the inference of the Markov Equivalence Class of the ground truth graph 14, while additional assumptions the structural equations generating effects from the cause ensure identifiability of a unique Directed Acyclic Graph (DAG) from observational data. In partic- ular, restrictions on the class of functions generating the data (linear or not) and on the distribution of the noise terms (i.e. additive noise models assumed in LINGAM and more) characterizing their non-deterministic relationships are necessary in order to infer causal directions [T]]
[3] Re-
quirements on the data collection process are also needed: although an error-free measurement model is commonly assumed, it has been recent subject of interest that measurement error in the observed values of the variables can greatly change the output of various causal discovery methods [15 16].
Real data hardly satisfy all of these assumptions at once, and it is often the case that these are impos- sible to verify, which calls for algorithms that demonstrate a certain degree of robustness with respect to violations of the model hypothesis. Previous work from Heinze-Deml et al. [[7 investigates the boundaries of robust graph inference under model misspecifications, on Structural Causal Models (SCM) with linear functional mechanisms. Mooij et al. [18] benchmark considers the case of additive noise models with nonlinear mechanisms, but only for datasets with two variables. Singh et al. [99] presents empirical evaluation limited to methods whose output is Markov Equivalence Class. Glymour et al. [4] review some of the existing approaches with particular attention to their required assumptions, but without experimental support to their analysis. Our paper presents an extensive pirical study that evaluates the performance of classical and recent causal discovery methods on obser- vational datasets generated from iid distributions under diverse background conditions. Notably, the effects of these conditions on most of the methods included in our benchmark have not been previously investigated. We compare causal discovery algorithms from the constraint and score-based literature, as well as methods based on restricted functional causal models of the family of additive nonlinear 13 20]. These include recent class of methods deriving connections between the score models [T] [T
matching 21 22 with the structure of the causal graph [3 24 255. Algorithms that focus on sequen-
tial data (26 27 g
33 36 37 are beyond the scope of this paper's benchmark- ing. Finally, we propose an experimental analysis of the stability of the benchmarked approaches with respect to the choices of their hyperparameters, which is the first effort of this type in the literature.
29 30 32
We summarise the contributions of our paper as follows:
We investigate the performance of current causal discovery methods in large-scale experimental study datasets generated under different background conditions with violations of the required background assumptions. Our experimental protocol consists of more than 2M experiments with 11 different causal discovery methods more than 60000 datasets synthetically generated.
We release the code for the generation of the synthetic data and a Python implementation of six causal discovery algorithms with a shared API. With this contribution, we aim at facilitating the benchmarking of future work in causal discovery on challenging scenarios, and the comparison with the most prominent existing baselines.
We analyze our experimental results, and present theoretical insights on why score aatching-based approaches show better robustness in the setting where assumptions on the data may be violated, compared to the other methods. Based on our empirical evidence, we suggest a new research direction focused on understanding the role of the statistical estimation algorithms applied for causal inference, and the connection of their inductive biases with good empirical performance.
2 The causal model
In this section, we define the problem of causal discovery, with brief introduction to the formalism of Structural Causal Models (SCMs). Then we provide an overview of SCMs for which sufficient conditions for the identifiability of the causal graph from observational data are known.

---

2.1 Problem definition
A Structural Causal Model M is defined by the set of endogenous variables X E R vertices of the causal graph 9 that we want identify, the exogenous noise terms U E R distributed according to pu. as well as the functional mechanisms F (fi...., fa), assigning the value of the variables X1,..., Xa as a deterministic function of their causes and of some random disturbance.
Each variable X; is defined by a structural equation:
X; = f(A,,, U.), =i= d.
where PA, X is the set of parents of X, in the causal graph G. and denotes the set of direct causes of X₁. Under this model, the recursive application of Equation (T entails joint distribution px. such that the Markov factorization holds:
d
px (X)
II p. (X,))).))
The goal of causal discovery is to infer the causal graph underlying X from set of observations sampled from px-
2.2 Identifiable models
In order to identify the causal graph of X Rd from purely observational data, further assumptions on the functional mechanisms in F and on the joint distribution pu of model ( are needed. Intuitively, having
condition between nonlinearity of the causal mechanisms and non- -ou-siassiyyyty of the noise terms is necessary to ensure the identifiability of the causal structure. Additionally, we consider causal sufficiency (Appendix A.3) of the model to be satisfied, unless differently specified.
Linear Non-Gaus Model (LINGAM). A linear SCM is defined by the system of structural equations
X = 1
B E RRd is the matrix of the coefficients that define X; as a linear combination its parents and the disturbance U,. Under the assumption of Gaussian distribution of the noise terms, the model is identifiable. This SCM is known as the LiNGAM (Linear Non-Gaussian Acyclic Model) [2]
Additive Noise Model. An Additive Noise Model (ANM) [ is defined by Equation ① when it represents the causal effects with nonlinear functional mechanisms and additive noise terms:
X, = (PA.) + U., Vi I,...,d,
with fo nonlinear. Additional conditions on the class FF of functional mechanisms and on the joint distribution of the noise terms are needed to ensure identifiability [13. In the remainder of the paper, we assume these to hold when referring ANMs.
Post Non Linear Model. The most general model for which sufficient conditions for the identifiability of the graph are known is the Post NonLinear model (PNL) [10. In this setting the structural equation (I) can be written as:
x; = (PAA)) + U.), 1......,,
(5)
where both Ii and fi are nonlinear functions and g; invertible. As for ANMs, we consider identifiability conditions defined in Zhang and Hyvärinen [10 to be satisfied in the rest of the paper.
3
Experimental design
In this section, we describe the experimental design choices regarding the generation of the synthetic datasets, the evaluated methods, and the selected metrics.
3.1 Datasets
The challenge of causal structure learning lies in the modeling assumptions of the data, which are often untestable. Our aim to investigate the performance of existing causal discovery methods in the setting where these assumptions are violated. To this end, we generate synthetic datasets under

---

DireaLiNGAA
GraN-DAG
PC
PCT
GES
RESYT
CAM
SCORE
DAS
NoGAM
DIIANN
Gaussian noise
Non-Gaussian noise"
Lincar mechanisms
Nonlinear mechanisms
Unfaithful distribution
X
Confounding effects
X
X
X
X
Measure
X
X
X
X
X
X
X
Output
CPDAG
PAG
CPDAG
DAG
DAG
DAG
DAG
DAG
DAG
DAG
DAG
GraN-DAG and GES optimize the Gaussian likelihood.
Table I: Summary of the methods assumptions and their output graph. The content of the cells
denotes whether the method supports (V) or not (X) the condition specified the corresponding row.
diverse background conditions, defined by modeling assumptions that do not match the working hypothesis of the evaluated methods.
Vanilla model. First, we specify an additive noise model with variables generated according to the structural equation
The exogenous terms follow Gaussian distribution U; (0, o;) with variance o; U (0.5, 1.0) uniformly sampled. We generate the nonlinear mechanisms fi by sampling Gaussian processes with a unit bandwidth RBF kernel (Appendix [.T]. We refer to this model as the vanilla scenario, as it is at one time both identifiable and compliant with the assumptions of the majority the benchmarked methods (see Table !
3.1.1 Misspecified scenarios
We define additional scenarios such that each specified model targets a specific assumption violation with respect to the vanilla conditions.
Confounded model. Let Z E Edd be a of latent common causes. For each pair of distinct nodes X; and Xj. we sample a Bernoulli random variable Cgj ~ Bernoulli(p) such that Ci = 1 implies a confounding effect between X; and Xj- The index & of the confounder Zμ is assigned at random. The parameter Measurement error model. Measurement errors in the process that generates the data regarded as source of mistakes for the inference of the causal graph [13 [6]. In order to account for potential errors induced by the measurements, specify a model in which the observed variables are:
{0.1, 0.2) determines the amount of confounded pairs in the graph.
8 = X4 ej Vi 11=11111.....: I,...,d,
(6)
a noisy version of the X;'s generated by the ANM of Equation (4). The e; disturbances are indepen- dent Gaussian random variables centered zero, whose variance parametrized by the inverse signal- to-noise ratio % = (ar(X;) Given that the total variance of X; is Var(X,) = Var(X,) + Var(e;), % controls the amount of variance in the observations that is explained by the error in the measurement. Each dataset with measurement error parametrized with > {0.2, 0.4, 0.6, 0.8). shared by all the ej.
Varie;)
Unfaithful model. To model violations of the faithfulness assumption (Appendix A.2). we tune the causal mechanisms of Equation 4 such that we induce direct cancellation of causal effects between some variables. In particular, for each triplet X; - Xx + X, +- X; in the graph, causal mechanisms are adjusted such that cancellation of the causal effect X, - X₁ occurs (for implementation details, see Appendix [B.4). This is a partial model of unfaithfulness, as only covers a limited subset of the scenarios under which unfaithful path canceling might occur, and must be viewed in the light that there is no established procedure to enforee unfaithful conditional independencies in the case of ANM with nonlinear relationships.
Autoregressive model. In order to simulate violations of the iid distribution of the data, we model observations as a stochastic process where each sample is indexed by time. In particular, we define the structural equations generating the data
(t) ax,(t f((A,(tt) + U. t= 1,2,3 a E.
(7)

---

Autoregressive effects are modeled with a time lag 1 = 1, whereas at 0 we define (0) with Equation ( The ground truth is the graph whose edges represent the connections between parents (0) and their direct effect (0).
Post Jonnr model. We replace nonlinear causal mechanisms of the additive noise models (4) with the structural equations defined in the PNI. model ( We select the post nonlinear function 9i such that (z()) = x², E IR, Vi ,,...,d.
LINGAM model. We define model with the linear system of structural equations 3 The non- non-Gaussian distribution of the noise terms is defined as a nonlinear transformation of a standard
normal random variable (see Appendix B.2). and the linear mechanisms are simulated by sampling the weighting coefficients of the parents of a node in the interval |-1,-0.00 -0.05] u [0.05, 1].
3.1.2 Data generation
For each specified model, we generate datasets that differ under the following characteristics: num- ber of nodes d E {5, 10, 20, 50), number of samples n E {100, 1000} and density of edges p E {sparse, dense}. We sample the ground truth causal structures according to different algorithms for random graph generation. In line with previous causal discovery literature
we gen- erate Erdos-Renyi (ER) 40) and Scale-free (SF) graphs [4]. Furthermore, we consider Gaussian Ran- dom Partitions (GRP)
25
38 39
42 and Fully Connected graphs (FC) (see Appendix BB.3). By considering all the combinations of the number of nodes, number of samples, admitted edge densities, and algorithms for structure generation, we define a cartesian product with all the graph configurations of interest. For each of such configurations and for each modeling scenario, we generate dataset 2> and its ground truth G with 20 different random seeds. Details on the generated data can be found in Appendix B3.
3.2 Methods
We consider 11 different algorithms and a random baseline spanning across the main families of causal discovery approaches: constraint and score-based methods, and methods defined under restrictions the structural causal equations. In the main text, provide detailed overview of the methods most relevant for the discussion of our key experimental findings. The remaining approaches are described in further detail in the Appendix C Table summarizes the algorithms' assumptions and the output object of their inference procedure.
Method outputs. Causal discovery algorithms output different graphical objects based on their underlying assumptions. If identifiability is not implied by the model requirements but faithfulness of the distribution is satisfied, one can instead recover the Markov equivalence class of the ground truth graph, that is, the of DAGs sharing the same conditional independencies. This is represented by a complete partially directed acyclic graph (CPDAG), where undirected edges X; - X, are meant encode conditional dependence between the variables, but uncertainty the edge orientation. If a method can identify a directed acyclic graph G = (X. E), one can define a partial ordering of the nodes π {X, ) J. T E {I...., d). such that whenever we have Xx, + X1, E. then X=, < <X,, (X n, is a successor of X_, in the ordering) [43] The permutation π is known as the topological order of G, and allows to disambiguate the direction of the edges in the graph. A topological order can be encoded in fully connected DAG with edges & = [Xx, T. + Xx, Xx,, x
Xz, Vi,j = I,....d], d). obtained connecting all nodes in the ordering ग with their successors.
Methods summary. A summary of all the methods included in the benchmark and their required p are limited identifying the Markov equivalence class of the DAG. Direct LiNGAM 44 is designed for inference on data generated by a linear Gan-Gaussian model whereas SCORE 23 NoGAM (25, DiffAN 45). DAS (24, RESIT (13). GraN-DAG (38) and CAM [46), are meant for inference on additive noise models: these methods perform inference in two steps procedure, first identifying topological ordering of the graph, and then selecting edges between those admitted by the inferred causal order. To enable fair comparison in our experiments, all methods (with the exception of DirectLiNGAM) are implemented with the same algorithm for edge detection, consisting of variable selection with sparse regression. This pruning strategy is known as CAM-pruning, being originally proposed in CAM paper (46). A detailed discus- sion of all the methods in the benchmark is presented in Appendix C In the Appendix [ we consider experiments on FCI which are not reported in the main text since we did not find metrics for a straightforward comparison of its output partial ancestral graph (PAG [47 with CPDAGs and DAGs.
assumptions is presented in Table T PC 4 and GES
5

---

Selected metrics To evaluate the output graphs we use the false positive and false negative rates, and the Fl score (details in the Appendix (D). the case of directed edges inferred with reversed direction, we count this error false negative. For methods that output a CPDAG with undirected edges, we evaluate them favorably by assuming correct orientation whenever possible, similar Zheng et al. [39] 48]. For the methods whose output also includes an estimate á of the topological order, we define the false negative rate of the fully connected DAG with edges &g {Xx, + Xa, Xs; - Xx, Vi,j = 1..... d),, denoted as FNR-î. If स is correct with respect to the ground truth graph, then FNR-#
This choice of metries reflects the implementation of most of the algorithms involved in the benchmark, which separates the topological ordering step from the actual edge selection. In particular, given that the majority of the methods share the same pruning procedure after the inference of the order, we expect that differences in the performance will be mostly observed in the FNR-# score.
3.2.1 Deepdive on SCORE, NoGAM and DiffAN
In this section, we review recent class of causal discovery algorithms, that derive constraints on the score function V log)(XX that uniquely identifies the directed causal graph of an additive noise model. Identifiability assumptions provide sufficient conditions to map a joint distribution px to the unique causal DAG 9 induced by the underlying SCM. Applying the logarithm to the Markov factorization of
the distribution of Equation
we observe that log px (X) ΣS log p(X, PA.). By inspection of [). it is possible to derive constraints mapping the score function to the
the gradient vector V log px
causal graph of ANM. Given node X, in the graph, its corresponding score entry is defined as:
8, X == 0x, Px XX = 0x, log (X, PA.)) +
Σ
0x, Pa PAa)-
(8)
KECH.
Instead, the rate of change of the log- log-likeldhood over a leaf node X₁ with the set of children CH2 = is:
ss (X) = 0x, log (x() = 0x, log (o.
99)
We see that, for a leaf node, the summation over the set of children of Equation 8 vanishes. Intuitively, being able to capture this asymmetry in the score entries enables the identification of the topological order of the causal graph.
SCORE. The SCORE algorithm (33 identifies the topological order of ANMs with Gaussian noise terms by iteratively finding leaf nodes as the argmin, Var [xx, Sx (X)], given that the following holds:
Var [0x,s; (X)] 0 + X, is a leaf, Vi = 1,...,d. d.
(10)
NoGAM. The NoGAM algorithm generalizes the ideas of SCORE on additive noise models with an arbitrary distribution of the noise terms. After some manipulations, it can be shown that for satisfies
leaf node X₁ the score entry of Equation 9
8t (X) dc, (u())
(11)
such that one could learn consistent estimator of st taking input the exogenous variable Uj. For an ANM, the authors of NoGAM show that the noise term of a leaf is equivalent to the residual defined as:
R4 = X; E
[
X \
{X,)], Vi = I,...,d.
(12)
Then, by replacing U2 with R2 in Equation (T). it is possible to find a consistent approximator of the score of a leaf using R₁ the predictor. Formally:
E
E (X)
R. - (x))
- -
X; is a leaf,
(13)
which identifies leaf node the argmin of the vector of the mean squared errors of the regression of the score entries 8; (X) on the corresponding residuals R;, for all i 1,...,d. d.
Connection of NoGAM with the post nonlinear model. is interesting to notice that, similarly to Equation (1T) for additive noise models, the score of leaf X, generated by a PNL. model can be defined as a function of the disturbance Ut.
Proposition 1. Let X Rdd be generated according to the post nonlinear model Then, the score function of a leaf node X₁ satisfies 81 (xx) = Oi log P(()).
This result suggests a connection with the NoGAM sorting criterion: indeed, one could hope iden- tify leaf nodes the graph by consistent estimation of the score a leaf from residuals equivalent to the noise terms. A more detailed discussion with the proof of Proposition iis presented in Appendix E

---

DAS. The DAS algorithm (acronym for Discovery At Scale, (24) identifies the topologiaa ordering with the same procedure defined in SCORE, while the two methods differ in the way they find edges in the graph. DAS edge selection procedure exploits the information in the non-diagonal entries of the Jacohian of the score. In particular, for ANM with Gaussian noise terms, it can be shown that:
(14) Exploiting Equation t4 define the inference rule for the edge selection in DAS provides a significant computational advantage with respect to SCORE, reducing the time complexity in the number of
E [ st [X]]] * 0 - X3 E (X), Vj E (1....,d) d} {
nodes from cubic to quadratic.
DiffAN. DiffAN 45 method finds the topological ordering of a DAG exploiting the same criterion of Equation (TO) of SCORE: the difference is in that it estimates the score function with probabilistic dif- fusion models, whereas SCORE, NoGAM, and DAS 24 rely on score matching estimation 21 22
4 Key experimental results and analysis
In this section we present our experimental findings on datasets generated according to the misspecified models of Section 3.I.IJ with theoretical insights into the performance of score matching- based approaches. We draw our conclusions by comparing the methods' performance I
against their accuracy in the vanilla scenario and against a random baseline [defined in Appendix
C.10). The results are discussed on datasets of size 1000 for Erdos-Renyí dense graphs with 20 nodes (ER-20 dense), and can be generalized different size and sparsity configurations, Due to space constraints, we include the plots only for the Fl score and FNR-f. whereas the false negative and false positive rates are discussed in Appendix []] In order to provide statistical significance to our conclusions, we repeat the experiments on each scenario over 20 datasets generated with different random seeds. To enable a fair comparison between the methods, we fix their hyperparameters to their optimal value with respect to each specific dataset, in the case where these can not be tuned without having access to the ground truth (see Appendix Oppendix[Gffor for a discussion on the tuning of GraNDAG and DiffAN learning hyperparameters) In the Appendix H the benchmarked methods with respect to different values of their hyperparameters.
we analyze the stability of
4.1 Can current methods infer causality when assumptions on the data are violated?
Our experimental findings suggest that score matching- baseed algorithms can robustly infer part of the causal information even in the case of misspecified ground truth data generation.
Post nonlinear model. Figure la (right) illustrates the accuracy of topological order estimates on post nonlinear model data. Among the selected methods, NoGAM shows better ability to generalize its performance to this scenario, with FNR- eгror rate significantly lower than the random bascline. Interestingly, we can interpret this observation in the light of Proposition [ which defines the score of leaf in the PNL model: our result indeed suggests that, similarly to the case of an additive noise model, it is possible to learn a consistent approximator of the score of leaf X2 from the exogenous variable Ui of a post nonlinear model. Notably, we also observe that RESIT order accuracy is better in the PNL scenario than in the vanilla case: Zhang and Hyvärinen [00 show that testing for independent residuals identifies the direction of causal relationships also under the PNL model.
LINGAM model. Figure Ib (right) shows that NoGAM can infer the causal order with remarkable accuracy in the case of ground truth data generated by a linear non-gaussian additive model. Together with our observations on the post nonlinear model, our empirical evidence corroborates the idea that the NoGAM algorithm is surprisingly robust with respect to the misspecification of the causal mechanisms. Notably, none of the other methods can infer the ordering with accuracy significantly better than the random baseline. This could lead to decreased performance in the realistic setting of mixed linear and nonlinear mechanisms. However, the Fll score in Figure [bb (left) shows that CAM-pruning is still able to correctly infer edges in the graph when these are admitted by the identified causal order. We note that, given that we observed high varsortability [50] for this model, display results on data standardized dividing by their empirical variance.
'We use the https: / gitb. com/ ct/c//////nnnaa implementations of RESIT and DirectLiNGAM, and the DoDicoover implementations of PC, GES, and FCI. For the remaining methods, we consider the GitHub official repositories of their papers and custom implementations.
2 Varsortabiliity of a dataset denotes partial agreement between the ordering induced by the values of marginal variance of the observed variables and the causal ordering of the underlying graphical model. Note that other simple baselines can still recover the causal order after rescaling of the data 44).
7

---

anden
das
Cam
Score
nogam
diffr
grandag
rest
pC
ges
misspecrfed
PNL model
BR-20 dense
LINGAM ER 20 dense
F1 score
FNR--
Flocore
FNR-T
1.0
0.
9.6
0.G
0.4
TA
0.4
0:2
0.0
(a) PNL
(b) LiNGAM
Confounded model- ER-20 dense
Measurement model ER-20 dense
0.0
0.0
0.2
22
0.2
0.2
1
(c) Latent confounders
(d) Measurement error
Unfaithful model ER-20 dense
Autoregressive model ER-20 dense
3.5
score
1.
PNR-6
1.0
score
1.0
FNBR-A
0.
0.6
6.2
2.0
0.0
(c) Unfaithful distribution
(f) Non-iid data distribution
Figure 1: Experimental results on the misspecified scenarios. For each method, we also display the violin plot of its performance on the vanilla scenario with transparent color. Fl score (the higher the better) and FNR-A - (the lower the better) are evaluated over 20 seeds on Erdos- Renyi dense graphs with 20 nodes (ER 20 dense). FNR-7 ft is not computed for GES and PC, methods whose output CPDAG. Note that DirectLiNGAM performance reported in Appendix [.2] on data under non-Gaussian distribution of the noise terms.
Confounded model. Spurious correlations, that occur when the causal sufficiency is violated, can not be handled by statistical tests for edge selection, as shown by the Fl score of Figure IC (left) (the amount of confounders is parametrized by 0.2). In this case, we are also interested to see whether the presence of latent confounders can disrupt the inference of the topological ordering when the observed variables have non-spurious connection in the causal graph. Figure Ic (right) indicates that the score matching- based approaches SCORE, DAS, and NoGAM can still be exploited to find a reliable ordering, while other methods fail to do so.
Measurement error. Given data generated under the model of Equation (6). we observe convergence
in distribution p(X; PAA)) S p(X; PA.)) for + 0. We are then interested in the boundaries of robust performance of each method with respect to increasing values of 7. Figure (Id) (right) illustrates FNR-A on datasets with 7 0.8 such that ~35% of the observed variance of each variable is due to noise in the measurements. Under these conditions, we see that score matching- approaches display robustness in the inference of the order where all the other methods' capability is comparable to that of the random bascline with statistical significance. This is also reflected in Figure (Td) (left), where SCORE, DAS, and NoGAM are the only algorithms whose Fl score (slightly) improves compared to the random baseline.
Unfaithful model. Figure Tel (right) shows that the ordering procedure of several methods, in particular SCORE, DAS, NoGAM, and GraN-DAG, seems unaffected by direct cancellation of causal effects, in fact displaying surprising decrease in the FNR-* performance with respect to the vanilla scenario. To understand these results, we note that under the occurrence of causal effect cancellations in the ground truth graph G, the unfaithful model defined in Section 3.1.I generates observations of X according to graph & whose causal order agrees with that of the ground truth: it
8

---

is indeed immediate to see that the causal order of X, + X X, also holds for the triplet X; + X, + Xg. Moreover, the set of edges of the graph ف is sparser than that of the ground truth, due to the cancellation of causal effects. Thus, given that inference on sparser graphs is generally easier, it can positively affect the empirical performance, in line with our observations.
Implications. Our experimental findings show that most of the benchmarked methods significantly decrease their performance the misspecified models. This is particularly problematic since the violations considered in this work are realistic and met on many real-world data. On the other hand, we observe surprising robustness in the inference of score matching-baser methods.
4.1.1 Discussion on PC and GES performance
The experimental results of Figure I show that the Fl score of GES and PC is consistently worse than random in the setting of Erdös-Reny' dense graphs with 20 nodes. We note that this pertains specifically to large graphs with dense connections, where the accuracy of both methods significantly degrades. In the case of PC, this is in line with previous theoretical findings in the literature: Uhler et al. 77 demonstrate that large and dense graphs are characterized by many quasivvioaations of the faithfulness assumptions, which can negatively impact the algorithm's inference ability. In Appendix M we discuss experiments on lower dimensional networks with sparse edge structures, showing that PC and GES report performance significantly better than random across different scenarios. Overall, we find that increasing the size and density of the graph negatively impacts the inference ability of PC and GES.
4.1.2 Discussion on score matching robustness
Our empirical findings indicate that score matching-basee methods are surprisingly capable of partial recovery of the graph structure in several of the misspecified scenarios. We connect this robust performance to the decomposition properties of the score function defined in Equations (8) and ( In particular, we argue that the common factor that enables leaf node identification in NoGAM and SCORE is that the score entry of leaf is characterized by smaller magnitude, compared to the score associated with a node that has children in the graph. To explain what we mean by this, we define simple condition under which it is possible to identify Icaf nodes and the causal order of the graph from the variance of the entries of the score function.
Definition Let X E IRd be a random vector defined by a structural causal model M T Let X be a leaf node of the causal graph g. We say that X₁ is score-ideniifable if I = argmin, Var(s; (X)).
Moreover, we say that the model is score-sortable if the recursive identification of ideetiidete leaf nodes in the causal graph and in the subgraphs defined by removing Icaf from the set of vertices up to a source node, yields a correct causal order. SCORE, NoGAM, and DAS present results for consistent inference of the structure of identifiable graph from properties of the score function and its second order partial derivatives. However, when these conditions are not satisfied, exploitation of score-sortability can heuristically estimate a causal ordering that partially agrees with the causal structure. Intuitively, the variance of the score of a non-leaf node 8; (X) of Equation (8) is proportional to the number of children in the summation. In particular, the total variance of 8; is the sum of the marginal variances of the two terms on the RHS of Equation (8), plus their covariance.
Errors in the ordering defined with score- sortabiiyyyyyyty are induced only if the variance associated with the score of a non-leat node can be smaller than the one relative to every leaf of the graph.
Proposition 2. Let X E TRd be a random vector whose elements X; are defined by a structural equation model M (1) that satisfies corte-sortabilty Then, for each subgraph of G defined by recursively removing a from the set of vertices up to source node, there exists a leaf X, such that Vi index of a node:
Var log (x PAA)) & Var(ūx, Ix, + Σ
Varlax, log PAL)] + C,
AECH.
with C E R accounting for the covariance term.
(See Appendix F for the proof.) Lemma 1 of SCORE defines a similar criterion of sortability of the causal variables on the variance of the second order partial derivatives of the log-likelihood which is always satisfied when X E gd is generated by an ANM with Gaussian distribution of

---

the noise terms, We can extend these considerations to the NoGAM algorithm, which identifies leaf nodes by minimizing the mean squared error of the predictions of the score entries from the residual estimators of the noise terms, as defined in Equation (T3). If we consider an uninformative predictor of the score function that maps every input residual to a constant value zero, the NoGAM algorithm is equivalent to a simple score-sortability heuristic criterion, identifying leaf nodes as the argmin; {|s ( In Appendix 1.3 we corroborate our considerations by comparing the empirical performance of score-sortability baseline with SCORE and NoGAM.
Implications. Score matching- based approaches SCORE, DAS, and NoGAM show empirical robust- ness in several scenarios included in our benchmark. We impute these results to the structure of the score function discussed in Section (3.2.1). and the algorithmic design choices of these methods that exploit different magnitude in the score of a leaf compared to other nodes with children in the graph.
4.2 Is the choice of statistical estimators neutral?
In the previous section, we motivated the empirical observations on the robustness of methods based on the score function. Given that the DiffAN algorithm differs from SCORE only in the score estimation procedure (where the former applies probabilistic diffusion models in place of the score matching), we can explain the gap in performance of DiffAN with the other approaches based on the score an effect of the different statistical estimation technique. From this observation, we suggest that score matching plays crucial role in connecting the gradient of the log-likelihood with effective causal inference.
Implications. The choice of modular statistical estimator for causal inference procedures is not neutral. We argue that inductive bias in statistical estimators may be connected with good empirical performance, and we think that this potential connection should be further investigated in future works.
5 Conclusion
In this work we perform large-scale empirical study on eleven causal discovery methods that provide empirical evidence on the limits reliable causal inference when the available data violate critical algorithmic assumptions. Our experimental findings highlight that score matching- approaches can robustly infer the causal order from data generated by misspecified models. It would be important have procedures for edge detection that display the same properties of robustness in diverse scenarios and to have better theoretical understanding of failure modes of CAM-pruming variable selection, given its broad use for causal discovery. Finally, we remark that this benchmarking is limited to the case of observational iid samples, and it would be of great practical interest to have equivalent empirical insights on the robustness of methods for causal discovery on sequential data in the setting of time series passively observed interventions.
Acknowledgements
We thank Kun Zhang and Carl-Johann Simon-Gabriel for the insightful discussions. This work has been supported by AFOSR, grant n. FA8655-20-1-7035. FM is supported by Programma Operativo Nazionale ricerca innovazione 2014-2020. FM partially contributed to this work during an internship at Amazon Web Services with FL. FL partially contributed while AWS.
References
[1] Jonas Peters, Dominik Janzing, and Bernhard Schlkopf. Elements of Causal Inference: Founda- lions and Learning Algorithms. The MIT Press, 2017. ISBN 0262037319.
[2] Judea Pearl. Causality. Cambridge University Press, Cambridge, UK, 2 edition, 2009. ISBN 778-0-521-895601 doi: 10. 10.1001//8097805111 180078000111011666 1803161.
[3] Peter Spirtes. Introduction to causal inference. Journal of Machine Learning Research, 11(54): 1643-1662, 2010. URL http::/jmmm http://aaers/////pppppssssssseee
html
[4] Spirtes, C. Glymour, and Scheines. Causation, Prediction, and Search. MIT press, 2nd edition, 2000.
10

---

[5] Naftali Harris and Mathias Drton. Pc algorithm for nonparanormal graphical models. Journal of Machine Learning Research, 14(69):3365-3333 2013. URL http: httpj//l...oo/ppp/eeeeeeeeeeeeee org/papers/
114/harris13) html
[6] Joseph Ramsey, Peter Spirtes, and Jiji Zhang. Adjacency-Faithfuunes and Conservative Causal Inference. In Proceedings of the 22nd Conference on Uncertainty in Artificial Intelligence, pages 401-408, 2006.
[7] Caroline Uhler, G. Raskutti, Peter Bühlmann, and B. Yu. Geometry of the faithfulness assump- tion in causal inference. The Annals of Statistics, 41, 07 2012. doi: 10.1214/12-AOS
[8] Hans Reichenbach. The Direction of Time. Mineola, N.Y.: Dover Publications, 1956.
[9] David Maxwell Chickering. Optimal structure identification with greedy search. J. Mach. Learn. Rex., 3005))))44) mar 2003. ISSN 1532-4435. doi: 10.112/153224303211877777
URL
2stps://ddi..org/100
1162/15524403332187771
[10] Kun Zhang and Aapo Hyvärinen. the identifiability of the post-nonlinear causal model. In Proceedings of the Twenty-Fift Conference Uncertainty in Artificial Intelligence. UAI '09, page 647-655, Arlington, Virginia, USA, 2009. AUAI Press. ISBN 9780974903958.
[11] Patrik Hoyer, Dominik Janzing, Joris M Mooij, Jonas Peters, and Bernhard Schölkopf. Non- linear causal discovery with additive noise models. In D. Koller, D. Schuurmans, Y. Bengio, and L.. Bottou, editors, Advances in Neural Information Processing Systems, volume 21. Cur- ran Associates, Inc., 2008. htt.httts://pproeeeeee
pdf
[12] Shohei Shimizu, Patrik O. Hoyer, Aapo Hyvärinen, and Antti Kerminen. A linear non-gaussian acyclic model for causal discovery. J. Mach. Learn. Res., 7:2003-2030, dec 2006. ISSN 1532-4435.
[13] Jonas Peters, Joris M. Mooij, Dominik Janzing, and Bernhard Schölkopf. Causal discovery with continuous additive noise models. Mach. Learn. Res., ([11):2009-05533 jan 2014. ISSN 1532-4435.
[14] Clark Glymour, Kun Zhang, and Peter Spirtes. Review of causal discovery methods based on graphical models. Frontiers in Genetics, 10, 2019. ISSN 1664-8021. doi: 10.3389/fgene. 2019. 00524. httt:s:://www.froiiirrssiiii org/ prg/articles//1 3389/fgene 2019.00524
[15] Kun Zhang, Mingming Gong, Joseph Ramsey, Kayhan Batmanghelich, Peter Spirtes, and Clark Glymour. Causal discovery in the presence of measurement error: Identifiability conditions. 06 2017.
[16] Richard Scheines and Joseph Ramsey. Measurement error and causal discovery. CEUR workshop proceedings, 1792:1-7, 06 2016.
[17] Christina Heinze-Deml, Marloes H. Maathuis, and Nicolai Meinshausen. Causal structure learning. Annual Review of Statistics and Its Application, 5(1):371-391, 2018. doi: 10.1 14.1 minuueesstaiistie--0111---000000
[18] Joris Mooij, Jonas Peters, Dominik Janzing, Jakob Zscheischler, and Bernhard Schölkopf. Distinguishing cause from effect using observational data: Methods and benchmarks. Journal of Machine Learning Research, 17:1-102, 04 2016.
[19] Karamjit Singh, Garima Gupta, Vartika Tewari, and Gautam Shroff. Comparative benchmarking of causal discovery techniques. 08 2017.
[20] Kun Zhang and Aapo Hyvärinen. Causality discovery with additive disturbances: An informaaiion-theorttical perspective. In Wray Buntine, Marko Grobelnik, Dunja Mladenić, and John Shawe-Taylon Taylor, editors, Machine Learning and Knowledge Discovery in Databases, pages 570-585, Berlin, Hcidelberg, 2009. Springer Berlin Heidelberg.
[21] Aapo Hyvärinen. Estimation of non-norma statistical models by score matching. Journal of Machine Learning Research. 6(24):695 2005. htLh:tpp: http:://mmmmr org wwwgg/paprs/vvee /v6/
hyvarinen05a. html
[22) Yingzhen Li and Richard Turner Gradient estimators for implicit models. 05 2017.
[23] Paul Rolland, Volkan Cevher, Mathaus Kleindessner, Chris Russell, Dominik Janzing, Bernhard Schölkopf, and Francesco Locatello. Score matching enables causal discovery of nonlinear additive noise models. In Kamalika Chaudhuri, Stefanie Jegelka, Le Song, Csaba Szepesvari,
11

---

Gang Niu, and Sivan Sabato, editors, Proceedings of the 39th International Conference on Machine Learning, volume 162 of Proceedings of Machine Learning Research, pages 18741- 18753. PMLR, 17-23 Jul 2022.
[24] Francesco Montagna, Nicoletta Noceti, Lorenzo Rosasco, Kun Zhang, and Francesco Locatello. Scalable causal discovery with score matching, In 2nd Conference on Causal Learning and Reasoning, 2023. URL httpps://openrewiew net/foorumidd-6vo000
[25] Francesco Montagna, Nicoletta Noceti, Lorenzo Rosasco, Kun Zhang, and Francesco Locatello. Causal discovery with score matching on additive models with arbitrary noise. In 2nd Conference on Causal Learning and Reasoning. 2023. URL https:://pperreviee. ffrrmmdd
rV00Bx90deu
[26] Atalanti A Mastakouri, Bernhard Schölkopf, and Dominik Janzing, Necessary and sufficient conditions for causal feature selection in time series with latent common causes, In International Conference on Machine Learning, pages 7502-7511. 7. PMLR, 2021.
[27] Andreas Gerhardus and Jakob Runge. High-recall causal discovery for autocorrelated time series with latent confounders. Advances in Neural Information Processing Systems, 33: 13:15-15-2 2020.
[28] Atalanti Mastakouri. Bernhard Schölkopf. and Dominik Janzing. Selecting causal brain features with a single conditional independence test per feature. Advances in Neural Information Processing Systems, 32, 2019.
[29] Niklas Pfister, Peter Bühlmann, and Jonas Peters. Invariant causal prediction for sequential data. Journal of the American Statistical Association, 1144527):1264-11 1276, 2019.
[30] Jakob Runge, Sebastian Bathiany, Erik Bollt, Gustau Camps-Valls, Dim Coumou, Ethan Deyle, Clark Glymour, Marlene Kretschmer, Miguel D Mahecha, Jordi Muñoz-Marr, al. Inferring causation from time series in earth system sciences. Nature communications, 10(1):1-13, 2019.
[31] Jakob Runge, Peer Nowack, Marlene Kretschmer, Seth Flaxman, and Dino Sejdinovic. Detecting and quantifying causal associations in large nonlinear time series datasets. Science Advances, (1 1)::aa44996, 2019.
[32] Krzysztof Chalupka, Pietro Perona, and Frederick Eberhardt. Fast conditional independence test for vector variables with large sample sizes. ArXiv, abs/ abs//880.00777 08747727777 2018.
[33] Daniel Malinsky and Peter Spirtes. Causal structure learning from multivariate time series in settings with unmeasured confounding. In Proceedings of 2018 ACM SIGKDD Workshop on Causal Disocvery, volume 92 of Proceedings of Machine Learning Research, pages 23-47, 2018.
[34] Dehua Cheng, Mohammad Taha Bahadori, and Yan Liu. Fblg: simple and effective approach for temporal dependence discovery from time series data. In KDD '14, 2014.
[35] Jonas Peters, Dominik Janzing, and Bernhard Schölkopf. Causal inference on time series using structural equation models. Advances in Neural Information Processing Systems, 07 2012.
[36] Alessio Moneta, Nadine Chlaß, Doris Entner, and Patrik Hoyer. Causal search in structural vector autoregressive models. In NIPS Mini-Symposium on Causality in Time Series, pages 95-114. PMLR, 2011.
[37] Doris Entner and Patrik O Hoyer. On causal discovery from time series data using FCL. Probabilistic graphical models, pages 121-128, 2010.
[38] Sébastien Lachapelle, Philippe Brouillard, Tristan Deleu, and Simon bacoste-Julien Gradient- based neural dag learning. In International Conference on Learning Representations, 2020. URL htsttaae//ppeerevieeeeeee///oore wwt.:eee/EEorumm.dd....AAAAAAA
[39]
Xun Zheng, Bryon Aragam, Pradeep K Ravikumar, and Eric P Xing. Dags with no tears: Contin- uous optimization for structure learning. In S. Bengio, H. Wallach, H. Larochelle, K. Grauman, N. Cesa-Bianchi, and R. Garnett, editors, Advances in Neural Information Processing Sys- tems, volume 31. Curran Associates, Inc., 2018 URL http:: httpep///oroooccciiggsss neurips. cc/ pdf,
paper_ htt.:///eppppeeeee///////////aaa
[40] Paul Erdos and Alfred Renyi. On the evolution of random graphs. Publ. Math. Inst. Hungary. Acad, Sci., 5:17-61, 1960.
12

---

[41] Albert-Laszlo Barabasi and Reka Albert. Emergence of scaling in random networks, Sci- ence, 286(5439):509-512, 1999. doi: 10.11166/scienoe.286. 285.5300.0.22254555500 URL http://www
[42] Ulrik Brandes, Marco Gaertler, and Dorothea Wagner. Experiments on graph clustering algorithms. In Giuseppe Di Battista and Uri Zwick, editors, Algorithmss ESA 2003, pages 568-579, Berlin, Heidelberg, 2003, Springer Berlin Heidelberg. ISBN 978-3-540-39658-11
[43] D. Koller and N. Fricdman. Probabilistic Graphical Models: Principles and Techniques. Adaptive computation and machine learning. MIT Press, 2009. ISBN 9780262013192. URL
https://books www.geoooog.....o. co in/bookss htt.boookssdddddd/pptttttaae
[44] Shohei Shimizu, Takanori Inazumi, Yasuhiro Sogawa, Aapo Hyvarinen, Yoshinobu Kawahara, Takashi Washio, Patrik Hoyer, and Kenneth Bollen. Directlingam: A direct method for learning a linear non-gaussian structural equation model. Journal of Machine Learning Research, 12, 01 2011.
[45] Pedro Sanchez, Xiao Liu, Alison Q O'Neil, and Sotirios A. Tsaftaris. Diffusion models for causal discovery via topologcall ordering. The Eleventh International Conference on Learning Representations, 2023. URL h/ttts:///ppneiii oppnreeveee net/forum?idd Idd-
[46] Peter Bühlmann, Jonas Peters, and Jan Ernest. CAM: Causal additive models, high-dimensional order search and penalized regression. The Annals of Statistics, 42(6). dec 2014. URL
https://doi..
/tttp:://oi..or/1....
[47] Jiji Zhang. Causal reasoning with ancestral graphs. Journal of Machine Learning Research, (47):1437-1474, 2008. http:
wwwpggaaprrs/vvee wttpeerrsvv///agg0ge5 www.g/aa.nnna.a. html
[48] Xun Zheng, Chen Dan, Bryon Aragam, Pradeep Ravikumar, and Eric Xing. Learning sparse nonparametric dags. In Silvia Chiappa and Roberto Calandra, editors, Proceedings of the Twenty Third International Conference on Arlificial Intelligence and Statistics, volume 108 of Proceedings of Machine Learning Research, pages 3414-3425. PMLR, 26-28 Aug 2020. URL https://proceediin
mml. ppwss/vv00//aiiie 20100////en20aa... html
[49] Alexander G. Reisach, Myriam Tami, Christof Seiler, Antoine Chambaz, and Sebastian Weich- wald. Simple sorting criteria help find the causal order additive noise models, 2023.
[50] Alexander Gilbert Reisach, Christof Seiler, and Sebastian Weichwald. Beware of the simulated DAG! causal discovery benchmarks may be easy to game. In A. Beygelzimer, Y. Dauphin, P. Liang, and Wortman Vaughan, editors, Advances in Neural Information Processing Systems, httpa::///peeeeeeeeeeee www.eeewwwweeeeeeuueeeeeiie t
2021. URL CtR.httpss://
[51] Thomas Verma and Judea Pearl. Causal networks: Semantics and expressiveness. Ross D. Shachter, Tod S. Levitt, Laveen N. Kanal, and John F. Lemmer, editors, Uncertainty in Artificial Intelligence, volume of Machine Intelligence and Pattern Recognition, pages 69-76. North- Holland, 1990. doi: https:///oi.org// http:://6biiioor0..11011116100.0
URL https
com/ http://c.......taeet.teettt////a
[52] Kun Zhang, Jonas Peters, Dominik Janzing, and Bernhard Schölkopf. Kernel-based conditional independence test and application in causal discovery. In Proceedings of the Twenty-Seventh Conference on Uncertainty in Artificial Intelligence, UAI' 11, pages 804-813, Arlington, Vir- ginia, USA, July 2011. AUAI Press.
[53] Anthony C. Constantinou. Evaluating structure learning algorithms with a balanced scoring function. CoRR, abs/ 905. 12666, 2019. http:a/a://rrr//.rrrr///////////
12666
[54] Jiji Zhang and Peter Spirtes. Strong faithfulness and uniform consistency in causal inference. In Proceedings of the Nineteenth Conference on Uncertainty in Artificial Intelligence, UAI'03, page 632-639, San Francisco, CA. USA, 2002. Morgan Kaufmann Publishers Inc. ISBN 0127056645.
13

---

Contents
Introduction
2
The causal model
Problem definition
2.2
Identifiable models
3
Experimental design
3.1
Datasets
3.1.1
Misspecified scenarios
3.1.2
Data generation
3.2
Methods
3.2.1
Deepdive SCORE, NoGAM and DifTAN
Key experimental results and analysis
4.1
Can current methods infer causality when assumptions the data are violated?
4.1.1
Discussion on PC and GES performance
4. 4.1.2
Discussion score matching robustness
2
Is the choice of statistical estimators neutral
10
5
Conclusion
10
Acknowledgements
10
Assumptions connecting causal and statistical properties of the data
15
A.1 Global Markov Property
15
A.2 Causal faithfulness
16
A.3 Causal sufficiency
16
B
Details on the synthetic data generation
16
Nonlinear causal mechanisms
16
2
Non-Gaussian distribution of the noise terms
16
B.3
Algorithms for random graphs simulation
16
B.4
Modeling of unfaithful distributions
17
B.5
Dataset configurations
17
Benchmark methods
18
CAM
18
C.2
RESIT
19
GraN-DAG
19
4
Direct
19
C.5
PC
20
14

---

C.6 GES
20
C.7 SCORE
20
C.8 NoGAM
C.9 DAS
21
C. 10 Random Baseline
21
D Metrics definition
22
Possible generalisation of NoGAM the PNL model
22
Proof of Proposition
22
Discussion
23
Proof of Proposition]
23
Tuning of the hyperparameters in the experiments
23
Stability with respect to hyperparameters choices
24
Other experimental results Erdos-Renyi graphs
25
The effect of non-iid distribution of the data
25
1.2
Experiments under arbitrary distribution of the noise terms
25
1.3
Experiments with score-sortability
25
J
Experiments on SF, GRP, and FC graphs
28
J.1
Experiments Scale-free graphs
28
J.2
Experiments on fully connected graphs
28
J.3
Experiments on GRP graphs
28
K
Other results
29
FCI experiments on confounded graphs
31
M Deep-dive in PC and GES experimental results
33
A
Assumptions connecting causal and statistical properties of the data
In this section, we describe detail several crucial assumptions for causal discovery.
A.1 Global Markov Property
Causal discovery from pure observations requires assumptions that connect the joint distribution of the data with their underlying causal structure, The Markov factorization of the distribution (2) allows interpreting conditional independencies of the graph induced by the model M as conditional independencies of the joint distribution px- This is known as the Global Markov Property of the distribution Px with respect the graph G.
Definition 2. A distribution Px satisfies the Global Markov Property with respect to DAG G if:
L Xs + XA px XB Xs.
(15)
15

---

In order to simulate nonlinear causal mechanisms of an additive noise model, we sample functions from Gaussian process, such that Vi = .,..., d, (XPA). = No, ((X PA XPA))). a mul- tivariate normal distribution centered at zero and with covariance matrix the Gaussian kernel with XB. Xs disjoint subsets of X. و denoting d-separation in the graph g, and Il px denoting independency in the joint distribution px.
A.2 Causal faithfulness
A distribution Px that satisfies the Global Markov Property, decomposes according to the Markov factorization of Equation 2 [5T]. If the inverse holds, then we can consider the conditional inde- pendencies observed in the distribution px be valid conditional independencies in the graph g:
X xx X
Xs X cXB Xs.
(16)
If (16) is satisfied, say that px is faithfil to the causal graph.
A.3
Causal sufficiency
Another fundamental assumption is the absence of unmeasured common causes in the graph. Re- ichenbach principle B defines a connection between statistical and causal associations. The principle states that given the statistical association between two variables X,Y, then there exists variable Z that causally influences both explaining all the dependence such that conditioning on Z makes them independent. Causal sufficiency assumes that Z coincides with one between X and Y: resorting to
the model Equation ( this means that for each pair X,.X, there are no latent common causes.
Under the assumption of causal sufficiency of the graph and faithful distribution. we can use condi- tional independence testing to infer the Markov Equivalence Class of the causal graph G from the data.
B Details on the synthetic data generation
B.1 Nonlinear causal mechanisms
K :PAAA;PAA)) where XPA, are the observations of the parents of the node X,.
B.2 Non-Gaussian distribution of the noise terms
We generate data with non- non-Gaussian noise terms as follows: for each node {L...., d}, we model noise terms following a Gaussian distribution U, ~ o0,, with variance ~ (0.5, 1.0). Those noise terms are then transformed via a random nonlinear function f, s.t. U; t(U,). In our experiments, we sampled three different functions f, modeled as multilayer perceptrons (MLPs) with 100 nodes in the single hidden layer, sigmoid activation functions, and weights sampled from U--a, ax), a E (0.5, 3.0), respectively (c.f. Fig. 2)
B.3 Algorithms for random graphs simulation
We use four random graph generation algorithms for sampling the ground truth causal structure of each dataset. In particular, we consider the Erdos-Renyi (ER) model, which allows specifying the number of nodes d and the average number of connections per node (or, alternatively, the probability p of connecting each pair of nodes). In ER graphs, pair of nodes have the same probability of being connected. Scale-free graphs are generated under a preferential attachment procedure [4]]. such that nodes with higher degree are more likely to be connected with a new node, allowing for the presence of hubs (i.e. high degree nodes) in the graphs. The Gaussian Random Partition model (GRP) 42 is created by connecting i subgraphs (i.e. partitions) generated by ER model. A parameter Pin specifies the probability of connecting a pair of nodes in the same partition, while Pout defines the probability of connections among distinct partitions. Clusters appear when Pin >> Pout (e.g. in our experiments we consider Pi.. 0.4, Pout = 0.05). Finally, we consider Fully Connected graphs, generated by sampling a topological order and connecting all nodes in the graph to their successors with a directed edge. Given a ground truth fully connected graph, the accuracy of the inference procedure is maximally sensitive to errors in the order.
16

---

U
1.50
0.05
1.25
0.05
0.00
1.25
1.50
1.75
X
(a) (-0.5, 0.5)
U
Ü
0.1
0.1
3
0.0
0.0
-5
0
5
5
0
5
3
5
X
(b) U U-5,,, 1.5)
U
t
Ü
10
0.1
0.10
5
0.05
0.0
0.00
5
-5
10
X
X
X
(c) (U3000,3300)
Figure 2: Gaussian noise (left) transformed via random nonlinear functions (center) to non- non-Gaussian iid noise (right). Weights of the MLP sampled from either dal U--0.5 0.5), (b) (15-1,,, 1.5), or @ ( U(-3.0,3.0),
B.4 Modeling of unfaithful distributions
Given a ground truth causal graph, we model an unfaithful distribution of the data by enforeing the cancellation of directed causal effects between pairs of nodes. In practice, we identify the fully connected triplets of nodes X; - X + X2 X, in the ground truth, and we adjust the causal mechanisms such that the direct effect of X, on X₁ cancels out. In order to clarify the implementation details of our model, we consider a graph corresponding to the fully connected graph with trivial topological order æ = (X,, X2, X3}. We allow for mixed linear and nonlinear edges with liner effect on the nodes, such the set of structural equations is defined as:
with vertices X1, X2, X3 and with the set of edges
Xi = U1,
X; = f(X1) U2,
(17)
X3 = f(Xi) - X2 Us.
with / nonlinear function. This definition of the mechanisms on X3 cancels out f(X,) in
structural equation. In the case of large graphs with the number of nodes d {5, 10, 20, 50} that use in our experiments, we verify the unfaithful independencies in the data via kernel-based test of conditional independence [52], in correspondence of the pairs of nodes whose causal effect cancels out. We use a threshold of 0.05 for the conditional independence testing.
B.5 Dataset configurations
In this section, we extend the discussion of Section 3.1.2 which presents overview of the parameters that define the different configurations for the generation of the synthetic datasets of our benchmark. We sample the ground truth structures according to four different algorithms for random graph
17

---

nodes
10 nodes
20 nodes
50 nodes
Sparse = 0.1*
m = 1 m
m 2
Dense p 0.4
2
8
Graphs are re-sampled such that they have at least edges.
Table 2:
Density schema for Erdos-Ree Renyi graphs. The parameter denotes the probability of an edge between each pair of nodes the graph, and m denotes the average number of edges for each node in the graph. We scale the parameter m with the number of nodes, such that the relative density (sparsity) is similar for all graph dimensions.
generation, and according to different specifications of density, number of nodes, and distribution of the noise terms. In the case of Erdos-Renyí (ER) generated graphs, define the density of the edges relative to the number of nodes, according to the schema defined in Table For the Scale-free (SF) model, we define the edge density in the graphs according to the same values of Table but we do not generate SF graphs of 5 nodes. Similarly, we generate fully connected (FC) and Gaussian Random Partition (GRP) graphs only for {10, 20, 50} nodes. FC generation does require specifying any parameter for the density. In the case of GRP graphs, we use Pin = 0.4 as the probability of edges between pair of nodes inside the same cluster, and Pout 0.1 as the probability of edges between pair of nodes belonging to different clusters.
For each of the graph configurations, we sample ground truth and dataset of observations generated according to one of the following scenarios (described in detail in Section 3.1.1.
Vanilla additive noise model.
PNL model, with invertible post nonlinear function g(x) 3 for each speficied structural equation.
LiNGAM model, where the number of structural equations with linear mechanisms is parametrized by 8 E {0.33, 0.66, 1.03. The first two values of 5 allow modeling mixed linear and nonlinear causal mechanisms, with respectively 33% and 66% of the structural equations being linear. Unless differently specified, we consider 8 = when referring the LiNGAM model.
Confounded model, where the number of confounded pairs is parametrized by P {0.1, 0.2), denoting the probability of two nodes having a common cause. Unless dif- ferently specified, we consider p 0.2 when referring to the confounded model.
Measurement error model, where the amount of variance explained by the additive error is parametrized by 7 E {0.2, 0.4, 0.6, 0.8), denoting the inverse signal to noise ratio y == Varic.] Var[X.] with X; and e; defined the structural equation (6). Unless differently specified, we consider y = 0.8 when referring to the measurement error model.
Unfaithful model, as discussed in the Appendix B.4
Autoregressive model, defined according to the structural equation in order to simluate non-iid samples in the data.
For each scenario, and for each parametrization that it admits, we generate a ground truth G according each of the graph configurations specified at the beginning of the section, and a corresponding pair of datasets 'D of size 100 and 1000. The dataset generation is repeated under four possible distributions of the noise terms. In particular, each dataset has exogenous variables that are either normally distributed, or following a randomly generated non- Gan-Gausiian distribution, as discussed in Appendix B.2 Finally, in order ensure statistically significant results, for each pair of graph and dataset conligurations we generate G, D according to 20 different random seeds.
C
Benchmark methods
C.1
CAM
CAM algorithm 46 infers a causal graph from data generated by an additive Gaussian noise model. First, it infers the topological ordering by finding the permutation of the graph nodes corresponding
18

---

to the fully connected graph that maximizes the logllikel of the data. After inference of the topological ordering, a pruning step is done by variable selection with regression. In particular, for each variable X, CAM fits a generalized additive model using as covariates all the predecessor of X; in the ordering, and performs hypothesis testing to select relevant parent variables. This known as the CAM-pruning algorithm. For graphs with size strietly larger than 20 nodes, the authors of CAM propose an additional preliminary edge selection step, known as Preliminary Neighbours Search (PNS): given order æ,, variable selection is performed by fitting for each j 1,... d an additive model of Xj versus all the other variables {X; X, > X; in J and choosing the K most important predictor variables possible parents of X3. This preliminary search step allows scaling CAM pruning to graphs of large dimensions. In our experiments, -AM-pruning is implemented with the preliminary neighbours search only for graphs of size 50, with K = 20.
C.2 RESIT
In RESIT (regression with subsequent independence test) B the authors exploit the independence of the noise terms under causal sufficiency to identify the topological order of the graph. For each
variable X,, they define the residuals R₁ = Xx E [X; X \ [X+]], such that for a leaf node X₁ holds that R2 U²- ElUi]. The method is based on the property that under causal sufficiency, the noise variables are independent of all the preceding variables: after estimating the residuals from the data, it identifies a leaf in the graph by finding the residual R that is unconditionally independent of any node X,, vi * in the graph. Once an order is given, they select a subset of the edges admitted by the fully connected graph encoding of the ordering. We implement this final step with CAM-pruning.
C.3 GraN-DAG
GraN-DAG 38] defines a continuous constrained optimization problem to infer the causal graph from an ANM with Gaussian noise terms. For each variable X, in the graph, the authors estimate the parameters of the conditional distribution p(X, | X (X.}) with neural network ø,. They define an adjacency matrix A E Rdxd representation of the causal DAG, by finding inactive paths in the neural network computations, where path is defined as the sequence of weights of the network from the input, j to the output k: if zero weights are encountered in the path, then the output k is independent of the input j. If this repeated for all paths from j to and for all outputs k. then all paths are inactive and X; is independent of X, conditional on the remaining variables, meaning that A. = O. Note that GraN-DAG in principle does not require a post-pree post-processing consisting of an edge selection procedure, given that its output is not fully connected encoding of topological order, but an arbitrarily sparse graph. However, in practice, it is the case that GraN-DAG output approximates a fully connected graph, with large number of false positives with respect to the ground truth edges (see Table of Appendix A.3 in Lachapelle et al. [38] for quantitative results). To account for this, the authors of the method apply the our experiments, in order to compare the goodness of the ordering encoded by the GraN- output before applying the CA M-M--uunnn step, sample one order random between those admitted by the output adjacency matrix and compute its FNR-#. Given that the order is selected at random, we consider unbiased solution, we show in Figure
step top of their neural network graph output. In
C.4 DirectLiNGAM
ICA-LiNGAM 44 formulates a causal discovery algorithm for the identifiable LiNGAM model, assuming linear mechanisms and non- Gaussian noise terms. The idea is that solving for X the system defined in Equation one obtains
X = AU,
(18)
where A = (I-B)- By standard linear ICA (independent component analysis) it is possible to find A, which is equivalent to finding the weighted adjacency matrix B. This intuition lies at the base of the DirectL LiNGAM algorithm [4]. variation of ICA-LiNGAM that uses pairwise independence measures to find the topological order of the graph, and covariance-basee regression to find the connection strengths in the matrix B.
19

---

12nden
uvtrege
Figure 3: In order to evaluate the goodness of the inferred ordering of GraN-DAG, we sample one topological order at random between those admitted by the adjacency matrix before the CAM-pruning step. In this figure. compare the empirical FNR-i of order randomly sampled between those admitted by the output, against average of the FNR-& computed on the set all possible orderings admitted by the output. We see that selecting an order at random gives an unbiased representation of the average order accuracy, between those admitted by GraN-DAC DAG output before the CAM-pruning The violin plots refer to the FNR-# evaluated on ER graphs with 10 nodes over 20 different random seeds.
C.5 PC
PC algorithm (Section 5 of Spirtes et al. (4) is causal discovery method based on conditional independence testing that finds a CPDAG from the data. First, starts from a fully connected undirected graph, and estimates the skeleton of the graph by removing edges between each pair of nodes X,, X, if it finds a subset c X X,, such that IL S. Then, it finds all the v-structures X, + X3 Xx along with their directions. Finally, additional orientation rules applied to direct as many edges as possible in the output CPDAG. In our experiments, we use a kernel-basee test of conditional independence [52]
C.6 GES
The GES algorithm (9] (Greedy Equivalent Search) defines discrete optimization problem over the space of all CPDAGs, and outputs the graph that maximizes the fit measure according to some score (e.g. the Bayesian Information Criterion (BIC) score). The algorithm defined as a two steps greedy procedure. It starts from empty graph, and in the forward step it adds directed edges one by one, each time selecting the directed edge that most increases the score. When edge addition doesn't improve the score any further, in the backward step it removes edges one by one until the score stops increasing. The DAG defined by this procedure is then transformed into a CPDAG, such that GES final output Markov equivalence class.
C.7 SCORE
Rolland et al. (33 defines a formal criterion for the identification of the causal order of a graph underlying an additive noise models with Gaussian distribution of the noise terms U2 (N00 a.). Under these assumptions, the score entry leaf node X₁ is st (X)
((AA) of
then easy to verify that (X) = such that the diagonal entry of the score's Jacobian associated to leaf node is constant. Based on this relation, formal criterion identifying leaf nodes holds:
Lemma 1 (Lemma of (23)). Let X be random vector generated according to an identifiable ANM with exogenous noise terms U, (N, of), and let X, X. Then
Var [0x, 8; (X) = X, a leaf, Vi = 1,.,,.,,dd
(19)
The Lemma is exploited by SCORE algorithm for estimation of the topological order, given dataset of i.i.d. observations matrix of the score Jxx(X)) via score matching [22]. Then, it identifies leaf in the graph as the argmin, (ar[0x,s(X) which removed from the graph and assigned position in the order vector. By iteratively repeating this two steps procedure up to the source nodes, all variables in X end up
Baxd, first estimates the diagonal elements of the Jacobian
20

---

being assigned position in the causal ordering. Finally, SCORE applies the CAM- algorithm select subset of the edges in the fully connected DAG encoding of the inferred topologiaal order.
C.8 NoGAM
Montagna et al. [25] proposes generalization of SCORE, defining a formal criterion for the identification of leaf nodes in graph induced by additive noise model without restrictions on the distribution of the noise terms. After some manipulations, can be shown that the score entry of leaf X, defined in Equation ( satisfies
st sx(X) = de, log (U()).
(20)
such that observations of the pair (U₁, 81 (X)) can be used to learn a predictor score entry. For an additive noise model, the authors show that the noise term of leaf is equal to the residual defined as:
R X₁ E [X, X \ Xi].
(21)
Then, possible to find a consistent approximator of the score entry of a leaf node using R& as the only predictor.
Lemma 2 (Lemma of (25)). Let X be random vector generated according to an identifiable
ANM, and let X; E X. Then
E EE s. (X)
R. - -x))) - X; is a leaf.
Similarly to SCORE, NoGAM algorithm defines procedure for estimation of the topological order by iterative identification of leaf nodes, which are found as the argmin; E E [xx(X) R.] 8; (x))
In practice, the residuals R2, 1..... d, can be estimated by any regression algorithm, whereas the score is approximated by score matching with Stein identity (22).
C.9 DAS
Montagna et al. 24 defines condition on the Jacobian of the score function that identifies the edges of the graph induced by an additive noise model with Gaussian distribution the noise terms, given a valid causal order.
Lemma 3 (Lemma of (24)). Let X be a random vector generated according an identifiable ANM with exogenous noise terms U+ ~ NO, of), and let X E X be a leaf node. Then:
E 0100,, (X)]] A + X, E PX((), Vj [I,...,d] \ 10)
(22)
In practice, off-diagonal elements of the Jacobian matrix contain information about conditional independencies of the variables in the DAG, such that they define a condition of identification of the graph edges. Given the ordering procedure of SCORE, DAS (acronym for Discovery At Scale) defines an algorithm that can map the score function to unique causal graph with directed edges: in
practice, condition [22) of Lemma 3 is verified via hypothesis testing for the mean equals to zero.
Note that, despite the fact that this approach provides consistent estimator of the causal graph, the authors of DAS retain a CAM.-pruning step on top of their edge selection procedure based on Lemma 3 in order reduce the number of false positives of the inferred output. The benefit of DAS preliminary edge detection is that it reduces the computational costs of CAM-pruning, which is cubic in the number of nodes in the graph, such that it doesn't scale well to high dimensional graphs. Overall, given an input dataset X E RR with n number of samples and d number of nodes in the graph, DAS computational complexity is O(d + d²), whereas, for the SCORE algorithm this is Odn" + nd
C.10 Random Baseline
In our experimental analysis of Section we consider the performance of a random baseline in terms of F1 score and FNR- accuracy of the order (Figure [ Our random baseline is defined as follows. Given a graph with d variables, we sample random topological order π as permutation of the vector of elements XI...., Xa. Then, given the fully connected graph admitted by the order, with the set of edges& & (X,, + Xx, Xx, < X Vi,j = 1...., d). for each pair of connected nodes we sample Bernoulli random variable with parameter p = 0.5, such that the edge is removed for Y 0.
21

---

D Metrics definition
For the evaluation of the experimental results of our benchmark, we consider the Fl score, the false positive (FP) and false negative (FN) rates of the inferred graph, and the false negative rate FNR-f of the fully connected encoding of the output topological order. In order to specify the Fl score, we need a definition of FP, FN, and true positive (TP), that applies to both undirected and directed edges, given that we evaluate both DAGs and CPDAGs.
We define as TP any predicted edge that is in the skeleton of the ground truth graph (i.e. the set of edges that doesn't take direction into account).
We define the FPs the edges in the skeleton of the predicted graph that not in the skeleton the ground truth graph. Note that this definition of FP doesn't penalize undirected edges edges inferred with reversed direction.
We define as FN pair of nodes that are disconnected in the predicted skeleton while being connected the ground truth. Additionally, we count as false negatives inferred edges whose direction is reversed with respect to the DAG ground truth.
Then, the Fl score is defined as the ratio TPFO
E Possible generalisation of NoGAM to the PNL model
Proposition of Section 3.2.1 suggests that it is possible to generalize Lemmapand accordingly, the NoGAM algorithm, to the case the post nonlinear model.
E.1 Proof of Proposition
Proof. Let X E Rd be random vector generated by the post nonlinear model of Equation (5). Given the Markov factorization of Equation (2). the logarithm of the joint distribution px satisfies the following equation:
log
log Px, (X.).
Then, for a node X, in the graph the score entry s, is defined according to Equation (8) whereas given a leaf node X, in the graph, si satisfies the following:
(x) = 0x, logpx(X) 0x, logP)(X))PAi).
(23)
Our goal is to show that 0x log pt (XI(PA)) = 0x, pogi) with Uj = (X)) - and g the inverse of the postnonlinear function g (which is invertible by modeling assumption). As notational remark, in what follows we will drop any sub-index on the distribution of the random variables, which we distinguish by their argument. Also, we denote realizations of random variables (or random vectors) with lowercase letters (e.g. 2; is the value of the random variable Xt). We rewrite the distribution of X₁ conditional on its parents by marginalizing over all values of Up:
p(a: pa,) =
p(r: pay, ())(())))))
(24)
[ p(ri ()pm)) (2) = g(fi(pa;) + u)) du
(25)
[ p(11 pay, ()p(()) [(w - - fi f(
(26)
with being the indicator function. Being an invertible function, the value of u equals to
9- (xz) fi (pa,) is unique, which implies that p(z; pa,, u)) 0 if i * g (ax) fi (( else
p(x) pay, u)) = 1. Let us denote ut := g (xx) )))) - fi(pa,). Then, the integral in Equation 26 simply
becomes:
p(æ1 pay) = / dp(u) = ut) = p(u)).
(27)
Thus, 0x, log (ogp(() 0x, log p(Ui.
22

---

E.2 Discussion
Proposition derives a connection between Lemma definne by Montagna et al. 25 for identifiable additive noise models to the case of a PNL model. Note that the authors define a consistent estimator of st score function of leaf node X; from the residual R; X E [Xx X \ XiJ. which satisfies R₁ = Ut in the case of an ANM with noise terms centered at zero. In general, the latter equality does not hold for post nonlinear model, meaning that regression of leaf variable against all the other variables of X does not guarantee consistent estimation of the disturbance on the leaf structural equation. This implies that, as is, NoGAM doesn't provide theoretical guarantees of consistent estimation of the topological order of a PNL model.
F Proof of Proposition 2
We define two lemmas preliminary to the proof of Proposition 2
Lemma 4. Let X E Rd be generated according an SCM M that satisfies score-sortability, and let G be the graph induced by the model. Then, there exists leaf node of 9 that is score crore-identffaale.
Proof. By contradiction, let's say that the node X₁ with / argmin, Varls, (X)) is not leaf node. Then, the causal order ı where X1 is a successor of all other nodes in the graph, is not correct ordering, implying that the model is not score-sortable.
Lemma 5. Let X E Rd be generated according an SCM M that satisfies score-sortability, and let be the graph induced by the model. Let My( the model defined removing the leaf node X₃ from the set of structural equations of M. Then, the model M. is score-sortable.
Proof. By contradiction, let's assume that Mi does not satisfy score-sortability. such that the node Xm with m = a...
Var[s;(X]] not a leaf node in the graph 93 induced by My Then, any topological order with Xm successor of all nodes 1,1-1. a.1-1,.... I,..., d, is a wrong topological ordering of the graph g. This implies that
M is not score-sortable.
Now, we present the proof of Proposition 2
Proof. (Proof of Proposition 2 By Lemma
4 being M a score-sortable model, there exists a leaf X1 such that == argmin, Var 88 [)). Then,
Varlo, log (X; PA.)] Var(a, (X, +
Σ
VarlOxlogpa [AA)]] C,
KECH,
for all 1,..., d. Moreover, by previous Lemma 5 the model M, score-sortable. Thus there exists an index m E (1, ( 1,
d) such that Xm, is leaf and Var[sm (Xy1))] S
Var[s;(X, (1) )]. Vi 1, .......6 1, 1, ..... d. Then, the topological ordering defined by iterative identification of leaf nodes with Lemma and Lemma 5 on the subgraphs resulting by removal of a leaf node, is correct with respect the model M.
G Tuning of the hyperparameters in the experiments
The methods included in the benchmark require the tuning of several hyperparameters for the inference procedure. In particular, PC, DAS, SCORE, NoGAM, RESIT, GraN-DAG, CAM, and DiffAN require a threshold a over the p-value of the statistical test used for the edge selection procedure. Instead, GES applies a regularization term weighted by its score, which penalizes the number of edges included in the inferred graph: the higher the value of X, the sparser the solution. Given that the tuning of both a and À requires prior knowledge about the sparsity of the ground truth, there is no established procedure for finding their optimal values in real-world settings, where the ground truth is not accessible. Thus, in order to enable fair comparison between all the methods, we always select the optimal value of a E {0.001, 0.01, 0.05, 0.1) and À (0.05, 0.5, 2, 5) over each benchmark dataset. Section we discuss the stability of the algorithms with respect to choices of these hyperparameters.
23

---

des
cam
Score
nogam
diffen
grandag
resit
PC
9e5
varilla BR-5-4parse
vanilla ER-1O-sparse
1.0
vanilla ER-20-sparse
vanilla ER.50-Gparse
3.8
68
3.6
3.0
varilla BR-S-daanse
ER-10-dense
vanill ER.20-dense
Vanilla
02
02
02
02
The violin plots in the figure represent the difference between the Fl score of a method running inference with hyperparameters optimized using the ground truth, versus the FI score of the same method using a default value of the hyperparameters. We denote this difference with [fl fll. In the case of GES, we define as 0.05. The violin datasets and graphs generated according to 20 different random seeds.
Figure 4:
default À = 0.5. For all the remaining methods, the default alpha threshold is defined as plots refer to the inference performance Results in the table are on data generated from the vanilla scenario. and we consider Erdos Kendy graphs with the number nodes in {5, 10, 20, 50) in the dense and sparse settings.
GraN-DAG and DiffAN both define a learning procedure over the data, which requires the tuning of several training hyperparameters, the most important of which is the learning rate 1). For each dataset, this is optimized over the loss function on a held-out validation set, without accessing the ground truth graph.
H Stability with respect to hyperparameters choices
Most causal discovery methods come with hyperparameters that alleviate minor assumption violations (e.g. sparsity regularization or higher thresholds on p-values in statistical tests). In the absence of background knowledge, tuning these hyperparameters is an art that often relies on pre-conceptions about reasonable solutions. In this section, we investigate the impact of these hyperparameters on the accuracy of the output graph. GES penalizes dense solution with regularization term in its score, weighted by a hyperparameter À that can not be tuned in the absence of the ground truth. Similarly, an threshold on p-values for statistical tests for edge selection is required by all the benchmarked methods (excluding GES and DirectLiNGAM) and can not be tuned without knowledge of the ground truth. In this section, we analyze the inference Fl score by fixing a and À to the commonly accepted default values of 0.05 and 0.5 respectively. In Figure we summarise the absolute value of the difference between the Fl score obtained with hyperparameters optimized on the ground truth of each dataset, against the F1 score yielded by inference with the default a and À values (denoted with fili in the plots). According our empirical findings, the case of graphs with at least 10 nodes, the median of this difference is in general lower than 0.1, and most of the time close to 0. Sparse graphs seem to be more affected in their performance by the hyperparameters choice: this means that using the default a and À causes an increase of false positives in the output graph. Bühlmann et al. [46] shows that under correct topological order, a graph whose set of edges a superset of the ground truth still provides consistent estimates of the causal effects, such that increasing the false positives doesn't t affect the outcome of downstream tasks, but only the statistical efficiency of the inference. Given that estimation of the topological ordering is not affected by the choice of a and À values, we suggest that the role of hyperparameters value this respect marginal with respect to the task of interest.
24

---

Implications. We observe remarkable stability of the benchmarked methods with respect to the choice of their hyperparameters. The biggest drops in Fl score are observed on sparse graphs. meaning that the default parameters cause an increase of false positives, which nevertheless does not affect the downstream task of interest of consistent estimation of causal effects.
I Other experimental results on Erdos-Renyi graphs
In this section, we present additional experimental results on Erdos-Renyi graphs.
L.1 The effect of non-iid distribution of the data
Figure TF (right) illustrates that all the methods included in our benchmark do not perform well on data sampled from non-iid distribution generated according to the autoregressive model of Equation 77) Fl score and FNR- are indeed similar to that of the random baseline. It clearly appears that none of the presented algorithms provide guarantees of good empirical performance in the setting of non-iid distribution of the data.
1.2 Experiments under arbitrary distribution of the noise terms
In Section 2.2 we discussed the effect of the distribution of the noise terms on the identifiability of the causal graph underlying an SCM. Given that the assumption of Gaussian distribution of the disturbances is often not satisfied in real datasets, it is important to provide empirical evidence on the performance of the benchmarked methods on data generated with an arbitrary distribution of the noise. In this section, we discuss experiments on data generated with the noise terms that are iid samples from the distribution of Figure 2c Similar to Section 4 we analyze results on ER graphs with 20 nodes, with experiments repeated over 20 random seeds. In this section, we include results of DirectLiNGAM, both linear and nonlinear SCMs.
Figure 3 iillustrates the FNR-7 score of the inferred topological order on data generated according to the vanilla model with non -on-Gaasiinn noise terms. Under these conditions, NoGAM and RESIT provide theoretical guarantees of consistent estimate of the causal ordering. Similarly, PC and GES do not make explicit assumptions on the distribution of the noise terms (despite the fact that GES optimizes a Gaussian likelihood). SCORE, DiFFAN, DAS, and CAM instead are limited by restrictions on the noise terms, which are required to be normally distributed. However, Figure 5 (right) shows that, except for CAM, they can estimate the order with accuracy comparable to that achieved in the case of vanilla generated data, with Gaussian distribution of the disturbances. These observations are in line with the experimental findings in Montagna et al. [25), which shows how the structure of the score entries of leaf nodes can still be exploited by SCORE for inference on data generated under arbitrary noise distribution. Our experimental results agree with this intuition: surprisingly, SCORE ordering ensures better FNR-t accuracy than RESIT, despite the latter being explicitly designed to be insensitive to the distribution of the noise terms. Interestingly, we notice that the median of the violin plot referred to DirectLiNGAM in Figure (right) is close that of RESIT and CAM: this suggests that in the realistic scenario of mixed linear and nonlinear mechanisms with non-Gaussian additive disturbances, we expect Directl.in .iGAM to give performance significantly (left), shows that the
better than several methods designed to perform on nonlinear ANM. in the case of methods whose ordering accuracy is comparable to the Gaussian case, the F1 score after pruning is also comparable to that on Gaussian data. This means that CAM-pruning is robust with respect arbitrary distributions of the noise terms. Additional experimental results on data generated according to the misspecified scenarios of Section 3.1.1 with non-Gaussian distribution of the disturbances, are presented in Figure 6
Implications. Most of the benchmarked methods are capable of robust inference on datasets generated by ANM with non-Gaussian noise terms. DirectL.iNGAM shows remarkable performance, comparable to that of several methods designed for inference on nonlinear additive noise models.
I.3 Experiments with score-sortability
In this section, we present the experimental results of a simple ordering algorithm, that we name ScoreSort, based on the score-sortability criterion defined in Section 4.1.2 Given the random vector
25

---

1 rdods
d86
com
Score
nogam
diffan
grandag
resit
lingam
pC
ges
gaussian
non gon-gausiin
Vanilla (non- Gaussian)
ER-20 dense
F1 score
FNR-
1.0
1.0
0.8
0.8
0.6
0.6
0.4
0.4
0.2
0.2
0.0
0.0
Figure 5: F1 score and FNR-i on data generated with non-Gaussian distribution of the noise (c.f. Figure 2c For each method, we also display the violin plot of its performance on the vanilla scenario with Gaussian noise terms, with transparent color. Fl score (the higher the better) and FNR-f (the lower the better) are evaluated over 20 seeds Erdos-Renyi dense graphs with 20 nodes (ER-20 dense). FNR-f is not computed for GES and PC methods, whose output is CPDAG.
lingam
gaussian
non-gaussia
model (non- Gaussian)
20 derse
LINGAM (non- Gaussian)
ER dense
score
FNR-
08
0.2
02
(a) PNL
(b) LiNGAM
Confounded model (non-Gaussian) ER-20 dense
Measurement nt error model Gounsanus ER-20 dense
F1 score
FNB-6
F1 score
FNII-f
1.0
1.0
0.6
0.6
6.2
3.0
8.0
0.0
(c) Latent confounders
(d) Measurement error
Unfaithfll model Inon-Gausse (ann-Gassiin)
ER-20 dense
Autoregressive model (non- Gaus-aaas
ER-20 dense
FNR.6
F1
FNR-i
0.8
08
(e) Unfaithful distribution
(f) Non-iid data distrihution
Figure 6: Experimental results on the misspecified scenarios with non- non-Gaussian distribution of the noise terms. For each method, also display the violin plot of its performance on the same misspecified scenario under Gaussian distribution of the noise terms, with transparent color. Fl score (the higher the better) and FNR-# (the lower the better) are evaluated over 20 seeds on Erdos-Renyi dense graphs with 20 nodes (ER-20 20 dense). FNR-# is not computed for GES and PC, methods whose output is a CPDAG.
26

---

nogam
Varilla Ei--3 denwe
model FR-20
LIGAM FR-2 Confounded model- ER-00 dense
3.2
PNR-F
HM--
1.5
PMl-
PAR-é
644
0.4
044
0.4
08
(a) Vanilla
(b) PNL
(c) LiNGAM
(d) Latent confounders
ressiceee
0.6
0.6
(e) Measurement error
(I) Unfaithful distribution
(g) Non-iid distribution
Figure 7: Experiments on score-sortabillty We compare the FNR-î accuracy of the simple ScoreSort baseline (c.f. Algorithm I with SCORE and NoGAM algorithms performance. The violin plots are evaluated over 20 seeds Pedos--Renyi dense graphs with 20 nodes (ER-20 dense).
X E Rd generated according to a structural causal model, the ScoreSort baseline identifies the index of a leaf node as the argmin,; Var[s;[X)] Then it removes the leaf node X₁ from the set of vertices of the graph, and identifies the next leaf with the argmin the variance of the score vector of the remaining (sub)graphs obtained by the iterative leaves removal yields a topological ordering æ for the graph underlying the SCM. If the model is score-sortable, then, according Proposition 2 the causal order π is correct with respect to the graph. Details on ScoreSort are presented in the Algorithm box [ In practice, ScoreSort estimates the score vector according to the same score-matching algorithm used by SCORE and NoGAM, which is based on the Stein identity (22).
of nodes. Identification of leaf nodes according this procedure over the
Figure 7 compares ScoreSort performance with NoGAM and SCORE ordering algorithms, data generated according to the vanilla and misspecified scenarios of Section [31.1 In line with considerations in the discussion that ScoreSort FNR-f accuracy is comparable, with statistical significance, to that of SCORE and NoGAM. This is true both in the case of data generated under vanilla and misspecified scenarios.
matching robustness in Section 4.1.2 we observe
Algorithm ScoreSort algorithm for inference of the causal order
Input: data matrix X B
+ []
nodes + [1.
for i = 1....,d do
$ mcore-matching(X)
Lindez argmin Var(s]
1 5- mrdesllinde..
+ [|
Remove linden-th column from X; Remove / from nodes end for
и + reverse(T) (first node a source, last node is a leaf)
return .
27

---

Experiments on SE, GRP, and FC graphs
In this section, we analyze the Fl score and FNR-0 accuracy of the benchmarked methods in the case when the ground truth graph is generated with Scale-free, Fully Connected, and Gaussian Random Partitions algorithms.
J.1 Experiments on Scale-free graphs
Figure 8 illustrates the F1 score and FNR-# on SF graphs. We see that similar to the case of ER networks, score matching- bbased methods show remarkable robustness in the inferred order in the case of several misspecified scenarios, particularly, on data generated by the PNL (Figure (a measurement error (Figure 8d right), and unfaithful models (Figure Se right). However, we notice two significant differences with respect to the conclusions that we derived in the case of ER graphs in Section Figure T in the case the LiNGAM model, SCORE, DAS and NoGAM display FNR-1 accuracy that is remarkably close to that vanilla data (Figure 8b rbight,,), whereas their decrease in performance in the case of latent confounders effects (Figure rightt, is worse than that observed on ER graphs. Interestingly, the results the F1 score show that DAS, SCORE, NoGAM, and DiffAN performance is surprisingly good (with respect to the random baseline) across all the misspecified scenarios, which suggests good performance of CAM-pruning on SF graphs. Moreover, we see that GraN-DAG and RESIT inference procedure is close to that of the random baseline in almost all the misspecified scenarios: this is also explained by the poor performance of these two methods on vanilla data and SF graphs (illustrated in the transparent violin plots of Figure (8)
based approaches show remarkable robustness even in the case of SF graphs. Interestingly, CAM-prunnng performance on SF graphs is generally better than the one relative to ER- -eggeenrate ground truths, such that the observed Fl score is often better than random. We also observe that RESIT and GraN-DAG ordering ability is negatively affected by the SF ground truth, in comparison the case of ER graphs.
Implications.
J.2 Experiments on fully connected graphs
In the case of fully connected graphs, the ground truth admits a unique topological ordering. This means that expect to observe increase in the false negative rate FNR-r, with respect to the results on ER graphs of Figure [ This is in line with our empirical evidence, as illustrated in Figure D However, we see that score matching- baannn--aaedddee approaches still show robust performance in the inference of the ordering with respect to misspecified scenarios, except for the case of data generated according to the LiNGAM ground truth model. Notably, the Fl score accuracy of GES consistently better than that of all the other methods, across every scenario. This is to be understood with the fact that the unpenalized BIC score optimized by GES always improves by increasing the number of edges in the graph. Given that we optimize the regularizer term À on each dataset, the optimal À value will naturally privilege the densest solutions, Different is the case for methods that rely on CAM- Maru-runnng. which display an Fl score consistently lower than the random baseline, except for the case of data generated by the unfaithful and LiNGAM models.
Implications. Score matching- based approaches are in general robust to misspecifications of the scenario in the case of fully connected ground truth. GES shows remarkable performance, that is partly explained by the optimization of the loss penalization term directly on the ground truth. Finally, we observe that the CAM-pruning step is negatively affected by the large density of the graphs.
J.3 Experiments on GRP graphs
In Figure 10 we see that score matching bbasdd methods and CAM algorithm display better robustness in the inference of the order than the remaining approaches, in reference all of the benchmarked scenarios. The FNR-# of RESIT, GraN-DAG, and DiffAN are significantly close to the random baseline for data generated according most of the ground truth models (with the exception of DiffAN on the LinGAM model and GraN-DAG on unfaithful samples). In terms of F1 score, most of the methods show good capability of inferring the ground truth graph, even in the case of data generated under assumption violations. Note that the Fl score of the random ground truth is remarkably bad, if compared the case of SF, FC, and ER graphs. This is line with the cluster structure of GRP graphs: given that the random baseline connects pair of nodes all with the same
28

---

das
Cam
nogam
diffr
grandag
rest
pC
ges
misspecried
PNL model SF-20 dense
UIGAAM- SF- 20 dense
F1 score
FNR-R
F1 Foore
FNR-T
1.0
1.0
0.
0.6
0.4
1A
0.4
00
(a) PNL
(b) LiNGAM
Confounded model-SF-22 dense
Measurement model SF-20 dense
1.0
2.0
0.0
0.2
0.
Latent confounders
(d) Measurement error
Unfaithful model
SF-20 dense
hutoregressive model SF-20 den
3.5
F1 score
1.
PNR-6
1.0
F1 score
1.0
PNBR-f
0.0
0.6
E
6.2
1.0
2.0
0.0
(c) Unfaithful distribution
(f) Non-iid data distribution
Figure 8: Experimental results on the misspecified scenarios. For each method, we also display the violin plot of its performance on the vanilla scenario with transparent color. Fl score (the higher the better) and FNR-f - (the lower the better) are evaluated over 20 seeds on Scale- free dense graphs with 20 nodes (SF-20 dense). FNR-6 - is not computed for GES and PC methods, whose output CPDAG. Note that DirectLiNGAM performance does appear, as both the linear mechanisms and non-Gaussian noise assumptions are violated.
probability 0.5, we expect large number of false positives due to edges between nodes of different clusters.
Implications. Score matching- natching-based approaches and CAM algorithm are remarkably robust to model misspecification both in terms of F1 score and FNR- accuracy.
K
Other results
Statistical efficiency. Figure 11 shows F1 score and FNR-f accuracy on datasets with sample size equals to 100 and 1000. Comparing the relative difference in performance with respect to different sample sizes, we get empirical idea of the statistical efficiency of the inference methods. In line with our expectations, the experimental results show that both metries are negatively affected by the reduction in sample size. Interestingly, in the case of SCORE, DAS, NoGAM, and DirectLiNGAM iiGAMM, we observe better stability of the FNR- compared to the other methods, with the score matching- based approaches that are in general significantly better than the random baseline also with datasets of size 100.
The effect of the graph size and density. Figure 1 illustrates the F1 score and the FNR- accuracy on datasets generated according to the vanilla scenario and ground truth graphs that differ in size and density. In particular, we consider the case of dense and sparse graphs, with {5, 10, 20, 50} nodes. Interestingly, we see good stability of the Fl score across different graph dimensions in the sparse case. The decrease in performance due to larger graph sizes is more evident in the case of dense graphs: this is particularly true for dense graphs with 50 nodes, where the preliminary neighbours search step (described in Section C.I before the CAM-pruning reduces the ability to
29

---

anden
das
Cam
Score
nogam
diffr
grandag
rest
pC
ges
misspecrfedd
PNL madel FC
LINGAM FC:20
F1 score
FNR-R
Floce
FNR-T
1.0
0.
0.
E
0.4
TA
SA
0.4
5.0
0.
(a) PNL
(b) LINGAM
Confounded model FC-20
Measurement me del FC-20
0.0
0.0
1
0.2
12
4
0.2
02
(c) Latent confounders
(d) Measurement error
Unfaithful model - FC-20
Autoregressive model FC-20
3.5
F1 score
1A
PNR-6
1.0
F1 score
1.0
FNBI-F
0.0
0.
0.6
6.2
1.0
2.0
0.0
(c) Unfaithful distribution
(f) Non-iid data distribution
Figure 9: Experimental results on the misspecified scenarios. For each method, we also display the violin plot of its performance on the vanilla scenario with transparent color. Fl score (the higher the better) and FNR- (the lower the better) are evaluated over 20 seeds on fully connected graphs with 20 nodes (FC-20). FNR- is not computed for GES and PC methods, whose output CPDAG. Note that DirectLiNGAM performance does appear, as both the linear mechanisms and non-Gaussian noise assumptions are violated.
infer true positives for most of the methods. Considering the FNR-# of the inferred orders, that, similarly to what we observed in the analysis of the Fl score, in the case of sparse ground truths most of the methods display stable results across different graph dimensions. Indeed, score matching- based approaches, as well as CAM, DiffAN, GraN-DAG, and DirectLiNGGAM do not display any clear evidence of degraded performance for larger graphs. In the dense setting, instead, we see that CAM and DirectLiNGAM accuracy in the inference of the order is negatively affected by larger dimensionality.
The balanced scoring function. Figure 12 illustrates the inference accuracy of misspecified scenarios in terms of the Balanced Scoring Function (BSF), proposed by Constantinou [53] This is defined as:
1 TP TN FP
FN
where: TP and FP denote the true and false positives, respectively, and TN and FN are the true and false negatives; a and represent the number of ares and independencies in the true graph respectively. The key difference between the BSF and the F1 score is given by the fact that the balanced scoring function accounts for the whole confusion matrix (i.e. TP, FP, TN, FN), whereas in the F1 score the true negatives are not included in the computation. The BSF ranges from -1 (worst) to 1 (best), with a value of 0 corresponding to the score of an empty graph. correlated to the F1 score. The main difference that we see is that the Fl score in general penalizes PC and GES performance more than the BSF (when comparing their accuracy to the random baseline), meaning that they tend to infer graphs with low true positive rates and large true negative rates.
general, we expected the BSF to be
30

---

Cam
nogam
diffr
grandag
rest
pC
ges
misspecried
PNL model GRP-20
LINGAM GRP-20
F1 score
FNR-R
F1
FNR-T
1.0
1.0
0.6
0.4
7.
14
0.4
0.0
(a) PNL
(b) LiNGAM
Confounded model GRP-20
Measurement model GRP
0.0
0.2
12
0.2
Latent confounders
(d) Measurement error
Unfaithful model GRP-20
Autoregressive model GRP-20
1.4
F1 score
1.0
PNR-6
1.0
F1 score
1.0
FNB--f
2.0
0.0
0.6
E
6.2
1.0
3.0
0.0
(c) Unfaithful distribution
(f) Non-iid data distribution
Figure 10: Experimental results on the misspecified scenarios. For each method, we also display the violin plot of its performance on the vanilla scenario with transparent color. Fl score (the higher the better) and FNR-A t (the lower the better) are evaluated over 20 seeds on Gaussian Random Partitions graphs with 20 nodes (GRP-20), FNR- is not computed for GES and PC methods, whose output is CPDAG. Note that DirectLiNGAM performance does not appear, both the linear mechanisms and non-Gaussian noise assumptions are violated.
L FCI experiments on confounded graphs
In this section, we describe the experimental setting for the FCI algorithm (Fast Causal Inference) 4 Given that the method can handle latent confounders, we focus experiments on data generated from graphs admitting latent common causes.
PAG. The output graphical object of FCI is a Partial Ancestral Graph (PAG) 47 It admits six types of edges. We denote the two ends of an edge as marks. The possible marks are tail (-), circle (c), and an arrowhead (>), which combined allow for six edges. These graphs represent an equivalence class for Maximal Ancestral Graphs, which are graphical objects that represents the presence of confounders effects and selection bias [2]
Metrics. For the evaluation of the FCI inferred output, we adopt the strategy proposed by Heinze- Deml et al. [17] (see their Section 4.2). We define true positives, false positives. and false negatives over three possible adjacency matrices, each one defined by a specific query.
IsPotentialParent query: the estimated adjacency matrix has Ay = 1 if there is an edge between X; 1 Xj. X, -x X; + Xy. X, o-x,. X2 a-o Xj. X, X, in the estimated PAG, else A; = 0. A.j 1 denotes the case in which X, is a potential parent of Xj-
IsAncestor query: the estimated adjacency matrix has A, = 1 if there is path from X; to X, with edges of type X; -ox,,, X; X,, X, o-x, X, in the estimated PAG, else Ag 0. Ag = denotes the case in which X; is an ancestor of Xj.
PPPtentialAnces/n query: the estimated adjacency matrix has Asj -1 if there a path from X₁ to X, with edges of type X; Xy. X₁ -oX,, X; + X,. X, o- Xg. X, o-o X,,
31

---

rdddos
das
Score
diffn
grandag
resit
lingam
100 samples
1000 samples
FR-10-sparse
score ERo20-sparse
score FR-S0-sparse
1.0
1.0
1.0
!
1
0.2
02
F] Score ER-S-dense
F1 Score ER-10-dense
F1 score ER-20-dense
F1 Score ER-SO-dense
0.4
8
0.0
0.0
0.0
PNR-h ER-S-sparse
FNR-G- ER-10 1--10-s
FNR--- ER-20-sparse
FNRR-G ER-SO-sparse
1.0
1.0
0.6
0.6
0.6
0.
0.2
FNR.A B..---RR-ddee..
FNR-8- E----ER------- ER-10-dense
FNR: E--EE-2--deee
EE-ERR-O0-. ER-50-dence
0.0
0.0
0.4
0.4
0.4
0.0
0.0
Figure 11: Experiments on the effect the sample size. We compare
score and FNR-i accuracy on datasets generated under the vanilla scenario with Gaussian noise, with different sample sizes. We remark that in the case of DirectLiNGAM, in order to provide meaningful results, we report the performance on datasets with non-Gaussian noise terms. Violin plots filled with color refer to datasets of size 100, and transparent violin plots refer to datasets of size 1000. The metrics reported on Prdos-Renyyi graphs of size (5, 10, 20, 50) both in the sparse and dense case (PC GES are not included graphs 50 nodes, as their computational demand is too high). Experiments are repeated over 20 different random seeds.
32

---

nogam
diffr
grandag
pC
ges
misspecriedd
PNL model
ER-20 dense
LINGAM
ER dense
Fl score
BSF
Fl score
BSF
1.0
1.0
10
3.0
0.5
(a) PNL
(b) LiNGAM
Confounded model- ER-20 dense
Measurement model ER-20 dense
BSF
25
0.5
-15
0.2
-05
Latent confounders
(d) Measurement error
Unfaithful model
ER-20 dense
Autoregressive model ER-20 dense
Fl score
BSF
F1 score
BSF
1.0
1.0
1.5
0.0
0.0
0.0
0.0
-10
(c) Unfaithful distribution
(f) Non-iid data distrihution
Figure 12: Experimental results on the misspecified scenarios. For each method, we also display the violin plot of its performance on the vanilla scenario with transparent color. BSF and F1 score (the higher the better) are evaluated over 20 seeds on Erdos Renyi dense graphs with 20 nodes (ER-20& dense).
Erdos-Renyi
X,0+ X, in the estimated PAG, else Ay- = A.j denotes the case which X; is potential ancestor of Xz-
For each adjacency matrix defined by one of the three queries, we define true positives, false negatives, and false positives follows:
A true positive (TP) is pair i.j with Au in both the inferred and ground truth adjacency matrices (with the ground truth DAG converted to PAG).
A false negative (FN) is pair i with Ai 0 in the inferred matrix, and Ai in the ground truth (with the ground truth DAG converted to a PAG).
A false positive (FP) is pair i, with A;j in the inferred matrix, and A, = 0 in the ground truth (with the ground truth DAG converted to PAG).
Given these definitions of TP, FN, FP, we define the F1 score as F1 TPHOMFP4FNN which we to present empirical results in Figure
TP
M Deep-dive in PC and GES experimental results
In this section, we analyze GES and PC performance in terms of false positive and false negative rates on graphs characterized by different numbers of nodes and edge density. Section 4.1.]] we discuss the case of inference with PC and GES on andös-Renyl dense and large graphs (20 nodes): Figure [ in the main text reports PC and GES Fl score be consistently and significantly worse than random across all of the tested scenarios. We argue that this is in line with previous findings in the literature
33

---

das
cam
nogam
grandag
lingam
pc
ges
ER-S-pparse
sparsee
sparse
soore
50 50-sparse
0.6
06
0.0
onte
F1 score ER:5 do-sedense
F1 score ER-1edense
F1 score ER-20-dense
F1 score ER 50 derse
1.0
0.8
0.6
0.5
02
0.2
52
FNR-G-ER-BO ER-S-aarse
FNR ER-10sparse
FNR-A- EN-20ssaasse
FAR-A BN-50-sparse
10-
10
10
0.
0.4
0.4
0.4
044
00
drten
FNR--- FNR---ER--S-E ER-S-dense
PAR-I- ER-10-derse
FNR-I- ER-20-dense
FNR-f- ER-S0-dense
0.6
0.6
0.6
0.6
02
0.2
02
Experiments the effect of the graph size and graph density. We compare the Fl and FNR-i accuracy on datasets generated under the vanilla scenario with Gaussian noise, on ground truth graphs with the number of nodes {5, 10, 20, 50) both in the sparse and dense case (PC and GES are not included for graphs of 50 nodes, as their computational demand is too high). We remark that in the case of DirectLiNGAM, in order to provide meaningful results, report the performance datasets with non- non- Gaussian noise terms. The metries reported Erdos-Remt Renoi graphs. Experiments repeated over 20 different random seeds.
Figure 13:
34

---

IsPotentis/Parent sporse
IsAncestor- ER sparse
IsFotentaalAncsstor sparse
1.
1.0
20
number of nodes
number of
Rumber of nodes
IsPotentislParent ER denie
IsAntestor- ER dente
IsPatentallAncistor - ER den
64
13
number nodes
number of
number nodes
Figure 14: FCT performance on dense and sparse ER graphs, on datasets generated under latent confounders effects.
diffan
pransag
lingam
missprec/iee iTeed
venilla
PNL model ER-20 derse
LINGAM ER-20 dense
FNR
FPR
1.0
FNR
1.0
FPR
6.8
6.1
0.0
(a) PNL
(b) LiNGAM
Confounded model ER-20 dense
Measurement error ER dense
1.0
0.6
7.5
0.
1.2
(c) Latent confounders
(d) Measurement error
Unfeithfull model - ER-20 dense
Autoregressive model ER-20 dense
FNR
FPR
3.0
FNR
FPR
0.4
2.A
2.A
0i
(c) Unfaithful distribution
(f) Non-iid data distribution
Figure 15: ENR (False Negative Rate) and FPR (False Positive Rate) of the experiments on the misspecified
scenarios, on Erdos-Reny. dense graphs with 20 nodes (ER-20 dense). For each method, we also display the violin plot of its performance on the vanilla scenario with transparent color. The noise terms normally distributed, except for the LiNGAM model, in which case we generate disturbances according to non-Gaussian distribution.
by Uhler et al. 77. which shows that the set of distributions that are not strong-faithff ([54) has non-zero Lebeasgue measure, in contrast to the set of unfaithful distributions that has zero measure (which justifies the common assumption of faithful causal models). In particular, the measure of the set of not strong- stronn-faithful distributions tends to increase for large and dense causal graphs, which we argue is key explain the degraded performance of PC in the ER-20 dense graphs of Figure [ Coherently with our claim, Figures [9]and] and 20 show FNR and FPR consistently better than random for both PC and GES, respectively on sparse and dense graphs of 5 nodes.
35

---

random
nogam
grancag
Irgam
ges
misspocifldd
vanilla
ER-20 sparse
LINGAM ER-20
08
O8
52
02
(a) PNL
(b) LiNGAM
Confounded model ER sparse
Meassrement error model ER sparse
PNR
FMR
RNR
FAR
3.0
10
2.0
10
0.4
5A
5.0
0.
(c) Latent confounders
(d) Measurement error
Unfaithful model ER-20
Autoregressive model ER-20
0.8
08
0.2
02
02
Unfaithful distribution
(f) Non-iid data distribution
Figure 16: FNR (False Negative Rate) and FPR (False Positive Rate) of the experiments the misspecified scenarios, Erdos-Remt Renyi sparse graphs with 20 nodes (ER-20 sparse). For each method, also display violin plot of its performance on the vanilla scenario distributed, except for the LiNGAM model, in which case we generate disturbances according to non- non-Gaussian distribution.
transparent color. noise terms are normally
36

---

random
nogam
grancag
Irgam
ges
misspoctflod
vanilla
model
LINGAM ER.
0.8
08
0.2
0.
(a) PNL
(b) LiNGAM
Confounded model ER-10 dense
Measurement error model ER dense
3.2
PNR
10
FPR
2.0
RNR
10
FAR
0.4
2.4
0.4
0.0
(c) Latent confounders
(d) Measurement error
Unfaithful model ER-10 dense
Autoregressive model ER-10 dense
0.08
08
0.2
2.2
1.2
02
Unfaithful distribution
(f) Non-iid data distribution
Figure 17: FNR (False Negative Rate) and FPR (False Positive Rate) of the experiments the misspecified scenarios, Frdos-Renyi dense graphs with 10 nodes (ER-10 dense). For each method, also display the violin plot of performance on the vanilla scenario distributed, except for the LiNGAM model, in which case we generate disturbances according to non- Gaussas distribution.
transparent color. noise terms are normally
37

---

random
nogam
grancag
Irgam
ges
misspociflod
vanilla
ER-10 sparse
LINGAM ER-10
08
08
0.2
02
(a) PNL
(b) LiNGAM
Confounded model ER-10 sparse
Micasurement error model ER-10 sparse
3.0
PNR
10
FMR
2.0
RNR
10
FAR
0.4
2.4
(c) Latent confounders
(d) Measurement error
Unfaithful model ER-10
Autoregressive model ER/10
0.88
08
0.2
12
02
Unfaithful distribution
(f) Non-iid data distribution
Figure 18: FNR (False Negative Rate) and FPR (False Positive Rate) of the experiments the misspecified scenarios, Erdos-Reny' sparse graphs with 10 nodes (ER-10 sparse). For each method, also display violin plot of its performance on the vanilla scenario distributed, except for the LiNGAM model, in which case we generate disturbances according to non Gan-Gassi distribution.
transparent color. noise terms are normally
38

---

diffen
misspeclfied
vanilla
PNL model ER-5 dense
LINGAM ER-S dense
FNR
FPR
RIR
FPR
1.0
64
2.4
0.4
0.0
(a) PNL
(b) LiNGAM
Confounded
0.
0.2
0.2
(c) Latent confounders
(d) Measurement error
Unfalthful model-ERe ER-S dense
Autoregressive model ER-5 dense
FNR
PPR
FNR
FPR
1.0
2.A
0.4
0.0
(e) Unfaithful distribution
(l) Non-iid data distribution
Figure 19: FNR (False Negative Rate) and FPR (False Positive Rate) of the experiments on the misspecified scenarios, on Erdos-Renyi dense graphs with 5 nodes (ER-5 dense). For each method, we also display the violin plot of its performance on the vanilla scenario with transparent color. The noise terms normally distributed, except for the LiNGAM model, in which
generate disturbances according non Gaussian distribution.
39

---

random
d5s
dittan
granzogg
restt
lingom
misspeciTed
varilla
PNL model
ER-5 sparse
UNGAM ER-S sparse
FPR
3.2
1.0
1.0
1.
0.
6.2
50
2.0
0.0
(a) PNL
(b) LiNGAM
Confounded model ER-S
Measurement model-
ER-5
FNR
FPR
FPR
0.3
3.
0.0
0.0
Latent confounders
(d) Measurement error
Unfaithful model ER-5 sparse
Autoregressive model ER-5 sper se
FPR
2.6
0.5
0.6
02
50
00
(c) Unfaithful distribution
(f) Non-iid data distrihution
Figure 20: FNR (False Negative Rate) and FPR (False Positive Rate) of the experiments the misspecified scenarios, on Erdos-Renyi sparse graphs with 5 nodes (ER-5 sparse). For each method, we also display the violin plot of its performance on the vanilla scenario with transparent color. The noise terms are normally distributed, except for the LiNGAM model, in which case we generate disturbances according to non- non-Gaussian distribution.
40


Reviews:
[{"review_id": "ZTPnK3oLp9", "reviewer": "Reviewer_QK5D", "summary": "The paper considers the important problem of what happens to the reliability of the causal model output of different types of causal discovery algorithms when the underlying assumptions behind the method are violated. It finds that score matching-based methods tend to be more robust and provides a theoretical explanation.\n", "questions": "(see also remarks above)\n- why do you try to emphasise the claim that ‘score matching procedures are *robust*’ whereas if anything in Fig.1 the FNR-\\pi performance of RESIT is almost always nearly identical between the vanilla and ‘misspecified’ cases, contrary to the score matching algorithms that do often show a drop in FNR-\\pi? (so RESIT may not be the ‘best’, but it is the most *robust* one)\n\nOther remarks (for lack of a better place):\n- typo title ‘casual’\n- 28: ‘These strong conditions are arguably necessary’ => No: there are a whole range of algorithms that can handle confounders (or even cycles) these days, and the FCI benchmark algorithm has been around for over 30 years now.\n- 57: in one sentence you mention 12 papers that are irrelevant to yours ….\n- Table 1: GES optimises the Gaussian likelihood => no: GES is a search strategy. It incorporates a scoring metric, but that does not need to be the Gaussian likelihood. (Although it might explain the observed performance).\n- 203: ‘no straightforward comparison between PAGs and CPDAGs’ => ?? both show edges and edge marks, and indicate when they don’t know … \n- 214: ‘differences .. mostly in the FNR-\\pi score’ => I have no idea why that should be the case, especially since you say all start from the same skeleton procedure? Also in FIg.1 RESIT shows no change in FNR-\\pi even though it does show significant differences in F1. Most importantly, however, is that it is very difficult as a reader to interpret this score in terms of how ‘good’ or ‘robust’ the actual causal model output is.", "limitations": "NA", "rating": 3, "confidence": 4, "soundness": 2, "presentation": 3, "contribution": 2, "strengths": "- very important and challenging problem that is often encountered in practice\n- good overview of causal discovery methods and common assumptions behind them\n- highly commendable attempt at providing a large suite of software + data sets for future standardised benchmark evaluations on (new) causal discovery algorithms\n", "weaknesses": "- Unfortunately, despite the laudable goal, the comparison is too limited and biased to be meaningful.\nThe generated models all have nonlinear interactions and/or non-Gaussian noise … which is a key assumption for e.g. score match-based methods. But then you should also see how they hold up in the standard linear Gaussian case (answer: terrible) before you start proclaiming their ‘remarkable robustness’. \nSimilarly, if ‘additive noise’ is key, then compare against e.g. multiplicative noise (not ‘post nonlinear’, as that can often be treated in the original way). ‘Faithfulness’ relates to independence based methods (CAM/PC/GES), so ‘unfaithful’ is not a misspecified model for the score match-based algorithms … so hardly ‘remarkable’ they are unaffected by that scenario. No consideration of acyclic vs. cyclic (in practice nearly all relevant systems have feedback mechanisms). \n\n- trying to get to a universal claim / quantification on what is ‘the’ most robust method is doomed to failure as it is an ill-posed problem. Instead what you can meaningfully do is try to assess how the performance of individual methods degrades as the degree of violation of assumptions increases: is there graceful degradation or sudden collapse beyond a certain threshold? Is a method more sensitive to some types of violations than others? Is the entire model ‘bad’, or is the skeleton still reliable and mainly the orientations that suffer? etc.\n\n- it is bad practice to take the outcome of a single run of a causal discovery algorithm as ‘the’ answer. At least one should do a form of bootstrapping to come to a balanced model interpretation, and use that as the basis for comparison between methods. Also, the FNR-\\pi metric used in the evaluation (Fig1(a-f)-right) is hard to interpret as a measure of model quality: instead something like the structural Hamming distance may be more appropriate.\n\n- main claim that ‘score matching-based methods are most robust’ is not justified based on the evidence presented in Fig.1 (see remarks below)\n\n- reported experimental results of some methods (PC/GES) are so laughably bad (much worse than random guessing) that I have no idea what happened, but something went wrong\n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for the valuable insights provided with their comment. In the response, we will follow the structure of the bullet list in the review.\n\n-   The reviewer comments that our \"comparison is too limited and biased to be meaningful”, supporting their claim with a series of points.\n\n\t> The generated models all have nonlinear interactions and/or non-Gaussian noise …\n\n\tOur models include:\n\n\t1. Linear mechanisms\n\n\t2. ANM with nonlinear mechanisms\n\n\t3. PNL model with nonlinear mechanisms\n\n\tThe choice of the noise terms is motivated by the fact that the distinction between Gaussian and non-Gaussian is what really matters in the assumptions of the methods involved, with CAM, SCORE, DAS, GraN-DAG, DiffAN only supporting Gaussian noise, while RESIT, NoGAM, and PC supporting both Gaussian and non-Gaussian noise, and DirectLiNGAM supporting only non-Gaussian noise terms. Given that these are the distinctions to bear in mind, our choice is consistent with the assumption we are willing to probe (i.e. effect of misspecification of the noise distribution on the performance).\n\n\t> You should also see how they hold up in the standard linear Gaussian case (answer: terrible)\n\n\tAs correctly remarked, there is no reason to expect that methods based on finding the topological ordering will be robust in the setting where ordering can not be identified for a pair of nodes (see L 98-99).\n\n\t> … ‘unfaithful’ is not a misspecified model ...\n\n\tWe agree with the reviewer, which is why we only define as “_surprising_” the “_decrease in FNR-pi performance with respect to the vanilla scenario_” (L 304-305), and immediately after we explain the observed results, connecting graph sparsity with the observed performance (L309-311). As noticed by reviewer eJF4, this is in line with previous literature: [1] Costantinou et al., 2019 point out that graph sparsity directly affects the f1 score. Our observation that associates unfaithful path canceling with this fact from previous literature is thus an original, pertinent and meaningful statement.\n\n\tWe are happy to add the citation to [1] to back our claim and to remove the word _surprising_ from the main text, to enhance better clarity.\n\n\t[1] _Evaluating structure learning algorithms with a balanced scoring function_, Costantinou et al., 2019\n\n\t> _No consideration of acyclic vs. cyclic_\n\n\tCyclic graphs are beyond the scope of this paper. We explicitly rely on topological ordering for evaluating the graphs (FNR-pi metric). This is clearly not possible with cyclic graphs. Moreover, despite the fact that in principle one would like to test any possible violation, this is not possible in the scope of a single paper, where our work already tackles important questions in a large scale effort, as recognized by the reviewers (”highly commendable attempt at providing a large suite of software + data sets”, R QK5D; “experiments are extensive.”, R eJF4).\n\n-   > trying to get to a universal claim/quantification on what is ‘the’ most robust method is doomed to failure as it is an ill-posed problem. ...: is there graceful degradation or sudden collapse beyond a certain threshold?\n    \n\n\tIn this regard, we note that Appendix B.5 provides a detailed discussion on the parametrization of the models with respect to the _amount_ of violation. In Appendix B.5 we also clarify that the results are presented in the setting where the parameter maximizes the amount of violation, e.g. for measurements error, we specify that we show results for the maximal noise-signal ratio: with the measurement error example in mind, we argue that if we observe robust performance of a method in this _extreme_ setting maximizing the noise-signal ratio, we expect good/better performance to be observed when noise-signal ratio is reduced, given that the distribution of the data would get closer to that of data generated by the vanilla model. Given that finer grained analysis of the experiments is not compatible with space limits, we will add a section discussing results degradation with respect to parameter changes, as we agree with the reviewer that this would be a valuable contribution.\n\n\t> _Is a method more sensitive to some types of violations than others?_\n\n\tThis information is already in the paper, and it is actually one of its key contributions. Section 4.1 is indeed structured in paragraphs, each one analyzing which methods are more sensitive and which are more robust (in the sense of performance under model misspecification) with respect to the assumption violation introduced in a specific data generating process. Similarly, Figure 1 (a-f) provides a comprehensive overview of the effect of each type of violation on the performance of each algorithm.\n\n\t> _Is the entire model ‘bad’, or is the skeleton still reliable and mainly the orientations that suffer?_\n\n\tWe agree that this is valuable information, and we commit to add such evaluation in the appendix of the paper. We add preliminary plots in the pdf attached to the rebuttal in the common response (Figure 1.)\n\n-   > it is bad practice to take the outcome of a single run of a causal discovery algorithm as ‘the’ answer ...\n    \n\n\tWe agree that when dealing with real data in applications, the standard practice is to use some type of bootstrap to quantify the uncertainty on edges. On the other hand, in the case of algorithms' evaluation and comparison, it is a well established and the most common practice in literature to consider the outcome of a single run as the answer (see, for instance, the papers of SCORE, DAS, NoGAM, GraNDAG, CAM, to cite a few of those that we used in our work). This difference is motivated by the fact that in practical applications the user is interested in inference on a single graph. In our setting, instead, adding bootstrapping on top of 2M experiments is clearly not a feasible option.\n\n-   For the last two points, please refer to the common answer."}, {"review_id": "uxM1DkFaE1", "reviewer": "Reviewer_eJF4", "summary": "The authors present the results of an extensive empirical study evaluating the performance of algorithms for learning the structure of causal graphical models from observational data. In particular, the study examines the robustness of these algorithms to violations of various untestable assumptions that are typically made by one or more of these algorithms. These violations are Gaussian noise, non-Gaussian noise, linear mechanisms, nonlinear mechanisms, unfaithful distributions, confounding effects, and measurement errors.", "questions": "(Repeating from above)\n\n1. Why does PC generally perform *worse than random* in terms of F1. Random or near-random performance is one thing, but performance that is consistently worse than random needs to be explained. \n\n2. Why do misspecifications such as measurement error further *decrease* the F1 score of PC?  We would expect that measurement error would push the F1 score more toward random rather that further decrease performance.\n\n3. Why does misspecification sometimes *improves* other scores (e.g., RESIT FNR-pi under Autoregressive, GES F1 score under LiNGAM). In one case (RESIT FNR-pi under PNL) you do explain (273-275), but the other anomalies are not explained.", "limitations": "(none)", "rating": 6, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "The focus of the paper is extremely important and overdue. There is only very weak understanding with in the research community about how existing methods for learning causal structure respond in practice to various types of misspecification. This is particularly true for several methods that have been devised fairly recently.\n\nThe experiments are extensive. As the authors state: “Our experimental protocol consists of more than 2M experiments with 11 different causal discovery methods on more than 60000 datasets synthetically generated.”\n\nThe authors release code for both the data generation and for the methods themselves.\n\nThe authors provide a random baseline, which is very useful for contextualizing the results. In some cases (e.g., confounding, autoregressive), the condition reduces almost all methods to essentially random performance.\n\nThe results are examined in some detail, and the authors provide a variety of high-level takeaways to help readers interpret the vast number of experiments that they have performed.", "weaknesses": "The authors use false-positive and false-negative rates (and F1) as their measures of error. This ignores a great deal of recent scholarship [1, 2, 3, 4] regarding the problems with these measures and proposals of better measures. Some of the measures require model parameters [3] or knowledge of likely intervention targets and outcome measures [4], but others do not [1, 2]. The paper would be improved by examining the results of these error measures as well as the simpler ones of FPR, FNR, and F1.\n\nThe authors observations regarding the results for unfaithful models provide a specific case of one reason for using these alternative measures. In an effort to explain why faithfulness violations had little effect on most methods, they say: “…the set of edges of the graph G ̃ is sparser than that of the ground truth, due to the cancellation of causal effects. Thus, given that inference on sparser graphs is generally easier, it can positively affect the empirical performance, in line with our observations.” As Constantinou [2] points out, graph sparsity directly affects the F1 score, and his proposed BSF accounts for this. \n\nThe simplest way to summarize the results (at least those in Figure 1) is that the relative ordering of performance of the methods is maintained under misspecification. That could reflect greater robustness of the score-matching based methods. However, it could equally well reflect just overall performance or some bias in the particular error measures used to evaluate the learned structures (that’s one reason for the comment above about error measures). Despite this, the authors claim that the results demonstrate that specific methods are surprisingly robust (“Our empirical findings indicate that score matching-based methods are surprisingly capable of partial recovery of the graph structure in several of the misspecified scenarios.”).\n\nIn the “Implications” section (312-315), the authors correctly point out that “Our experimental findings show that most of the benchmarked methods significantly decrease their performance on the misspecified models. This is particularly problematic since the violations considered in this work are realistic and met on many real-world data.” This is correct and, to my mind, this is one of the primary results of the paper. Almost without exception, each of these misspecifications reduces the performance of nearly all methods. However, the authors then go on to note: “On the other hand, we observe surprising robustness in the inference of score matching-based methods.” This is much less well supported by the empirical results. The paper would be greatly improved by providing empirical results that distinguished between the hypotheses of “consistent proportional decrease from varying initial performance” and “special robustness of some methods despite their initial high performance”. I’m seeing strong evidence for the first, but not for the second.\n\nThere are at least two aspects of the results that are not sufficiently explained and should be. First, the paper should explain why PC generally performs *worse than random* in terms of F1. Random or near-random performance is one thing, but performance that is consistently worse than random needs to be explained. Related to this, why do misspecifications such as measurement error further *decrease* the F1 score of PC?  We would expect that measurement error would push the F1 score more toward random rather that further decrease performance.\n\nSecond, the authors should explain why misspecification sometimes *improves* other scores (e.g., RESIT FNR-pi under Autoregressive, GES F1 score under LiNGAM). In one case (RESIT FNR-pi under PNL) you do explain (273-275), but the other anomalies are not explained. \n\nThe results of the hypothesis tests against the random baseline should be displayed on Figure 1 somehow.\n\n*References*\n\n1. Constantinou, A. C., Liu, Y., Chobtham, K., Guo, Z., and Kitson, N. K. (2021). Large-scale empirical validation of Bayesian Network structure learning algorithms with noisy data. *International Journal of Approximate Reasoning*, vol. 131, pp. 151-188.\n\n2. Constantinou, A. C. (2019). Evaluating structure learning algorithms with a balanced scoring function. *arXiv preprint arXiv:1905.12666*.\n\n3. Gentzel, A., Garant, D., & Jensen, D. (2019). The case for evaluating causal models using interventional measures and empirical data. *Advances in Neural Information Processing Systems*, *32*.\n\n4. Peters, J., & Bühlmann, P. (2015). Structural intervention distance for evaluating causal graphs. *Neural computation*, *27*(3), 771-799.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for their effort spent in providing valuable feedback on our work. In the response, we will follow the paragraph listing structure of the review.\n\n-   > The authors use false-positive and false-negative rates (and F1) as their measures of error. This ignores a great deal of recent scholarship [1, 2, 3, 4] ...\n    \n\tWe agree that [1, 2] Balanced Score Function (BSF) proposal seems appropriate for the evaluation of causal graphs and it is appealing in the sense that requires no extra knowledge about model parameters or interventions (which are out of scope in our work). Still, we argue that F1 is a valuable choice for the evaluation of the causal graph outputs. In the rebuttal's `.pdf` file (Figure 2), we propose example plots that show a strong correlation between F1 and BSF. This is in line with the results found in Figure 1 of [2]: in the comparison between F1 and the normalized BSF, we see perfect consistency between the two scores in all of the proposed scenarios, except for one (scenario 3.3 on the x-axis). This suggests that using BSF rather than F1 score would not revert any of the conclusions about the _robustness_ of the involved methods with respect to model misspecification, which is the core of our work. We support our claim with empirical evidence by comparing the BSF and the F1 score (Figure 3 in the rebuttal `.pdf`) in the case of unfaithful data (mentioned in the second paragraph of the _Weaknesses_ section by the reviewer): BSF and F1 score follows a consistent pattern in the performance, meaning that our conclusions in this setting are not biased by the choice of F1 as one of the metrics. We notice that the normalized BSF is always larger than F1, which can be explained by the inclusion of the TN in the count of the score, and that is consistent with the observations in Figure 1 of [2].\n\n\tAs the choice of the scoring function in causal discovery is far from trivial, we agree that BSF results would be a valuable addition to our work, as they provide a summary metric accounting for the entire confusion matrix. We commit to add a discussion on BSF results in the appendix.\n\n-   > The simplest way to summarize the results is that the relative ordering of performance of the methods is maintained under misspecification. That could reflect greater robustness of the score-matching based methods. However, it could equally well reflect just overall performance or some bias in the particular error measures used ...\n    \n\tPlease refer to the comment above (first point in the response)\n\n-   > ... However, the authors then go on to note: “On the other hand, we observe surprising robustness in the inference of score matching-based methods.” This is much less well supported by the empirical results. ...\n    \n\tPlease refer to the second point in the common response.\n\n-   > The paper would be greatly improved by providing empirical results that distinguished between the hypotheses of “consistent proportional decrease from varying initial performance” and “special robustness of some methods despite their initial high performance”.\n    \n\tPlease refer to the second point in the common response, especially the remark we make about the use of the word robustness: as it is common in machine learning, we use it as a synonym for performance in misspecified scenarios. If the use of the word robustness sounds confusing here, we are willing to replace it with the expression _performance under misspecification_.\n\n-   > … the paper should explain why PC generally performs _worse than random_ in terms of F1.\n    \n\tPlease refer to the first point in the common response. We commit to add the discussion proposed also in the main text, as this is indeed of interest.\n\n-   > …Related to this, why do misspecifications such as measurement error further _decrease_ the F1 score of PC? …\n    \n\tWe point the reviewer to the closed form example in Section 2 of _Causal discovery in the presence of measurement error: Identifiability conditions,_ Zhang et al., 2017 : given the causal graph X_1 ← X_2 → X_3 with variables subject to measurement error and to linear mechanisms, it is possible to show that the correlation between X_1 and X_2 vanishes as the noise-signal ratio tends to infinity (implying an increase in the false negatives rate when using correlation to test for conditional independence). Inversely, checking the correlation coefficient between X_1 and X_3 conditional on X_2, as the noise-signal ratio goes to infinity this tends to approximate the correlation coefficient between X_1 and X_2 in the absence of measurement error (i.e. we would observe an increase in the false positives rate).\n\n\tThis intuition is corroborated by the observations in Figures 14d to 18d in Appendix, which report the FNR and FPR.\n\n-   > authors should explain why misspecification sometimes _improves_ other scores (e.g., RESIT FNR-pi under Autoregressive, GES F1 score under LiNGAM).\n    \n\tIn the case of GES, the LiNGAM model satisfies the assumption of the method, thus we have no reason to expect the performance to be negatively affected.\n\n\tIn the case of RESIT in the autoregressive setting, the discrepancy in performance on vanilla and non-iid data is substantially equal to that observed in all of the other misspecified scenarios (Figure 1 b,c,d,e right), while the only remarkable gap is observed for RESIT on PNL data. The violin plot of the FNR-pi on non-iid data completely overlaps with that on vanilla data, and there is only a slight shift in the median (< 0.05), thus we do not think it is of scientific relevance as it might be due to the samples fluctuations: given that evaluation is done on newly sampled data for each scenario, we are not interested in absolute statements about each sample, but instead about statements that appear to be statistically significant with respect to our experiments."}, {"review_id": "AbEqXDjHt3", "reviewer": "Reviewer_JGo2", "summary": "The authors focus on evaluating causal structure learning algorithms.  These algorithms, in order to claim that their outputs are causal, require various assumptions, most of which can't be tested easily in observational data.  This work aims to test how robust a set of causal structure learning algorithms are to violations of these assumptions.  The authors pay particular attention to a group of additive noise model-based methods that use score matching and find that they are generally the most robust to assumption violations.", "questions": "You discuss SCORE, NoGAM, and DiffAN as being part of the same class of methods that derive constraints on the score function (in Section 3.2.1).  However, in the post-results discussion (Section 4.1.1), you discuss the group of SCORE, DAS, and NoGAM instead.  Despite operating similarly to NoGAM and SCORE, why do you think DiffAN was not as robust?  And similarly, what makes DAS comparable here, despite not being part of the original score-based discussion?\n\nFollowing on to my comment about using a Gaussian process as the nonlinear mechanism, what led to this choice?  Did you consider and reject other mechanisms, or do you have reason to believe that this one is realistic enough to be a reasonable choice?", "limitations": "In general, this work is focused on addressing limitations in the field (since most work merely specifies the assumptions without assessing empirically how necessary these assumptions are).  I think the results of this work are a bit limited, as discussed above, due to the limited diversity of methods assessed.\n\nIn terms of social impact, I only see positives to deeper exploration of the robustness (or brittleness) of methods in the face of assumption violations.", "rating": 8, "confidence": 4, "soundness": 3, "presentation": 4, "contribution": 3, "strengths": "I appreciate the authors' focus on evaluating assumptions - this is rarely done and is very important for the field.\n\nThe assumption-violation settings are well-chosen, and Table 1 provides a helpful reference.  In general, the way you chose to implement these seems reasonable.\n\nThe random baseline, as well as showing both the misspecified and vanilla results side-by-side, is a great choice and helps readability of otherwise dense graphics.", "weaknesses": "The diversity of method is lacking.  In Table 1, FCI and DirectLiNGAM are included, but neither are present in the actual evaluation.  This means that of the 9 methods presented in the evaluation, 7 are additive noise model-based methods with nearly the same assumptions.  While it says that experiments on FCI are in Appendix I.5, no section I.5 exists.  Appendix L looks like it starts discussing FCI results but never actually does so.  Similarly, the Figure 1 caption says that DirectLiNGAM results are in Appendix I.4, but again, no section I.4 exists.  I understand that it's much easier to compare algorithms that all make the same assumptions, since they can be compared against the same type of baseline.  But the result is that it's more a comparison of which additive noise-based method is more robust to assumption violations than a true comparison across structure learning methods.  I otherwise think this work is great, but for a paper focused solely on evaluation and pitched as a comparison across different types of methods (since the conclusion is that score matching methods are best), it really feels like it's lacking substantial representation from those other types of methods.\n\nI understand you're already manipulating a lot of settings in combination, but I wish the baseline datasets were more functionally varied.  From what it looks like, since all of the methods (apart from LiNGAM, which was omitted from the evaluation) require nonlinearity, you defined that nonlinear mechanism as a Gaussian process.  However, there are other nonlinear mechanisms that exist, and some additional variations there would help round out the analysis.\n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "\nWe thank the reviewer for considering our work and for their valuable comments and insights.\n\n-   > The diversity of method is lacking. In Table 1, FCI and DirectLiNGAM are included, but neither are present in the actual evaluation.\n    \n\n\tIn the case of FCI, when confounders are not present in the model, FCI is reduced to simple PC. The only interesting information comes in the case of data generated under confounders effects: in this setting, FNR, FPR, and F1 definitions of section 3.2 (paragraph starting at L204) don’t allow for straightforward comparison of PAGs with DAGs and CPDAGs. This is why we rely on the metrics used in Heinze et al., 2017. In this sense, a comparison with the other methods’ robustness is not discussed, as we would be comparing different metrics, which is not immediate.\n\n\tIn the case of DirectLiNGAM, we omit the results from the main manuscript as most of the data are generated from nonlinear models with Gaussian noise, which is heavily incompatible with the method's assumptions (in practice we observed DirectLiNGAM performance to be significantly worse than the random baseline in this setting, which is not surprising). For this reason, we moved the discussion on DirectLiNGAM to Appendix I.2 _Experiments under arbitrary distribution of the noise terms_.\n\n\tConcerning the reviewer's remark that _diversity of methods is lacking_, we agree that our analysis is biased towards methods based on nonlinear additive noise models, and we agree to be more explicit about this fact in the camera ready version of the paper. Nevertheless, our benchmark also includes the most prominent constraint and score-based methods in literature that also have an available and well maintained Python implementation. In this regard, please consider the well established library for causal discovery _causal-learn_: the only algorithms available for constraint and score based methods are PC, FCI, GES, and additionally DirectLiNGAM. Similarly, for the _Dodiscover_ library that we used in our work. Limitations in the available Python implementations are substantial and thus reflected in our benchmark, where existing implementations were privileged. In this sense, we believe that our software contribution is an additional important value of our work.\n\n-   > While it says that experiments on FCI are in Appendix I.5, no section I.5 exists … Similarly, the Figure 1 caption says that DirectLiNGAM results are in Appendix I.4, but again, no section I.4 exists.\n    \n\n\tWe now notice that we introduced a mismatch in the section’s numbers between the file with only the main manuscript with respect to the file with the supplementary material. We are sorry for the confusion this has generated and ask please to consider appendix I.2 of the pdf with supplementary material for the discussion on DirectLiNGAM performance. The analysis of FCI is in Appendix L (as already found out by the reviewer) rather than Appendix I.5. Again, we apologize for the inconvenience.\n\n-   > you defined that nonlinear mechanism as a Gaussian process. However, there are other nonlinear mechanisms that exist, and some additional variations there would help round out the analysis.\n    \n\n\tand again, in the _Questions_ section:\n\n\t> Following on to my comment about using a Gaussian process as the nonlinear mechanism, what led to this choice? ...\n\n\tAs correctly hinted by the reviewer, this choice is mostly due to the combinatorial nature of the number of experiments (“I understand you are already manipulating a lot of settings\"), which _explodes_ each time a new _dimension_ for data generation is added. In addition to this, we point to a large amount of established and recent literature using our same proposal in the data generation (e.g. _Distinguishing cause from effect using observational data: methods and benchmarks_, Mooij et al. 2009; _Causal discovery with reinforcement learning_, Zhu et al., 2016, + the papers introducing SCORE, DAS, NoGAM, GraN-DAG, RESIT, CAM, ...).\n\n-   > … Despite operating similarly to NoGAM and SCORE, why do you think DiffAN was not as robust?\n    \n\n\tWe discuss this point in L 246-248 “… the difference is in that it [DiffAN] estimates the score function with probabilistic diffusion models, whereas SCORE, NoGAM, and DAS rely on score matching estimation”. The main novelty of DiffAN with respect to SCORE, is that it estimates the score function with probabilistic diffusion models, which appears to be less effective than score-matching in the finite samples regime. To make this point, we dedicate the content of Section 4.2, with title _Is the choice of statistical estimators neutral?_. In particular, see L 355-358.\n\n-   > ... what makes DAS comparable here, despite not being part of the original score-based discussion?\n    \n\n\tThe details of DAS are discussed in Appendix C.9, in which we clarify that DAS can be seen as an accelerated version of SCORE, where they share the same ordering procedure (based on identifying leaf nodes on the variance of the Jacobian of the score function estimated via score matching), while they differ in the pruning step (DAS also has a preliminary pruning step that exploits information in the Jacobian of the score function to reduce the computational cost of running CAM pruning afterward). We will include this discussion in the main paper in Section 3.2.1."}]

Author rebuttal:
We thank the reviewers for their effort in reading our paper also providing valuable insights. In general, there is agreement that this is _very important_ work (R QK5D , eJF4, JGO2), both in terms of the subject of the paper and in the scale of the experiments (_highly commendable attempt at providing a large suite of software + data sets for future standardized benchmark evaluations_, R QK5D; _experiments are extensive_, R eJF4), with meaningful practical impact. Moreover, there is appreciation in the presentation of the results, with the use of a random baseline that _is a great choice and helps readability_ (R eJF4, R JGO2), and in the choice of the tested violations (R JGO2).

In this common response, we address comments shared by more than one reviewer. We are responding to each reviewer independently on the remaining points.

-   R QK5D expresses concerns with respect to the performance of PC and GES methods with respect to the random baseline, mentioning “_reported experimental results of some methods (PC/GES) are so laughably bad (much worse than random guessing) that I have no idea what happened, but something went wrong_”, whereas R eJF4 suggests that the paper should explain why PC generally performs worse than random in terms of F1.

	We respectfully disagree with the comment from R QK5D that “_something went wrong_”. Figures from 14 to 18 illustrate the methods’ performance (FNR and FPR) on ER graphs with 5 and 10 nodes, both sparse and dense, and on sparse graphs of 20 nodes: under these conditions, PC and GES perform consistently better than random in the vanilla case, whereas for the misspecified scenarios, we observe performance reasonably close to the remainder of the benchmarked methods. Additionally, we remark (as already reported in the paper, footnote 1 of page 6) that we used PC and GES implementation from the well-established and maintained library _DoDiscover__._ Moreover, we report that our findings are in line with those in the literature: for example, [1, 2, 3, 4, 5, 6] all explicitly state that PC and/or GES are excluded as baselines for the evaluation of their proposed method, as they perform “_significatly worse_”. In light of these considerations, we have no reason to believe that something went wrong or that there are bugs in the code.

	Finally, we provide a more detailed explanation of the observed performance of PC and GES. In the dense and large graphs case (which is the one under consideration in Figure 1), a lot of quasi-violations of faithfulness are expected to happen (see _Geometry of the faithfulness assumption in causal inference_, Uhler et al., 2012 ). Thus, we would expect in this setting PC and GES to be systematically wrong (which is what we observe), and coherently to this, we do not expect any significant drop in performance of these two methods switching from the vanilla to the unfaithful setting (which is coherent with the observations of Figure 1e). If the reviewers agree that this is a valuable point, we are willing to include this discussion in Section 4, in order to provide better clarity with respect to the observed results for PC and GES.

-   R QK5D and R eJF4 express concern about our claim that score matching-based methods provides robust performance with respect to model misspecification. In relation to this, R QK5D adds ”_if anything in Fig.1 the FNR-\pi performance of RESIT is almost always nearly identical between the vanilla and ‘misspecified’ cases, contrary to the score matching algorithms that do often show a drop in FNR-\pi? (so RESIT may not be the ‘best’, but it is the most robust one)”_. We merge our response to the two points (score match and RESIT robustness) here, as these two points are related.

	In the case where the vanilla performance is comparable to the random baseline (which is the case for RESIT, Fig1 a-f) the method is not only _not the best_ (as remarked by the reviewer in reference to RESIT), but in general unable to make inference better than random with statistical significance in _any_ of the settings. If we don’t account for this systematic inability of the method of doing reliable inference, we should conclude that the random baseline is the most robust causal discovery algorithm under assumption violation, but this would clearly be an ill-defined notion of robustness.

	On the other hand, score matching-based methods are the only ones that show performance better than random with statistical significance across almost all of the considered data generating scenarios. In this regard, we point out that in machine learning robustness is generally intended as the ability of the model to perform well in misspecified scenarios. If the use of the word _robustness_ sounds confusing in this setting, we are happy to replace it with the expression _performance under misspecification_. In this sense, the best performance is consistently provided by score matching-based methods, especially in terms of their ability to infer the topological ordering (measured by FNR-pi score).

_References_

[1] Score matching enables causal discovery of nonlinear additive noise models, Rolland et al., 2022

[2] Causal discovery with score matching on additive models with arbitrary noise, Montagna et al., 2023

[3] Scalable causal discovery with score matching, Montagna et al., 2023

[4] CAM: Causal additive models, high-dimensional order search and penalized regression, Buhlmann et al., 2014

[5] Gradient-based neural dag learning, Lachapelle et al., 2019

[6] Dags with no tears: Continuous optimization for structure learning, Zheng et al., 2018

Meta-review: Two of the reviewers recommend acceptance and were further strengthened in their recommendations by the authors' replies.  One reviewer disagrees, so we could not reach consensus.  I am recommending acceptance based on all the comments taken into consideration in total, but the comments of that last reviewer should be addressed in the final paper by the authors.  Of particular importance are the comments in by that reviewer in reply to the author rebuttal, reproduced here:

"I am afraid my main concern still stands: laudable goal, ok method, but overstated claims based on biased evaluations. The method relies strongly on nonlinearity of interactions or non-Gaussianity of the noise. Other methods can have other assumptions (e.g. 'linear Gaussian'). If you want to claim your method is much more robust than others, you cannot present a range of scenarios where each case satisfies your key signal for inference (all scenarios involve strong nonlinearities or strongly non-Gaussian noise), but violates those of one or more others. To then claim 'look how much more robust our method is in general' is then at best 'not sufficiently supported'.

The authors acknowledge the method will perform terrible in the linear Gaussian setting. So then show how the performance behaves in systems that vary in degrees of nonlinearity + non-Gaussianity. Does the algorithm still perform well in small to moderate amounts of nonlinearity + non-Gaussianity or not? How strong do they need to be? What if parts of the model have interactions that are close to 'linear Gaussian'? etc.

As a result, I think the claims in the paper could be potentially misleading to users that are not aware of these restrictions, and therefore I am afraid I have to stick to reject."
Decision: Accept (poster)