A Scale-Invariant Sorting Criterion to Find a Causal Order in Additive Noise Models
Authors: Alexander Gilbert Reisach
Myriam Tami
Christof Seiler
Antoine Chambaz
Sebastian Weichwald
Venue: NeurIPS 2023 poster (2023)
Keywords: Causal Discovery
Directed Acyclic Graph
Varsortability
Additive Noise Model
Structural Causal Model
Simulation
Benchmark

Abstract: Additive Noise Models (ANMs) are a common model class for causal discovery from observational data. Due to a lack of real-world data for which an underlying ANM is known, ANMs with randomly sampled parameters are commonly used to simulate data for the evaluation of causal discovery algorithms. While some parameters may be fixed by explicit assumptions, fully specifying an ANM requires choosing all parameters. Reisach et al. (2021) show that, for many ANM parameter choices, sorting the variables by increasing variance yields an ordering close to a causal order and introduce ‘var-sortability’ to quantify this alignment. Since increasing variances may be unrealistic and cannot be exploited when data scales are arbitrary, ANM data are often rescaled to unit variance in causal discovery benchmarking.

We show that synthetic ANM data are characterized by another pattern that is scale-invariant and thus persists even after standardization: the explainable fraction of a variable’s variance, as captured by the coefficient of determination $R^2$, tends to increase along the causal order. The result is high ‘$R^2$-sortability’, meaning that sorting the variables by increasing $R^2$ yields an ordering close to a causal order. We propose a computationally efficient baseline algorithm termed ‘$R^2$-SortnRegress’ that exploits high $R^2$-sortability and that can match and exceed the performance of established causal discovery algorithms. We show analytically that sufficiently high edge weights lead to a relative decrease of the noise contributions along causal chains, resulting in increasingly deterministic relationships and high $R^2$. We characterize $R^2$-sortability on synthetic data with different simulation parameters and find high values in common settings. Our findings reveal high $R^2$-sortability as an assumption about the data generating process relevant to causal discovery and implicit in many ANM sampling schemes. It should be made explicit, as its prevalence in real-world data is an open question. For causal discovery benchmarking, we provide implementations of $R^2$-sortability, the $R^2$-SortnRegress algorithm, and ANM simulation procedures in our library CausalDisco at https://causaldisco.github.io/CausalDisco/.

A Scale-Invariant Sorting Criterion to Find a Causal Order in Additive Noise Models
Alexander G. Reisach" CNRS, MAPS Université Paris Cité F-75006 Paris, France
Myriam Tami CentraleSupélec Université Paris-Saclay F-91 190 Gif-sur- Yvette, France
Christof Seiler Department of Advanced Computing Sciences, Maastricht University, The Netherlands
Antoine Chambaz CNRS, MAP5 Université Paris Cité F-75006 Paris, France
Sebastian Weichwald Department Mathematical Sciences and Pioneer Centre for AI, University Copenhagen, Denmark
Abstract
Additive Noise Models (ANMs) are common model class for causal discovery from observational data. Due lack real-world data for which an underlying ANM is known, ANMs with randomly sampled parameters are commonly used to simulate data for the evaluation of causal discovery algorithms. While some parameters may be fixed by explicit assumptions, fully specifying an ANM requires choosing all parameters. Reisach al. (2021) show that, for many ANM parameter choices, sorting the variables by increasing variance yields an ordering close to a causal order and introduce 'var-sortability" to quantify this alignment. Since increasing variances may be unrealistic and cannot be exploited when data scales are arbitrary, ANM data are often rescaled to unit variance causal discovery benchmarking.
We show that synthetic ANM data are characterized by another pattern that is scale-invariant and thus persists even after standardization: the explainable fraction of a variable's variance, as captured by the coefficient of determination Rt2 tends increase along the causal order. The result is high "R?-sortability' meaning that sorting the variables by increasing R² yields an ordering close to a causal order. We propose a computationally efficient baseline algorithm termed 'R². SortnRegress' that exploits high It2 -rrsttbbbyyy and that can match and exceed the performance of established causal discovery algorithms. We show analytically that sufficiently high edge weights lead to a relative decrease of the noise contributions along causal chains, resulting in increasingly deterministic relationships and high synthetic data with different simulation
Rt2 We characterize R² -Foottabbyyytty parameters and find high values in common settings. Our findings reveal high It- -ortaabbbyyy as assumption about the data generating process relevant to causal discovery and implicit in many ANM sampling schemes. It should be made explicit, its prevalence in real-wol world data is an open question. For causal discovery benchmarking, we provide implementations R² sortaabiil, the R2-SortnRegress algorithm, and ANM simulation procedures in our library CausalDisco.
"Correspondence to alexander.reisach @ math.cnrs.ft
37th Conference on Neural Information Processing Systems (NcurIPS 2023).

---

Introduction
Causal Discovery in Additive Noise Models. Understanding the causal relationships between variables is a common goal across the sciences (Imbens and Rubin 2015) and may facilitate reliable prediction, the identification of goal-directed action, and counterfactual reasoning (Spirtes, Glymour, and Scheines 1993; Pearl 2009; Peters, Janzing, et al. 2017). Causal reasoning using Structural Causal Models (SCMs) consists of discovering graph encoding the causal structure between variables, learning the corresponding functions of the SCM, and using them to define and estimate causal quantities of interest. Discovering causal structure requires interventional data or assumptions on the functions and distributions to restrict the class of possible SCMs that describe a given data- generrtnng process (see Glymour al. 2019). Since collecting interventional data may be costly, infeasible, or unethical, we are interested in suitable assumptions to learn causal structure from observational data. Additive Noise Models (ANMs) encode a popular functional assumption that yields identifiability of the causal structure under various additional assumptions on the noise distributions, given the causal Markov (Kiiveri et al. 1984) and faithfulness (Spirtes. Glymour, and Scheines 1993) assumptions. In ANMs, the causal graph can be identified by testing the independence of regression residuals (see Hoyer et al. 2008; Mooij al. 2009; Peters, Mooij. et al. 2011). Another approach to use assumptions on the relative magnitudes of noise variances for identifiability (Bühlmann et al. 2014; Loh and Bühlmann 2014; Ghoshal and Honorio 2018; Chen al. 2019; Park 2020). A range of causal discovery algorithms use such assumptions to find candidate causal structures (for an overview see Heinze Deml et al. 2018; Kitson et al. 2023). Well-known structure learning algorithms include the constraint- based PC algorithm (Spirtes and Glymour 1991) and the score-based fast greedy equivalence search (FGES) algorithm (Meek 1997; Chickering 2002). Using a characterization of directed acyclic graphs as level set of differentiable function over adjacency matrices paved the way for causal discovery employing continuous optimization (Zheng et al. 2018), which has inspired numerous variants (Vowels et al. 2022).
Implicit Assumptions in Synthetic ANM Data Generation. The use of synthetic data from simulated ANMs is common for the evaluation of causal discovery algorithms because there are few real-world datasets for which an underlying ANM is known with high confidence. For such benchmark results to be indicative of how causal discovery algorithms (and estimation algorithms, see Curth et al. 2021) may perform in practice, simulated data need to plausibly mimic observations of real-world data generating processes. Reisach et al. (2021) demonstrate that data in ANM benchmarks exhibit high var-sortability, a measure quantifying the alignment between the causal order and an ordering of the variables by their variances (see Section 3 for the definition). High var-sortability implies a tendency for variances to increase along the causal order. As a consequence, sorting variables by increasing variance to obtain a causal order achieves state-of-the-art results on those benchmarks, on par with or better than those of continuous structure learning algorithms, the PC algorithm, or FGES. Var-sortability also explains why the magnitude of regression coefficients may contain more information about causal links than the corresponding p-values (Weichwald et al. 2020). However, an increase in variances along the causal order may be unrealistic because it would lead to downstream variables taking values outside plausible ranges. On top of that, var-sortability cannot be exploited on weal-world data with arbitrary data scales because variances depend the data scale and measurement unit For these reasons, following the findings of Reisach et al. (2021), Kaiser and Sipos (2022), and Seng et al. (2022), synthetic ANM data are often standardized to control the increase of variances (see, for example, Rios et al. 2021; Rolland et al. 2022; Mogensen et al. 2022; Lorch et al. 2022; Xu al. 2022; Li et al. 2023; Montagna et al. 2023), which deteriorates the performance of algorithms that rely on information in the variances. As first demonstrated for var-sortability (Reisach et al. 2021), the parameter choices necessary for sampling different ANMs constitute implicit assumptions about the resulting data-generating processes. These assumptions may result unintended patterns in simulated data, which can strongly affect causal discovery results, Post-hoc standardization addresses the problem only superficially and breaks other scalendependent assumptions on the data generating process, such as assumptions on the magnitude of noise variances. Moreover, standardization does not alter any scale-ine independent patterns that may inadvertently arise in ANM simulations. This motivates our search for scale-invariant criterion that quantifies such scale-independent patterns in synthetic ANM data. Assessing and controlling their impact on causal discovery helps make ANM simulations more realistic and causal discovery benchmark performances more informative.

---

Contribution. We show that in data from ANMs with parameters randomly sampled following common practice, not only variance, but also the fraction of variance explainable by the other variables, tends to increase along the causal order. We introduce R to quantify the alignment between a causal order and the order of increasing coefficients of determination R². In R Rertoors,,,gssss, we provide baseline algorithm exploiting R Ro-sortabyy that performs well compared to established methods on common simulation-bain based causal discovery benchmarks. In R2, we propose a sorting criterion for causal discovery that is applicable even when the data scale is arbitrary, as is often the case in practice. Unlike var-sortability, one cannot simply rescale the obtained variables to ensure desired R ssosttbbiyyyyyy in the simulated ANMs. Since real-world R² -sorabiilities are yet unknown, characterizing R² sortability and uncovering its role as a driver of causal discovery performance enables one distinguish between data with different levels of R2- -orttabbbyy when benchmarking or applying causal discovery algorithms to real-world data. To understand when high RR-sortability may be expected, parameter sampling schemes. We identify a criterion on the edge weight distribution that determines the convergence of R² to 1 along causal chains. An empirical analysis of random graphs confirms that the weight distribution strongly affects -fesstrtbill,ty, as does the average in-degree. The ANM literature tends focus on assumptions on the additive noise distributions (see for example Shimizu et al. 2006; Peters and Bühlmann 2014; Park 2020), leading to arguably arbitrary choices for the other ANM parameters not previously considered central to the model class (for example, drawing edge weights iid from some distribution). Our analysis contributes to closing the gap between theoretical results on the identifiability of ANMs and the need to make decisions on all ANM parameters in simulations. Narrowing this gap is necessary to ensure that structure learning results on synthetic data are meaningful and relevant practice.
analyze the emergence of R2 sortttiabyyyty in ANMs for different
2 Additive Noise Models
We follow standard practice in defining the model class of linear ANMs. Additionally, we detail the sampling process to generate synthetic data.
Let X = [X1, ..., Xa]T E Rd be a vector of d random variables. The causal structure between the components of X can be represented by a graph over nodes (X1..... Xa} and a set of directed edges between nodes. We encode the set of directed edges by a binary adjacency matrix B E (0. 1 where B,, = 1 if X, + Xx, that is, if X, is direct cause of X,, and Bs,: Throughout, we assume every graph to be a directed acyclic graph (DAG). Based on a causal DAG, we define an ANM. Let [ol..... ..... aa] be vector of positive iid random variables. Let PN(6) be a distribution with parameter & controlling the standard deviation. Given a draw of noise standard deviations, let N [N1, N₂) T be the vector of independent noise variables with Ne ~ (o,) For each f, let ft be measurable function such that X, = fi (Pa()))) + Ne. where Pa(X,) is the set of parents of X, in the causal graph. We assume that fi. .... fa and (6) are chosen such that all X1, .... Xa have finite second moments and zero mean. Since we can always subtract empirical means, the assumption of zero means does not come with any loss of generality; our implementations we subtract empirical means or use regression models with an intercept.
otherwise.
In this work, we consider linear ANMs and assume that fi...., fa are linear functions. For linear ANMs, we define a weight matrix W E [ddxd where Wa, is the weight of the causal link from X, to X, and W., 0 if X, is not a parent of X,. The structural causal model can be written as
X = X
2.1 Synthetic Data Generation
ANMs with randomly sampled parameters are commonly used for generating synthetic data and benchmarking causal discovery algorithms (see for example Scutari 2010; Ramsey et al. 2018: Kalainathan et al. 2020). We examine a pattern in data generated from such ANMs (Re-sortabiliiy. see Section 3.1) and design an algorithm to exploit it (e=-SornnRegssss, see Section 3.2). Here, we describe the steps to sample ANMs for synthetic data generation using the following parameters:
3

---

Parameter
Description
Number of nodes
Pg
Graph structure distribution Average node in-degree Parametrized noise distribution
(ő)
Po
Distribution of noise standard deviations Distribution of edge weights
Pw
1. Generate the Graph Structure. We generate the causal graph by sampling DAG adjacency matrix B {0, 1] for given number of nodes d. In any DAG, the variables can be permuted such that its adjacency matrix is upper triangular. We use this fact to obtain a random DAG adjacency matrix from undirected random graph models by first deleting all but the upper triangle of the adjacency matrix and then shuffling the variable order (Zheng et al. 2018). In our simulations, we the Lends--ééii (ER) (Erdös, Rényi, al. 1960) and scale-free (SF) (Barabási and Albert 1999) random graph models for Pg with an average node in-degree ?/. We denote the distribution of Erdös-Rényí random graphs with nodes and yd edges GFR (FR(d, 7d), and the distribution of scale-free graphs with the same parameters as GsF(d, 7d).
2. Define Noise Distributions. In linear ANMs, each variable X, is a linear function of its parents plus an additive noise variable N,. We choose distributional family for (d) (for example, zero-mean Gaussian or Uniform) and independently sample standard deviations ".... dd from Pa. We then set N₂ to have distribution (a.) with standard deviation d₁.
3. Draw Weight Parameters. We sample 3a.t for 8, I = 1,..., d independently from Pw and define the weight matrix W B [Batle,=1..... where denotes the element-wise product.
4. Sample From the ANM. To sample observations from the ANM with given graph (step 1), edge weights (step 2), and noise distributions (step 3), we use that W T invertible W the adjacency matrix of DAG to re-arrange Equation (1) and obtain the data-generating equation:
X= xIdd-WWT]
N.
We denote dataset of observations of X X X ...., X(m)) Budd The i-th observation
of variable in X is ( and (+) E Rd is the i-th observation vector.
Defining and Exploiting R²-Sortability
We describe a pattern in the fractions of explainable variance along the causal order and propose to measure as "R"-sortabllity.. RF-sortability is closely related var-sortability (Reisach al. 2021), which measures pattern in the variances of variables in causal graph that can be exploited for causal discovery. Similarly to var-sortability. R² -P-sortability can also be exploited to learn the causal structure of ANMs by use of simple algorithm. In contrast to var-sortability, R² If?-sortability invariant to rescaling the variables and the presented algorithm, R²2 ff=-SortRegrress recovers the causal structure equally well from raw, standardized, rescaled data.
3.1 From Var-Sortability to R²-Sortability
Var-sortability measures the agreement between ordering by variance and causal ordering. In many common ANM simulation schemes, the variances of variables tend to increase along the causal order, leading to high var-sortability. This accumulation of noise along the causal order motivates investigation of a related pattern. The variance variable X, is given as
Var(X; Var(W.,, X) + Var( (N) = Var(W., X) o²
If the variance Var (X.) increases for different X, along the causal order, but each noise variance of is sampled iid, then the fraction of cause-explained variance
Var(W.., XX
Var(W., X) o3

---

is also likely to increase for different X₁ along the causal order. Unlike the variance, we cannot calculate (nor obtain an unbiased estimate of) a variable's ' fraction of cause-explained variance without knowing its parents in the graph. Instead, we propose to use upper bound as a proxy for the fraction of cause-explained variance by computing the fraction of explainable variance of variable X, given all remaining variables, known as the coefficient of determination
R = Var (x E Xe X ((...))))))))
Var (X,)
(Glantz et al. 2001). In practice, we need to choose regression model and regress the variable onto all other variables to estimate this quantity. Here, we choose linear models Me.s TR|S| - I2 Xs + (0,Xs) Xs) for the regression of Xx onto Xs with S C {I,..., d} \ {1} and RISs, We denote the least-squares fit by Mf.s and estimate the R² coefficient of determination of X, given X3 via
(MM,s x) Var (X, - (Xs((ss)
for S [l.....d]] d3 (1)
Var (Xt)
Per definition, 1²2 is scale-invariant: it does not change when changing the scale of components of X.
To find a common definition for sortability by different ordering criteria, we introduce a family of r-sortability criteria v, that, for different functions т, assign scalar in [0, 1] to the variables X and graph G (with adjacency matrix Bg) of an ANM as follows:
ΣL E(--n)BBB En E(s-nenn)) 1
1 a < b
v. (X
where incr(a, b)
1/2 a b
(2)
a > b
and By is the i-th matrix power of the adjacency matrix Bg of graph G, and (s + f) E B3 if and only if at least one directed path from X, X; of length exists in G. In effect, ()) is the fraction of directed paths of unique length between any two nodes in G satisfying the sortability condition in the numerator. Note that other definitions of sortability that count paths differently are possible, see Appendix D for a comparison to two alternative definitions. We obtain the original var-sortability for T r(X,t) t) Var (X₂) and denote it as VVg. We obtain R²-sortability for T(X,t) R² (M ()(()))))))())> X) and denote it as V2a.
If ga is 1, then the causal order is identified by the variable order of increasing R² (decreasing if VR² = 0). A value of yæ 0.5 means that ordering the variables by R² amounts to random guess of the causal ordering. R2 -fetsoibayyy is related but not equal to var-sortability, and an ANM that may be identifiable under one criterion may not be so under the other. We present quantitative comparison across simulation settings in Appendix C, showing that when R?L-sortability is high, it tends be close to sortabibyyy (as well as sortability by cause-explained variances). The correspondence may be arbitrarily distorted by re-scaling of the variables since var-sortability is scale sensitive, whereas R -ortaabbiyy is scale-invariant In the following subsection we show that high R² if--sortbillity can be exploited similarly to high var-sortability.
3.2 Exploiting R2 -r-sortability
We introduce Tr=-SortnRegrss'' an ordering-based search algorithm (Teyssier and Koller 2005) that exploits high R2 -fessottaillity. This procedure follows Var-SortnRegress,, a baseline algorithm that obtains state-of-the-art performance in causal discovery benchmarks with high var-sortability (Reisach et al. 2021). Var-SortmRegress is computationally similar to the Bottom- Up and Top- Top-Down algorithms by Ghoshal and Honorio (2018) and Chen al. (2019), which rely on an equal noise variance assumption. However, unlike the assumptions underlying these algorithms, the R² criterion ceale-nnvrriant such that H² Rorrrreeessssss is unaffected by changes in data scale.
In RP-SortnRegress ( (Algorithm 1), the coefficient of determination is estimated for every variable to obtain a candidate causal order. Then, regression of each node onto its predecessors in that order is performed to obtain the candidate causal graph. To encourage sparsity, one can use penalty function
Reisach et al. (2021) refer to Rergressssssress as SortnRegress. to H² Sortmnegress...
We add the Var to emphasize the contrast
5

---

X(f). In our implementation, we use an L1 penalty with the regularization parameter chosen via the Bayesian Information Criterion (Schwarz 1978).
As before, we denote by My. 7. parameterized model X, with covariates Xg and parameters 0. For an example, consider the linear model Mf,s: Xs (0, Xs) with 0, -RISI, To simplify notation, define the mean squared error (MSE) of model MSE(X,, MY.s) n - Mf,s(X')) ))) )²
Algorithm 1: R SortoRReeessss
Data: X E Baxd
Input: À /* Penalty function, for example X(0) = =|
Result: Weight matrix estimate W E Baxd
Candidate causal order
f 0 R
** Estimate R using d regressions (with intercept; cf, Section 2)
for = I,..., d do
arg ming MSE (Xt, (i..i...
1²² (MM (M
Find permutation p that sorts x ascending
/s Estimate weights using d 1 regressions (with intercept)
W- gxxd
for i= 2.......
t= [j: ((j =
S = {j: p(j) < {
W... ming MSE (X.,M..) + A(0)
return W
4 Utilizing -2-Sortabylity for Causal Discovery
We compare Re=SortnRegresss to a var-sortability baseline and established structure learning algo- rithms in the literature across simulated data with varying (f-ssortaillities, and evaluate the algorithms on the real-world protein signalling dataset by Sachs et al. (2005).
4.1 Experiments on Simulated Data
We evaluate and compare Re=SorstResress on observations obtained from ANMs simulated with the following parameters:
Pg (GER (20, 40), (20
(Erdōs-Rényí and scale-free graphs, 20, 7 (Gaussian noise distributions)
(Ø) (0, x²)
P. Unif (0.
(noise standard deviations)
Pw
(0.5,1))
(edge weights)
In this Gaussian setting, the graph is not identifiable unless additional specific assumptions on the choice of edge and noise parameters are employed. We sample 500 ANMS (graph, noise standard deviations, and edge weights) using ER graphs, and another 500 using SF graphs. From each ANM, we generate 1000 observations. Data are standardized to mimic real-world settings with unknown data scales. We compare R²- Re-SorsnRegress to the PC algorithm (Spirtes and Glymour 1991) and FGES (Meek 1997; Chickering 2002). Representative for scale-sensitive algorithms, we include Var-SortnRegress, introduced by Reisach al. (2021) to exploit var-sortability, as well as their RandomRegress baseline (we refer to their article for comparisons to further algorithms). For details on the implementations used, see Appendix A. The completed partially directed acyelic graphs (CPDAGs) returned by PC and FGES are favorably transformed into DAGs by assuming the correct direction for any undirected edge there is a corresponding directed edge in true graph, and breaking any remaining cycles. We evaluate the structural intervention distance (SID) (Peters and Bühlmann 2015) between the recovered and the ground-truth DAGs. The SID counts the number of incorrectly estimated interventional distributions between node pairs when using parent- adjunttadjustment

---

in the recovered graph compared to the true graph (lower values are better). It can be understood as measure of disagreement between a causal order induced by the recovered graph and causal order induced by the true graph.
Structural Intervention Distance
FC
J Sorinflegrew
Structural Intervention Distance
200
PC
Re SortnRearce
FIES
175
FGES
Rardorffagres drr))
000
150
rats)
125
1.3
100
76
60
25
0.4
0.6
0.5
0.5
0.50
0.70
0.70
0,50
0.50
0.90
630
soortaability
(a) Erdōs Rényi graphs
(b) Scale-free ffee-feee graphs
Figure 1: Performance comparison on standardized synthetic data in terms of SID (lower is better) using moving window averages. Rº- R7e-Sortneegress (blue) performs well if RR-sortability is high, achieving results competitive with established methods. For reference, we show the performance achieved by Var-SortnRegress on raw data (gray dashed line), which worsens to that of Random- Regress (gray) when the raw data scale is unavailable, as is the case here after standardization.
The results shown in Figure 1. To analyze trends, windows of width of 0.1 centered around 0.05, SID from window window, and the shaded areas indicate the 95% confidence interval of the mean. The means of non-empty windows start at 0.25 for ER graphs and at 0.65 for SF graphs. For both graph types, there are many instances with R -orttabbbly well above 0.5. This is notable since we simulate standard Gaussian non-equal variance setting with parameters chosen in line with common practice in the literature (see for example the default settings of the simulation procedures provided by Squires 2018; Kalainathan et al. 2020; Zhang al. 2021). We present an analogous simulation with lower absolute edge weights in (-0.5, 0.1) (0.1, 1, 0.5) Appendix A. show results in terms the structural Hamming distance (SHD) in Appendix A.2, and provide a visualization of the distribution of -ortaabiily for higher and lower edge weights in Appendix A.3.
window averaging with overlapping 0.95. The lines indicate the change in mean
R2 Re-Sresnneeress successfully exploits high Il"-sorabiliity for causal discovery. On ER graphs (Figure la), Re=SSortRegress outperforms the RandomRegress baseline when 1?*-sortability is greater than 0.5. For R -o-strtabilisies close to 1, RRe-oortnegrrsss outperforms PC and matches or exceeds the performance of FGES. On SF graphs (Figure 1b), -F-sortabiliitis are generally higher, and R² Re-SoosnRegress outperforms PC and FGES for all but the lowest values, even coming close to the performance of Var- Var-Sosrmeegress on raw data for very high R2 sotaabbbbslssisss On raw synthetic data, Var-SortnRegress successfully exploits the high var- cartsorabyylity typical for many ANM simulation settings. Its performance improves for high R Re-sortability showing the link between R2 -sortabiility and var-sortability. However, this performance relies on knowledge of the true data scale to exploit var-sortability. The same is true for other scale-sensitive scores such as the one used by NOTEARS (Zheng et al. 2018) and many other algorithms that build upon it. Without information in the variances (as the case after standardization), Var-Sorttneerress is reduced to RandomRegress, which consistently performs poorly across all settings. FGES outperforms PC, which may be explained by FGES benefitting from optimizing a likelihood specific to the setting, thus using information that is not available PC. Further causal discovery results for settings with different noise and noise standard deviation distributions are shown in Appendix A.4.
Real- World Data
We analyze the observational part of the Sachs et al. (2005) protein signalling dataset to provide an example of what 12- -o"sortabilities and corresponding causal discovery performances may be expected in the real world. The dataset consists of 853 observations of 11 variables, connected by 17 edges in the ground-trutt consensus graph. We find the RR-sortability to be 0.82, which above the neutral value of 0.5, but not as high as the values for which R R2-SortnRegress outperforms the other methods
7

---

in the simulation shown in Section 4.1. We run the same algorithms as before on this data and evaluate them in terms of SID and SHD to the ground-truth consensus graph. To gain insight into the robustness of the results, we resample the data with replacement and obtain 30 datasets with a mean R² sortsaibbyy of 0.79 and mean var-sortability of 0.65. The results (mean, [min, max]) are shown in the table below. We find that none of the algorithms perform particularly well (for reference: the empty graph has an SID of 53 and an SHD of 17), Although the other algorithms perform better than Ragres pronounced than in our simulation (Section 4.1).
average with PC performing best, the differences between them are less
Algorithm
SID
SHD
H2 Sortoneeessssss
48.13
[43. 51]
14.43
[12, 17]
PC
43.70
[37.49]
11.40
9, 16]
FGES
48.43
[45,
12.77
[11, 15]
RandomRegress
52.57
[42, 591
16.20
[14, 19]
Var-Sertn
45.13
[43.47]
15.53
[13, 18]
These results showcase the difficulty of translating simulation performances to weal-world settings. For benchmarks to be representative of what to expect in real-world data, benchmarking studies should differentiate between different settings known to affect causal discovery performance. This includes R - and we therefore recommend that benchmarks differentiate between data with different levels of R -orttabbiyy and report the performance of Re-SortnRegress baseline.
5 Emergence and Prevalence of R -R-Sttbbbyyyyy in Simulations
The Re-sortability of an ANM depends the relative magnitudes of the R2 coefficients of determi- nation of each variable given all others, and the R coefficients in turn depend on the graph structure, noise variances, and edge weights. Although one can make the additional assumption of an equal fraction of cause-explained variance for all non-root nodes (see, for example, Sections 5,1 of Squires al. 2022; Agrawal et al. 2023), this does not guarantee R -o-soraability close to 0.5, and the R² may still carry information about the causal order as we show in Appendix E. Furthermore, one cannot isolate the effect of individual parameter choices on RR-sortability, nor easily obtain the expected -oo-sotabyllity when sampling ANMs given some parameter distributions, because the R² values are determined by complex interplay between the different parameters. By analyzing the simpler case of causal chains, we show that the weight distribution plays crucial role in the emergence of Ro-sortability. We also empirically assess R² sortaiibyyyy for different ANM simulations to elucidate the connection between R² -F"sootabilityy and ANM parameter choices.
5.1 The Edge Weight Distribution as a Driver of R -Sottabiiity in Causal Chains
Our goal is to describe the emergence of high R along causal chain. We begin by analyzing the node variances. In a chain from X, to Xp of given length p > 0, the structure of & is fixed. Thus, sampling parameters consists of drawing edge weights Wa,1 .... Wp-1.p Pw and noise standard deviations 0o ..., Pp Po (cf. Section 2. Importantly, these parameters are commonly drawn iid, and we adopt this assumption here. As we show in Appendix B. the variance of the terminal node p in the resulting chain ANM is lower bounded by the following function of the random parameters
Xp
log [Wis+al-
Y=U
We assume that the distribution of noise standard deviations P. has bounded positive support, as is commonly the case when sampling ANM parameters (see Squires 2018; Kalainathan et al. 2020; Zhang et al. 2021, as well as our simulations). We introduce the following sufficient condition for a weight distribution Pw to result in diverging node variances along causal chains of increasing length
0 < E [log []] < +oo, with V Pw.
If Equation (4) holds, the lower bound on the node variances given in Equation (3) diverges almost surely by the strong law of large numbers
log W Witt1l a.s.
+00,
p-000
8

---

since the j,1+1 are sampled iid from Pw. This implies that the variance of a terminal node in chain of length p also diverges almost surely to infinity as p goes to infinity. Under the same condition, the R2 of the terminal node given all other nodes converges almost surely to 1 as p goes to infinity. We can show this using the fraction of cause-explained variance as lower bound on the R² of variable Xp and the divergence of the variances to infinity:
Var (Xp E Xp p ((oo....N((pt 2 21-
Var (Xp 13 [Xx 1x-11)
R²2 =
Var (Xp)
Var (Xp)
as.
1.
Var
p-200
For Pw with finite E [log [V] with V Pw. it is thus sufficient that E log [V] > 0 for the terminal node's R² in a p-chain to converge to almost surely as p - +00, under iid sampling of edge weights and noise standard deviations. The fact that Equation (4) only depends on Pw highlights the importance the weight distribution as a driver of R².
Increasingly Deterministic Relationships Along the Causal Order. Our argument for the almost sure convergence of R along causal chains relies on a) the divergence of the total variance of the terminal node, given Equation (4) holds, and b) independently drawn noise variances o² from bounded distribution. Together, these two observations imply that the relative contributions of the noise variances to the total variances tend to decrease along a causal chain. If data are standardized to avoid increasing variances (and high var-sortability), the noise variances tend vanish² along the causal chain, resulting in increasingly deterministic relationships that effectively no longer follow an ANM
5.2 Impact of the Weight Distribution R it--Sortayylity in Random Graphs
Understanding the drivers of R² -Posttabyyyyty in random graphs is necessary to control the level of R²- sortability in data simulated from ANMs with randomly sampled parameters. Given the importance of the weight distribution in our analysis of -r=ssotabillity in causal chains, we empirically evaluate the impact of the weight distribution -ossoraability in random graphs. We provide additional experiments showing how the weight distribution and average in-degree impact Re-sortabilityy in Appendix B.2.
We analyze Rº -?2soortbilityy for different graph sizes. We use ER and SF graphs, that is Pg E (G 2d), Gsr(d, ] with d nodes and 2, and simulate graphs with d 2 5 such that 2d does not exceed the maximum number possible edges. Let V be a random variable with V ~ Pw. We test different edge weight distributions Pw Unif((-a, -0.1) ) U (0.1, a)) with such that the E In [V] are approximately evenly spaced in [-1, 2.5). We otherwise use the same settings in Section 4.1. For each graph model, weight distribution, and number of nodes, independently simulate 50 ANMs (graph, noise standard deviations, and edge weights).
0.1 chosen
Weight
T porttabilltty
Pi-sortabiitty
P1]-25
2,9. []]=20
08-
1.11 oep 80h K0--111-5
0.-
0.7
=/0.1600) Bll [11/P11110
+49.1400) Rih
0.6
0.6
6,9 [V]]-o0
=(0.1.130) Els [7]]-5
0.8
+/9.1.273)
R0x (00
10
30
6
37
Nodes
Nodes
(a) Erdōs
Erdōs-Rééyi Rényi graphs
(b) Scale-free graphs
Figure 2: -ortaabbyy at different graph sizes for range of weight distributions with different
[V] with V Pw.
The results are shown in Figure 2. We observe that R² -orttbbbbyy is consistently higher for edge weight distributions Pw with higher E [In [V]]. For each weight distribution, observe slowing
"The noise contribution standardized variable X (wy-1.) X + N,)/ Var (X;) in causal chain, given Var(N,/ / Var (X;)), vanishes for Var (X,) +00.

---

increase in R sorttabbiy as graph size increases. Compared to ER graphs (Figure 2a), R -ortabiilyyy in SF graphs tends increase faster in graph size and is consistently higher on average (Figure 2b). In all cases, mean R2 sortabiiiyy exceeds 0.5, and weight distributions with E [In [V] > 0 reach levels of -orttabbbiyy for which R Rerteeenee achieves competitive causal discovery performance in our simulation shown in Section 4. Settings with E [log [V]] < 0 may result in lower R -ossottabiliitis, but also make detecting statistical dependencies between variables connected by long paths difficult as the weight product goes to O. In summary, we find that the weight distribution strongly affects R2 -ortsaibayyl and can be used to steer -i--sortaillity in simulations.
How the weight distribution affects -sstttbbbyy also depends on the density and topology of the graph, giving rise to systematic differences between ER and SF graphs. In Appendix B.2 we show the R² sortabiiiy for graphs of 50 nodes with different E [In [V]] and average in in-degree In Erdōs-Rényi graphs, we find that -ossortabilities are highest when neither of the parameters take extreme values, as is often the case in simulations. In scale-free graphs we observe -sssortabilities close to across all parameter settings. These results show that the graph connectivity interacts with the weight distribution and both can affect R -orsoriabylyity and thereby the difficulty of the causal discovery task.
Discussion and Conclusion
We introduce R² as simple ocale-invariant sorting criterion that helps find the causal order in linear ANMs. Exploiting R² -orstabbbllity can give competitive structure learning results if -o-sortability is sufficiently high, which we find to be the case for many ANMs with iid sampled parameters. This shows that the parameter choices necessary for ANM simulations strongly affect causal discovery beyond the well-understood distinction between identifiable and non-identifiable settings due to properties of the noise distribution (such as equal noise variance or non- non-Gaussianity, High R². sortability particularly problematic because high R² are closely related to high variances. When data are standardized to remove information in the variances, the relative noise contributions can become very small, resulting in a nearly deterministic model. Alternative ANM sampling procedures and the R² -o-sortabilliy of nonlinear ANMs and structural equation models without additive noise are thus of interest for future research. Additionally, determining what sortaabbilsss may be expected in real-world data is an open problem. A complete theoretical characterization of the conditions sufficient and/or necessary for extreme R could help decide when to assume and exploit it in practice. Exploiting -tassorbbiylity could provide way of using domain knowledge about parameters such as the weight distribution and the connectivity of the graph for causal discovery. On synthetic data, the prevalence of high R -f?-sttabilityy raises the question of how well other causal discovery algorithms perform relative to the H² If"-sortbillity of the benchmark, and whether they may also be affected by R² -f*-sortaillity. For example, the same causal discovery performance may be less impressive if the s?e-sortability on a dataset is 0.95 rather than 0.5. In this context, the R² oe-SortmRegress algorithm offers strong baseline that is indicative of the difficulty of the causal discovery task. Our empirical results on the impact of the weight distribution and graph connectivity provide insight how steer R² -ortabbbiyyy within existing linear ANM simulations. Simulating data with different so=sortabilittes can help provide context to structure learning performance and clarifies the impact of simulation choices, which is necessary evaluate the real-world applicability of causal discovery algorithms.
Acknowledgements
We thank Brice Hannebicque for helpful discussions earlier draft of this paper. AGR received funding from the European Union's 's Horizon 2020 research and innovation programme under the Marie Sktodowska-Curd grant agreement No 945332
10

---

References
[1] Raj Agrawal, Chandler Squires, Neha Prasad, and Caroline Uhler. The DeCAMFounder: nonlinear causal discovery in the presence of hidden variables. In: Journal of the Royal Statistical Society Series B: Statistical Methodology (July 2023), qkad071 (cit. on pp. 8, 22).
[2] Albert-L.ászle Barabási and Réka Albert. Emergence of Scaling in Random Networks. In:
Science 286.5439 (1999), pp. 509-512 (cit. on p. 4)..
[3] Peter Bühlmann, Jonas Peters, and Jan Emnest. CAM: Causal Additive Models, High- Dimensional Order Search and Penalized Regression. In; The Annals of Statistics 42.6 (2014),
pp. 2526-2556 (cit. on p. 2).
[4] Wenyu Chen, Mathias Drton, and Y Samuel Wang. On causal discovery with an equal-variance assumption. In: Biometrika 106.4 (2019), pp. 973-980 (cit. on pp. 2, 5, 14).
[5] David Maxwell Chickering. Optimal Structure Identification With Greedy Search. In: Journal of Machine Leurning Reseurch 3 (2002), pp. 507-554 (cit. pp. 2, 6).
[6] Gabor Csardi, Tamas Nepusz, et al. The igraph software package for complex network research. In: InterJournal Complex Systems 1695.5 (2006), pp. 1-9 (cit. p. 14).
[7] Alicia Curth, David Svensson, Jim Weatherall, and Mihaela van der Schaar. Really Doing Great Estimating CATE? A Critical Look at ML Benchmarking Practices in Treatment Effect Estimation. In: Proceedings of the Neural Information Processing Systems Track on Datasets and Benchmarks. Vol. 1. 2021 (cit. on p. 2).
[8] Paul Erdōs, Alfréd Rényi, et al. On the Evolution of Random Graphs. In: Publ. Math. Inst. Hung. Acad. Sci 5.1 (1960), pp. 17-60 (cit. on p. 4).
[9] Asish Ghoshal and Jean Honorio. Learning linear structural equation models in polynomial time and sample complexity. In: International Conference on Artificial Intelligence and Statistics. PMLR, 2018, pp. 1466-1475 (cit. on pp. 5)..
[10] Stanton A. Glantz, Bryan K. Slinker, and Torsten B. Neilands. Primer of Applied Regression & Analysis of Variance. Vol. 654. McGraw Hill Inc., 2001 (cit. on p. 5).
[11] Clark Glymour, Kun Zhang, and Peter Spirtes. Review of Causal Discovery Methods Based on Graphical Models. In: Frontiers in Genetics 10 (2019), 524 (cit. on p. 2).
[12] Charles R. Harris, K. Jarrod Millman, Stéfan J. van der Walt, Ralf Gommers, Pauli Virtanen, David Cournapeau, Eric Wieser, Julian Taylor, Sebastian Berg, Nathaniel J. Smith, Robert Kern, Matti Picus. Stephan Hoyer, Marten H. van Kerkwijk, Matthew Brett, Allan Haldane, Jaime Fernández del Río, Mark Wiebe, Pearu Peterson, Pierre GGraadd-archant Kevin Sheppard, Tyler Reddy, Warren Weckesser, Hameer Abbasi, Christoph Gohlke, and Travis E. Oliphant. Array programming with NumPy. In: Nature 585.7825 (Sept. 2020), pp. 357-362 (cit. on p. 14).
[13] Christina Heinze-Deml, Marloes H. Maathuis, and Nicolai Meinshausen. Causal Structure Learning. In: Annual Review of Statistics and Its Application 5 (2018), pp. 371-391 (cit. on
p. 2).
[14] Patrik O. Hoyer, Dominik Janzing. Joris M. Mooij, Jonas Peters, Bernhard Schölkopf. et al. Nonlinear causal discovery with additive noise models. In: Advances in Neural Information Processing Systems. Vol. 21. Citeseer. 2008, pp. 689-696 (cit. on p. 2).
[15] Guido W. Imbens and Donald B. Rubin. Causal Inference for Statistics, Social, and Biomedical Sciences. Cambridge University Press, 2015 (cit. on p. 2).
[16] Marcus Kaiser and Maksim Sipos. Unsuitability of NOTEARS for Causal Graph Discov- ery when Dealing with Dimensional Quantities. In: Neural Processing Letters 54.3 (2022),
pp. 1587-1595 (cit. on p. 2).
[17] Diviyan Kalainathan, Olivier Goudet, and Ritik Dutta. Causal Discovery Toolbox: Uncover- ing causal relationships in Python. In: Journal of Machine Learning Research 21.1 (2020),
pp. 1406-1410 (cit. on pp. 3, 7, 8).
[18] Maurice G. Kendall. A new measure of rank correlation, In: Biometrika 30.1/2 (1938), pp. 81- 93 (cit. on p. 21).
[19] Harri Kiiveri, Terry P. Speed, and John B. Carlin. Recursive causal models. In: Journal of the Australian Mathematical Society 36.1 (1984), pp. 30-52 (cit. on p. 2).

---

[20] Neville K. Kitson, Anthony C. Constantinou, Zhigao Guo, Yang Liu, and Kiattikun Chobtham.
A survey of Bayesian Network structure learning. In: Artificial Intelligence Review (2023).
pp. 1-94 (cit. on p. 2).
[21] Chunlin Li, Xiaotong Shen, and Wei Pan. Nonlinear Causal Discovery with Confounders, In: Journal of the American Statistical Association (2023) (cit. on p. 2).
[22] Po-Ling Loh and Peter Bühlmann. High- figh-dimensional learning of linear causal networks via inverse covariance estimation. In: Journal of Machine Learning Research 15.1 (2014), pp. 3065-3105 (cit. on p. 2).
[23] Lars Lorch, Scott Sussex, Jonas Rothfuss, Andreas Krause, and Bernhard Schölkopf. Amor- tized Inference for Causal Structure Learning. In: NeurIPS 2022 Workshop on Causality for Real-world Impact. 2022 (cit. on p. 2).
[24]
Christopher Meek. Graphical Models: Selecting causal and statistical models. PhD thesis, Carnegie Mellon University, 1997 (cit. on pp. 2, 6).
[25] Phillip B. Mogensen, Nikolaj Thams, and Jonas Peters. Invariant Ancestry Search. In: Interna- tional Conference Machine Learning. PMLR, 2022, pp. 15832-15857 (cit. p. 2).
[26] Francesco Montagna, Atalanti A. Mastakouri, Elias Eulig. Nicoletta Noceti, Lorenzo Rosasco, Dominik Janzing, Bryon Aragam, and Francesco Locatello. Assumption violations in causal discovery and the robustness of score matching. In: arXiv preprint arXiv:2310.13387 (2023) (cit. p. 2).
[27] Joris M. Mooij, Dominik Janzing, Jonas Peters, and Bernhard Schölkopf. Regression by dependence minimization and its application to causal inference in additive noise models. In: Proceedings of the 26th annual international conference on machine learning. 2009, pp. 745- 752 (cit. on p. 2).
[28]
Gunwoong Park. Identifiability of Additive Noise Models Using Conditional Variances. In: Journal of Machine Learning Research 21.75 (2020), pp. 1-34 (cit. on pp. 2, 3).
[29]
Judea Pearl. Causality. Cambridge University Press, 2009 (cit. on p. 2).
[30] Fabian Pedregosa, Gaël Varoquaux, Alexandre Gramfort, Vincent Michel, Bertrand Thirion, Olivier Grisel, Mathieu Blondel, Peter Prettenhofer, Ron Weiss, Vincent Dubourg, et al. Scikit- learn: Machine learning in Python. In: Journal of Machine Learning Research 12 (2011), pp. 2825-2830 (cit. on p. 14).
[31] Jonas Peters and Peter Bühlmann. Identifiability of Gaussian structural equation models with equal error variances. In: Biometrika 101.1 (2014), pp. 219-228 (cit. on p. 3).
[32] Jonas Peters and Peter Bühlmann. Structural Intervention Distance for Evaluating Causal Graphs. In: Neural computation 27.3 (2015), pp. 771-799 (cit. on pp. 6, 14),
[33] Jonas Peters, Dominik Janzing, and Bernhard Schölkopf. Elements of Causal Inference: Foundations and Learning Algorithms. The MIT Press, 2017 (cit. on p. 2).
[34] Jonas Peters, Joris M. Mooij, Dominik Janzing, and Bernhard Schölkopf. Identifiability of Causal Graphs using Functional Models. In: Proceedings of the Seeeny- Conference on Uncertainty Artificial Intelligence. AUAI Press, 2011, pp. 589-598 (cit. on p. 2).
[35] Joseph D. Ramsey, Kun Zhang, Madelyn Glymour, Ruben S. Romero, Biwei Huang, Imme Ebert-Uphoff, Savini Samarasinghe, Elizabeth A. Barnes, and Clark Glymour. TETRAD- AA Toolbox for Causal Discovery. In: 8th International Workshop on Climate Informatics. 2018 (cit. on pp. 3, 14).
[36] Alexander G. Reisach, Christof Seiler, and Sebastian Weichwald. Beware of the Simulated DAG! Causal Discovery Benchmarks May Be Easy To Game In: Advances Neural Infor- mation Processing Systems. Vol. 34. 2021, pp. 27772-27784 (cit. on pp. 1, 2, 4-6, 14, 19,
20).
[37] Felix L. Rios, Giusi Moffa, and Jack Kuipers. Benchpress: A Scalable and Versatile Workflow for Benchmarking Structure Learning Algorithms. In: arXiv preprint urXiv:2107.03388 (2021)
(cit. on p. 2).
[38] Paul Rolland, Volkan Cevher, Matthäus Kleindessner, Chris Russell, Dominik Janzing, Bern- hard Schölkopf, and Francesco Locatello. Score Matching Enables Causal Discovery of Nonlinear Additive Noise Models. In: International Conference on Muchine Learning. PMLR, 2022. pp. 18741-18751 (cit. on 2).
12

---

[39] Karen Sachs, Omar Perez, Dana Pe' er, Douglas A. Lauffenburger, and Garry P. Nolan. Causal signain-ignaling networks derived from multiparameter single-cell data. In: Science 308.5721 (2005), pp. 523-529 (cit. pp. 6, 7).
[40] Gideon Schwarz. Estimating the Dimension of a Model. In: The Annals of Statistics 6.2 (1978), pp. 461-464 (cit. on p. 6).
[41] Marco Scutari. Learning Bayesian Networks with the bnlearn Package. In: Journal of Statistical Software 35.3 (2010) (cit. on p)
[42] Skipper Seabold and Josef Perktold. Statsmodels: Econometric and Statistical Modeling with Python. In: Proceedings of the 9th Python in Science Conference. Vol. 57. 61. 2010, pp. 10-
25080 (cit. on p. 14).
[43] Jonas Seng, Matej Zečević, Devendra Singh Dhami, and Kristian Kersting, Tearing Apart NOTEARS: Controlling the Graph Prediction via Variance Manipulation. In: arXiv preprint arXiv:2206.07111 (2022) (cit. on p. 2).
[44] Shohei Shimizu, Aapo Hyvärinen, Patrik O. Hoyer, and Yutaka Kano. Finding a causal ordering via independent component analysis. In: Computational Statistics & Data Analysis 50.11 (2006), pp. 3278-3293 (cit. on p. 3).
[45] Peter Spirtes and Clark Glymour. An Algorithm for Fast Recovery of Sparse Causal Graphs. In: Social Science Computer Review 9.1 (1991), pp. 62-72 (cit. on pp. 2, 6).
[46] Peter Spirtes, Clark Glymour, and Richard Scheines. Causation, Prediction, and Search. 1993.
In: Lecture Notes in Statistics 81 (1993) (cit. on p. 2).
[47] Chandler Squires. causaldag: creation, manipulation, and learning of causal models. 2018
(cit. on pp. 7, 8).
[48] Chandler Squires, Annie Yun, Eshaan Nichani, Raj Agrawal, and Caroline Uhler. Causal Structure Discovery between Clusters of Nodes Induced by Latent Factors. In: Proceedings Causal Learning and Reasoning. Vol. 177. PMLR, Nov. 2022,
of the First Conference
pp. 669-687 (cit. on pp. 8, 22, 23).
[49] Mare Teyssier and Daphne Koller. Ordering- Baasei
Based Search: A Simple and Effective Algo- rithm for Learning Bayesian Networks. In: Proceedings of the Twenty-Fin Fiirtt Conference on Uncertainty in Artificial Intelligence. AUAI Press, 2005, pp. 584-590 (cit. on p. 5)..
[50] Matthew J. Vowels, Necati Cihan Camgoz, and Richard Bowden. D'Ya Like DAGs? A Survey on Structure Learning and Causal Discovery. In: ACM Computing Surveys 55.4 (2022) (cit. on p. 2).
[51] Sebastian Weichwald, Martin E. Jakobsen, Phillip B. Mogensen, Lasse Petersen, Nikolaj Thams, and Gherardo Varando. Causal structure learning from time series: Large regression coefficients may predict causal links better practice than small p-values. In: Proceedings of the NeurlPS 2019 Competition and Demonstration Track, Proceedings of Machine Learning Research. Vol. 123. PMLR. 2020, pp. 27-36 (cit. on p. 2).
[52] Sascha Xu, Osman A. Mian, Alexander Marx, and Jilles Vreeken. Inferring Cause and Effect in the Presence of Heteroscedastic Noise. In: International Conference on Machine Learning. PMLR, 2022, pp. 24615 2-1630 (cit. on p. 2).
[53] Keli Zhang, Shengyu Zhu, Marcus Kalander, Ignavier Ng, Junjian Ye, Zhitang Chen, and Lujia Pan. gCastle: A Python Toolbox for Causal Discovery. In: arXiv preprint arXiv:2111. 1515555 (2021) (cit. on pp. 7,8, 8).
[54] Xun Zheng, Bryon Aragam, Pradeep K Ravikumar, and Eric P Xing. DAGs with NO TEARS: Continuous Optimization for Structure Learning. In: Advances Neural Information Process- ing Systems. Vol. 32. 2018, pp. 9472-9488 -973-9433--- (cit. on pp. 2, 4, 7).
13

---

Appendix
A Causal Discovery Experiments
We provide an open-source implementation of - de-SorsnRegress, and an ANM sam- pling procedure in our library CausalDisco. In our experiments we make extensive use of the NumPy (Harris et al. 2020), seikit-learn (Pedregosa 2010) Python packages, as well as the Python interface to the igraph (Csardi, Nepusz, et al. 2006) package. For the PC and FGES algorithms we use the implementation by Ramsey et al. (2018). The CPDAGs returned by PC and FGES are favorably transformed into DAGs by assuming the correct direction for any undirected edge if there is corresponding directed edge in the true graph. To break any remaining cycles, we iteratively pick one of the shortest ones and remove one of the edges that break the most cycles until no cycles remain. For var-sortability and Var-SortnRegress we use the implementations provided by Reisach al. (2021). In the experiments shown in Appendix A.4 additionally show the Top-Down algorithm by Chen et al. (2019) using their implementation. To evaluate structure learning performance in SID (Peters and Bühlmann 2015) run the authors' implementation using version 3.6.3 of the R programming language.
2011), and statsmodels (Seabold and Perktold
A.1 Causal Discovery Performance and R² -orttaibyyyyty Low Absolute Edge Weights
Given the impact of high absolute edge weights on R² as shown in Section 5, it may seem tempting to simulate ANMs with low absolute edge weights to avoid high levels of - To assess the potential of this idea, simulate setting with low absolute edge weights. The simulation in Sec- tion 4.1 draws from Pw = Unif((-2, (05) 2)), resulting in E [log [V] 0.16 with ~ Pw. Here, we draw from Pw Unif((-0.5, 0.1) U (0.1, 0.5)), resulting in E [log [V]] - 1.29 with
V ~ Pw, while keeping the remaining settings the same as in Section 4.1.
Structural Intervention Distance
300-
SartnRogrrss
Structural Intervention Distance
PC
Sartrileara
FGES
IGES
Rander Regess
2014
Rayfow
2004
(ce van data)
Vie SommSmrreee
175-
[en raw data)
Vo- Vy--Sonillggee
200-
150
125
190
100-
1004
75
50
65
65
0'
0.6
0.0
0.65
0.50
0.55
So-sortabliity
Ro-ttraaiility
(a) ER graphs
(b) SF graphs
Figure 3: Performance comparison on simulated data with low absolute edge weights in terms of SID (lower is better). oe-SorsnRegress (blue) performs well R² -o-sotaabyl is high, achieving results competitive with established methods. For reference. show the performance achieved by Rar-SootnRegres on raw data (gray dashed line), which worsens to that RandomRegress (gray) when the raw data scale is unavailable, is the case here after standardization.
The results are shown in Figure Compared to the setting with higher absolute edge weights shown in Section 4. we observe lower values of R?²-sortability overall. On ER graphs, as shown in Figure 3a, Rt-SSortReegres performs worse than RandomRegress for R²- If2-sottaillitie below 0.5 and outperforms it for RR- -o-sttaabllsstes greater than 0.5, as expected from the definition of R² -ortobibilyy which yields 0.5 to be the neutral value, For R² sortaabiillsies closer to 1, R -SrrtnReggesss outperforms PC and narrows the gap to FGES. Var-SortnRegress performs similar to PC and worse than FGES in this setting, likely because settings with higher edge weights (Reisach et al. 2021). On SF graphs, as shown in Figure 3b, R² -o-soraabilities are generally above the neutral value of 0.5 despite the lower absolute edge weights, and R22-SornResress performs well, even outperforming all other algorithms many instances. The
lower edge weights mitigates the variance explosion that characterizes
14

---

excellent absolute and relative performance of R² -SrrtnReggrsss and the high values of R -ortaabilyyy indicate that the SF graph structure may be particularly favorable for an exploitable R² pattern arise. Choosing lower absolute edge weights results in lower -ottabiabil and degrades the comparative performance of Te-SoosnRegress on ER graphs and, lesser extent, also on graphs. However, Ra Rertersnessssss still performs excellently on SF graphs and can achieve good results on ER graphs with high -f-sortabillit, which do occur despite the low absolute edge weights. Overall, extreme R2 -sortabbilities are less frequent for low absolute edge weights, but they are not impossible and can impact causal discovery performance. therefore recommend reporting R2 -sortabbiyy regardless of the edge weight range.
A.2 Evaluation of Causal Discovery Results Using Structural Hamming Distance
Here we show causal discovery performances in terms of the structural Hamming distance (SHD). Note that SHD, although frequently used as a performance measure in causal discovery, cannot be directly interpreted terms of the suitability the discovered graph for effect estimation, except if it is zero. For example, excess edges to causal predecessors to causal successors in the true graph despite their different implications for effect estimation. For the same reason, SHD results very sensitive to the density of the true graph. As result, the absolute performance of regression-based methods such as fR-SSortnRegresx, Var-SortnRegress, and RandomRegress strongly depends the the sparsity penalty used. For additional context, recall that the empty graph has constant SHD) the product of average in-degree and number of nodes ydd, which is equal to 40 for the settings shown here.
treated the same excess edges
S-SSatnilaagee
Structural Hamming Distance
Strictural Hamming Distance
PCFS Rarde
Regean
80
Van SoniPogris
VaeSuni Rogess (en raw data)
033
0.5
ULT
0.9
0.50
0.70
0,70
0.5
010
RC
sortability
RP sartability
ER graphs
SF graphs
(a) Higher edge weights: setting from Section 4.1 with Pw Unif([-1, -05) (0.5, 1))
p':Sootafiingree
It' Sertnitegree
Structural Hamming Distance
Structural Hamming Distance
PC
20-
FCES Dasteni
Rendorrllggees
lagese
Van-Sorrliggee
10
30-
10
0.5
0.6
0.6
0.7
0.5
0.55
0.65
0.70
0.75
0.50
0.56
0.00
0.06
R sosttrabillity
R2
Ro-ttrrabilityy
ER graphs
SF graphs
(b) Lower edge weights: setting from Appendix A.1 with Pw = (.0)-.00
Figure 4: Performance comparison simulated data with high and low absolute edge weights in terms of SHD (lower is better). It2-SortnRegres (blue) performance improves with higher I22. sortability, in trend similar to that observed for SID (Section 4.1 and Appendix A.1).
As can be seen in Figure 4, the trends performance for Re-SortnRegress are the same as in terms of SID (Section 4.1 and Appendix A.1). The fact that the results of R² SorrnRReeessss comparatively
15

---

better in terms of SID than in terms of SHD indicates that the algorithm finds a good causal order and inserts excess edges that do not affect the correctness of adjustment sets as measured by SID, The same logic applies to FGES, which consistently performs better than PC in terms of SID (Section 4.1 and Appendix A.1) but does not perform better consistently in terms of SHD (Figure 4).
A.3 Distribution R2 R
In Figure 5 we provide the R sortaabbiyyyyy histograms for the settings described in Section 4 and Appendix A. I as reference points for the distribution of R² sortabibbyy in random graphs sampled with common settings. We observe that the distribution of -o-sortabillty is close to symmetric in ER graphs, while for SF graphs has a strong left skew. The mean values are higher when edge weights are high (Figure 5a) than when they are low (Figure 5b), as is expected from our reasoning about the impact of the weight distribution outlined in Section
200
Count
Count
100
0.2
0.4
0.6
0.8
0.05
0.70
0.75
0.80
0.85
0.90
0.95
1.00
R+-soorability
sorttiiablylty
(a) Higher edge weights: setting from Section with Pw Unif((-1,, u (0.5, 1))
200
Count
Count
200
100
0.2
0.4
0.6
0.8
0.6
0.7
0.8
0.9
1.0
-orttiiability
sortaiibbyy
(b) Lower edge weights: setting from Appendix A.1 with Bv
Unif((-0.5, -0.1) U (0.1,
Figure 5: Histograms and kernel density estimates of the distribution of 1tth /t*=sortabilltyin different simulations.
A.4 Causal Discovery Results For Different Noise Settings
Many identifiability results in the ANM literature focus on the distribution of the additive noise and its variance, and the distinction between different distributions is common in the evaluation of causal discovery algorithms. Here, we present R² RE-SortnRegres results and R² --ortabbilltyy separately for different noise distributions and variances and show results in terms of SID and SHD (note that ff2-SortnRerress is regularized using the Bayesian Information Criterion, which effectively strikes a balance between SID and SHD). We also provide an evaluation of the performance in terms of the discovery of the Markov Equivalence Class (MEC). We transform recovered DAGs into corresponding CPDAGs before calculating the MEC-SID, and additionally transform true DAGs into corresponding CPDAGs before evaluating the MEC-SHD. For each setting we independently sample 30 graphs using the edge weight distribution Pw Unif((-2 -0.5) U (0.5, 2)) and run our algorithms on 1000 iid observations from each ANM.
The causal discovery performances and Rº sortabiilisss of each setting are shown in Figure 6. We find that -?2-sortaillity is consistently high on average (*s 0.8). except for the setting with exponentially distributed noise standard deviations (+= 0.65) shown in Figure 6b. Unlike the uniform distribution used in the other settings, the exponential distribution is positively unbounded, meaning that there is always a chance that the noise standard deviations are large enough to substantially impact R². sortability. In the settings with high average R² --sortbbility, R² ---rrnnngggesss performs similarly to FGES in terms of SID, and somewhat worse in terms of SHD. In Figure 6b, the setting with exponentially distributed a and average Vg² R 0.65, R2 Re-onnnnRessss expectedly performs worse than in the other settings. On raw data, the Top-Down algorithm performs well when the equal noise variance assumption is approximately fulfilled but poorly otherwise (e.g. Figure 6c), and consistently performs similarly worse than RandomRegress on standardized data. The performance of PC and
16

---

FGES is consistent across DAG recovery settings, with only a slight dip for FGES in the case of non- Gaussian noise shown in Figure 6c where its likelihood is misspecified. Conversely, DirectLiNGAM performs well in the setting with non-Gaussian noise, but poorly in the other settings. We find qualitatively similar results for MEC recovery shown in Figure 6d, although the RandomRegress baseline appears to be somewhat stronger, likely because the direction of edges is less impactful for MEC comparisons.
DAG Recovery
Randonnaaraa
NGMM
DreeLLIGMM-
IPC-
PC--
Algoritho
Algorithm
Topp Down
Top Down
Van Pogiie
Van Regrre
10
do
150
550
16o
150
Structural Intervention Distance
Structural Hamming Distance
(a) ANM with Pg GER (20, 40), Px = Ni0, and Po = Unif(0.5, 2); avg. V 122 PS 0.83.
Rendorvilgrres
WCAM
DiossLiNCMM-
Algorithm
Algorthm
PGES
FGES-
Tap
Doen
Van Serttrrrrragrns
100
200
500
20
60
Structural Intervention Distance
Stractural Hamming Distance
(b) ANM with Pg GER 40), (x) = 00, and Po Exp(1.54); avg. Vga 0.65.
landoo
london/aggm
DesULNCAAM-
DisutiNNAM
PC-
PC
Algorithm
Algorithm
FGES
FGEs
Thp Down
Tay Down-
SortSortr
gi Seeteaa
2001
110
Structural Interventian Distance
Stractural Hamming Distance
(c) ANM with Ps = GER (20,40), 40), (o) = Unif(- Lnif(-v3c) 330, v30). and Po Un.f 2); avg. Vgaa 2 0.83.
MEC Recovery
Toadoeellagee
Dewtl INGAMM-
DradLNGAM-
Algorithm
Algorithm
PGES
FGES
Top Down
Top Dovn
Van Sortriiggnna
2o. SertiRegess
MEC Structural Intervention Distance (upper bound)
MEC Structural Hamming Distance
(d) ANM with Pg = GFa (10, 20), PN(0)= = N(0,, )²), and P. = Unif(0.5,, 2): avg. Vn= Ro 0.80.
Figure Causal discovery performance and -ortttbbayy for different Pi and Po.
B R2 -e-sortability for Different Simulation Parameter Choices
B.1 Terminal Node Variance in a Causal Chain Given Weights and Noise Standard Deviations
# Xp)) of length p > with fixed edge weights w for j 0, .... 1 and standard deviations o, for j = 0, ..., P of the independent
In a causal chain (Xo w0.2, X, WL. X₂ 122,3
17

---

noise variables No...., Np. we have that
Var (Xp) = Var (Wip 1, xp-1) + Var (Np) (unfolding the recursion)
2
p-1
p-1
p-1
Var (N₂)
πI Wyst2
+ Var (Np) = Σ
of
W II j
+ top
1=0
p-1
II ==
u Ju+1
2 σ
Σ ==
log |uys+1|--
The application of the log in the final step allows us lower-bound the product by sum and employ the law of large numbers. In Figure we present a simulation of the convergence of R² and the cause- explained variance fraction alongside the increase of variances along the causal order to illustrate our theoretical results given in Section 5.1. We sample edge weights Pw Unif(( 1 -0.5), (0.5, 2)) and noise standard deviations Pa ~ Unif(0.5, 2). We sample 30 chains independently in this way to obtain confidence intervals.
1.0
-0.8
RP CEV-fraction
-0.
R
CEV-fraction
log Var(aanneee)
11
3
12
15
18
Position causal chain
Figure 7: R2, fraction of cause-explained variance (CEV), and variances in simulated causal chains.
B.2
Empirical Evaluation of R²-Sortability in ANMs
In this section we present the Re-sortability obtained for graphs with range of average node in- degrees 7 and weight distributions Pw across random graph models, noise distributions, and noise standard deviation distributions. For each combination of parameters, 122 -fe-sortabllity vaz for 20 independently sampled graphs with d = 50 nodes and 1000 iid observations per ANM.
compute the mean value of
ER and SF Graph Models
20-
20-
15-
15:
05
30
1.0
1.0
V.00
0.1
8.0
00:
02
05-
002
64 13
4.
5.5
75
16
70
50
4.0
50
K00
Average node degree
Average node degrre
(a) ER graphs
(b) SF graphs
Figure 8: Sensitivity of R² -ossorbbbly to parameters in 50-node graphs with Px (O, o²).
P, ~ Unif(0.55, 2) to y and Pw with different E [In [V] with V ~ Pw in ER and SF random graphs.
18

---

Different Noise Distributions and Noise Standard Deviations
10
25-
08
2.0
35
15-
1.5
0.
1.0
100
05-
8.0
.05
-00
155
Average node in-degree
Average node in-degree
(a) Pv (6) (0, x²), P. Onif(0.5, 2)
(b) (o) Unif( v3d, v30). P. Unif(0.5, 2)
100
25
25-
2.0
20-
15
1.9
0.5
V22
05
0.0
032
1.0
10-
04 15
4.5
7.5
0.0
15
45
000
Average node in-degree
Average in-degree
(c) (x) (0, Exp(1.54)
(d) Pv Unif(-V30, V36,, Exp(1.54)
Figure 9: Sensitivity of -fstsoraaillity in (50, 750) graphs to y and Pw with different IE [In [V]] with V Pw, for different noise distributions (r) and noise standard deviation distributions P..
Figure 8 shows comparison of -oottabiiyyyyy in ER and SF graphs. We choose slightly different values of y for the two settings because the SF graph generating mechanism requires integer values. In ER graphs. observe that R sortaabbyyy is high unless E [log [V] or very low. In SF graphs we observe Ro-sortabilities close to across all settings, the relationship between the parameters and R² -sssoraability observed for ER graphs is only faintly visible.
Figure 9 shows the impact of the noise and noise standard deviation distributions on H2 -f--sotabilitty We observe consistent trend of high R2 -Sorttaabyyyty unless E' [log [V] o very low. R stttsiibbss also appear somewhat lower when edge weights and 7 are simultaneously very high. However, find that high E [log [V] and y yield R² extremely close to 1 for many variables. What looks like lower R² -ortaabbyyyy may therefore be due to finite samples and limited numerical precision. R² -orttabbbll does not differ visibly between different noise distributions Pv (ő), but comparing the first to the second row, we observe that exponentially distributed noise standard deviations a leads to lower R²-sortability, provided the E. [log] [ are not too high. This observation is in line with the comparatively low Rd R'-sortability observed for exponentially distributed a in Appendix A.4. It may be explained by the positive unboundedness of the exponential distribution, which enables the noise variances to disrupt R²2 patterns even when cause-explained variance is high.
C Relationship Between the Different 7-Sortabilities
Our use of R2 is motivated by reasoning about increase in the fraction of cause-explained variance (CEV) along the causal order for data with high var-sortability (see Section 3.1). We obtain CEV- sortability using the definition of 7-sortability in Equation (2) with r(X,t) R ( (P(a X),
and denote it as vcEv. Here, analyze the relationship between var-sortability, R and CEV-sortability to gain insight into the co-occurrence of var-sortability (which been shown to be prevalent in many simulation settings, Reisach et al. 2021) and R²- -zo-sttabill well as
19

---

the match between Rº sortabbiliy and CEV-sortability"" For each combination of parameters, we compute the mean sortability of 20 independently sampled graphs with d 50 nodes and 1000 iid observations per ANM.
C.1 Erdós-Rényí Graphs
14
Y
Veey
Veav
44
04
Average node degree
Average node in degree 7
Averge node in degreg 7
(a) Var- and CEV-sortability
(b) I and CEV-sortability
(c) R²- and Var-sortability
Figure 10: Ratio between sortabilities for Pw with different E [In [V] with V Pw and different in (50, 250) graphs with (N (0, Φ²) and P. Unif(0.5, 2).
We see in Figure 10a that CEV-sortability indeed tracks var-sortability closely in most settings as indicated by our analysis of the weight distribution in Section 5. The two measures disagree most when the parameters E [log V and both take very low values. In Figure 10b, R Ro-sortability is generally good proxy for CEV-sortability except when E [log V]] or y are very low. Figure 10c shows that -sotorbbbiy tracks var-sortability well across many settings, unless one of the parameters is very low. The disagreement for very high values of 7 and E [log [V] may be due the difficulty of distinguishing between R (or CEV-fractions) numerically in finite samples as they converge to 1.
that
C.2 Scale-Free Graphs
V Vegy
Vid
Vesy
12
24 50 Average nade in-dggro
10 20 13 49 50 Average node in-cagres
24 32
50 70
Avesage indegrro
(a) Var- and CEV-sortability
(b) R and CEV-sortability
(c) R2- and Var-sortability
Figure 11: Ratio between sortabilities for Pw with different [In [V] with V Pw and different ? in Gsr 750) graphs with PN = (0, Φ²) and P. = Unif(0.5, 2).
We choose slightly different values for compared to Figure 10 because the scale-free graph generating mechanism requires integer values. Figure 11 shows that var-sortability, R2 sortabiility, and CEV-sortability are almost perfectly aligned in all settings. The patterns visible in Figure 10 for ER graphs are only faintly visible for SF graphs, and the alignment between the measures is consistently close.
D Definition of Sortability
The original definition of var-sortability by Reisach al. (2021), given in Equation (2) for r(X,t) Var (Xz), measures the fraction of all directed paths of different length between cause-effect pairs
Note that the range of values and the interpretation of the color scale differs between figures.
20

---

that are correctly ordered by variance. We investigate two alternative definitions of sortability to explore the impact of counting the directed paths between cause-effect pair differently.
D.1 Path-Existence (PE) Sortability
The first alternative sortability definition measures the fraction of directed cause-effect pairs (meaning that there exists least one directed path from cause effect) that are correctly sorted by criterion Since this definition treats all connected pairs equally, irrespectively of the number or length of the connecting paths, we refer as Path-existence sortability' It given as
(PPD (X. E(.-iDEE By ixcr(r(X,s)),.. T(X, where incr( b) = 1/2 a = V
(5)
>
and Σi= Bi is the sum of the matrices obtained by taking the binary adjacency matrix Bg of graph
the i-th power. (s +) t) path from X, to X, in G.
Σ Bi is true for s, è if and only if there is least one directed
D.2 Path-Count (PC) Sortability
The second alternative sortability definition measures the fraction of all directed paths between cause- effect pairs that are correctly sorted by a criterion т. This definition can be understood as weighing cause-effect pairs by the number of connecting paths, and therefore refer it as 'Path-count sortability". is given as
Ei (e (B²),, (x,9))
a<
where incr(a, b) 1/ a=b
(6)
Ei S E( (B%) s.t
and (B%)s,, is the entry at position 8, of the binary adjacency matrix Bg of graph G taken to the i-th power, which is equal to the number of distinct directed paths of length i from X, to X, in g. +) E By is true for s, t and only if there is at least one directed path from X, X, of length d in g.
D.3
Comparison of Sortability Definitions
We conduct experiment compare the definitions of sortability outlined above to the original definition presented in Equation (2) using r(X, t) = Var (X₂) to obtain var-sortability and =(X, t) R2 Ie (M ((Vd) 8) to obtain R²*ssortbbliiyy, as in Section 3. We independently sample 500 ANMs with the parameters Pg = GER (20, 40), (o) N/00, Φ²).. Pa Uni -0.5) (0.5, 2)), and Pw Unif((-0.5, -0.1) U (0.1, 0.5)). From each ANM, we sample 1000 iid observations.
08
Pall-oant
Parh eeieeeeeee
05
60
Origial
Figure 12: Rénói graphs: comparison of the sortabilities obtained as defined by Equation (2) (following the original definition by Reisach et al. 2021) alternative definitions that count directed paths between cause-effect pair differently.
Figure 12 shows comparison of the different definitions of sortability. We can see that the alternatives are highly correlated with the original definition for both Re-sortabilitty and var-sortability in terms of the Kendall rank correlation coefficient (Kendall 1938). The rank correlation coefficient between
21

---

path-existence and path-count sortability is somewhat lower, indicating that the original definition of sortability strikes balance between the two alternatives. The point clouds are also close the diagonal, meaning that the differences in sortability for individual graphs tend to be small. However, for graph structures with different connectivity patters, the differences between the definitions may be more pronounced. To investigate this aspect, we run another simulation in the same fashion, this time using 500 scale-free graphs (Pg = Gsp(20, 40)), and the same parameters otherwise.
Palsoan
Pat
(Poood by: fanda nak comaati owrders 355 Tavorsa rat snvini E90
Rommalin Kendal an constasio orricur D.
Prontd taccal use convatto cotficer 303 Yeoeora kancal rab comaii corfders 179
Vesamality Honall canein cslicau
6.0
Origial
Orignal
Path Fath-easstence
Figure 13: Scale-free graphs: comparison of the sortabilities obtained defined by Equation (2) (following the original definition by Reisach et al. 2021) to alternative definitions that treat the connectivity of cause-effect pair differently.
As can be seen in Figure 13, for scale-free graphs the original definition of sortability is similar the path-count definition, and somewhat less so to the path-existence definition. As before, we observe the biggest difference in terms of rank correlation between path-count and path-existence sortability. While the rank correlations are close to those observed in Figure 12, there is stronger trend of path-count sortability tending be higher, and path-existence sortability tending be lower than those of the original definition. In addition, the trends appear to be different for R -or-sabbayll and var-sortability, particular when comparing path-existence sortability to the other two definitions. Overall, we find that the original definition strikes a balance between the alternatives. Nonetheless, the deviation from the diagonal and difference between R2 sortability and var-sortability may indicate that some definitions could be more appropriate than others depending on the choice of 7 and the assumptions about the graph structure.
E
Controlling CEV Does Not Guarantee Moderate R R²-sortability
One strategy for avoiding var-sortability that may be interesting in the light of R -R-sortability is to enforce equal node variances by assuming constant expla-expl variance and constant exogenous noise variance, thereby also fixing the fraction of cause-explained variance (Squires et al. 2022; Agrawal et al. 2023, Sections 5.
In a causal chain, for the variances and exogenous noise fractions to be equal, the edge weights also have be of equal magnitude. For this reason, in such chain and last node correctly using R², resulting in a -?2-sortability close to if the chain is sufficiently long.
best hope order the first
In general however, the effect of such a sampling scheme on R² -Posorabbyyyty is more complex. Using the parametrization described in Sections 5.1 of Squires et al. (2022), which enforces unit variances and a noise contribution of half of the total variance (except for root nodes, which have only exogenous variance), provide a counterexample show that controlling the fraction of cause-explained variance does not guarantee n22 -Pasottbbilityy close to 0.5. Consider the causal model
X1 = N₁
X2 = X1
+
V2
X3 Xi + X2
N₂
V4+ 2/2
X1 X2 X3
X4
+ No
66 4. +8 1
v2
V4+2v2
22

---

where all N, are independent and the root cause noise Ni ~ ~V, 1), while N2. Ns. N4 ~ NN,, 2) such that all variables other than the root cause have cause-cxplained variance fraction of +. and all variables have a marginal variance equal to I, as in Section 5.1 of Squires et al. (2022). We are interested in the nº² (2,,t) for all & E 21, 4J. We can use partial correlations to analytically compute
R (X,t) =1 (2 (2-
where >' the covariance matrix of X, which can be obtained analytically by computing the pairwise covariances using the definition of the random variables. Using the definitions
(1 a)) (1
1 2+v2
1+a
V442V2
V4+2V2
B=
7
V2
44+2V2
v6 +
83
v6 + 1a 88
we can write the covariance matrix of (X1, X2, X4) as
[
8
7
=
a
8 1
7
3
3
8
2
7
8
1
Note that a > 3> 7 > 8. which already hints at a possible ordering of the variables by R². Computing the R² yields (x, = 0.59. R2 (X, 2) x 0.59, (x, 3) 0.53, and R²(X, 4) 0.5. This means that X, (the root node) and X2 have the same R² while, conversely and perhaps surprisingly, the variables X2, X3. and X4. which all have the same fraction of cause-explained variance, have different and descending R². Even X2 and X3 have different π² despite having the same fraction of eapainex variance, and both being neither root cause -Rssorabbility of this example is 1/22. In finite sample settings where the equality of the T22 of X, and X2 is not exact, sorting by increasing R² is thus likely give causal ordering that is either equal very close to the inverse of the true causal order. The inverse causal order is markedly different from random ordering (v 0.5), and shows that the R² carry information about the causal ordering. As this example shows, controlling the CEV does not guarantee a -f@-sortability close to 0.5, and we recommend evaluating and reporting R -orttabbbyy even when node variances or the cause-explained variance fractions are controlled.
leaf node. The
23


Reviews:
[{"review_id": "tCkjD10dqf", "reviewer": "Reviewer_j7gp", "summary": "This paper proposes an interesting extension of the var-sortability approach proposed by Reisach et al. (2021), denoted R2-sortability. While var-sortability measures the agreement between the causal order and the order of increasing marginal variance, R2-sortability measures the agreement between causal ordering and the order of the explainable fraction of a variable’s variance (as captured by the coefficient of determination R2). Contrary to var-sortability, the R2-sortability approach is a scale-invariant approach which can be applied even when the data scale is arbitrary. \n\nThe paper proposes a baseline algorithm, “R2-SortnRegress” which can match, and even exceed, the performance of established causal discovery methods such as the PC and FGES algorithms. Comparisons against these baselines (and the Var-SortnRegress and RandomRegress baselines) were performed using synthetic data generated from random and scale-free graphs using Gaussian noise data (with noise standard deviations and edge weights i.i.d. sampled from uniform distributions). The paper also presents a brief real data comparison using a single dataset.\n\nThe paper provides a mathematical analysis of the influence of the weight distribution on R2 in the simple setting of causal chains, as well as, an empirical evaluation on more complex random graphs, underscoring the role of the weight distribution as a driver of R2.  \nThe analyses presented in the paper illustrate the need to make decisions on all ANM parameters in simulations. \n", "questions": "Why did the paper only select the PC and FGES as traditional causal discovery algorithms in their comparisons? (The Reisach et al. (2021) paper includes other popular methods.)\n\nAlso, why the synthetic data experiments only report SID, while the real data illustration reports SID and SHD? \n", "limitations": "Yes.", "rating": 7, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "This is a well written and interesting paper. The approach is novel and its need is well motivated. The proposed approach is sound and well described in the text. ", "weaknesses": "I only have a few very minor comments. The notation is sometimes a little inconsistent. For instance, $\\sigma$ and $\\phi$ are used interchangeably to represent the standard deviation. Also, in Algorithm 1, $\\sigma$ is also used to represent the permutation that sorts $\\pi$ in ascending order.\n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank you for your thorough and positive review.\nWe propose to make the following edits in response to your review and answer your questions individually below.\n\n---\n\n### Edit summary\n\n* **Edit 1** We will include DirectLiNGAM (Shimizu 2021) as an additional algorithm in our comparisons in Appendix A3, Figure 5 where we consider non-Gaussian noise, among other settings.\n* **Edit 2** We will add SHD results for the settings shown in Figure 1 and Figure 3 as a new subsection in Appendix A2 and refer to them in Section 4.1.\n* **Edit 3** We will replace $\\sigma$ in the algorithm and align the use of $\\phi$ and $\\sigma$.\n\n---\n\n### Answers to your questions\n\n**Comparison to other algorithms**\n\nFollowing the findings by Reisach [2021] we decided to include Var-SortnRegress as a baseline representative of the performance of scale-sensitive algorithms given the large and growing number of such algorithms (Vowels 2021).\nIn Section A3 we also include the _Top-Down_ algorithm as an additional method that exploits another type of assumption (equal noise variances), which performs largely similar to _Var-SortnRegress_. \\\nWe will further add the DirectLiNGAM algorithm as described in **Edit 1**. \nDirectLiNGAM highlights the specificity of the non-Gaussian setting by outperforming the other methods for non-Gaussian noise but performing badly otherwise, in line with theory.\nWe may be able to include another algorithm for the camera-ready version of the paper if there is one that would provide a substantially different aspect or perspective and thus help highlight the main takeaway of the comparison?\n\n**Addition of SHD results (qualitatively the same as SID results)**\n\nThank you for the suggestion. We agree that the SHD results may be interesting to readers and will add them as described in **Edit 2**.\nThe trends in terms of SHD are the same as in the SID results. $R^2$-SortnRegress improves with increasing $R^2$-sortability, and absolute performances are better on SF graphs than on ER graphs.\n\n---\n\n### Additional points\n\n**Notation** — Thank you for pointing this out, we will address it (cf. **Edit 3**)."}, {"review_id": "FQH9V9UKDY", "reviewer": "Reviewer_iqwe", "summary": "This paper is concerned with synthetic data generation for causal discovery. The paper explores the scale invariant pattern given by the coefficient of determination $R^2$ that potentially exists in synthetic data benchmarks. The authors present analysis in the case of linear ANMs. They also find out that prior over linear parameters might have significant influence on the $R^2$ sortability.", "questions": "1. Is it possible to generate the synthetic data by first fixing the desired $R^2$, and generate the synthetic data which exactly matches this value?\n2. Could you comment on whether this is a potential issue beyond ANMs and linear models?", "limitations": "The authors have done a clear job of stating the assumptions under which their conclusions are valid.", "rating": 6, "confidence": 4, "soundness": 3, "presentation": 4, "contribution": 3, "strengths": "(++) There are many new algorithms being developed for causal discovery, and all of them evaluate on synthetic data as evaluation on real world data is very hard. Hence, any new insights which helps bridge the gap between synthetic data and real data is very timely and welcome.\n(++) This work generalises the results of Reisach et al (2021) in giving a general `sortability` criterion. Further, the $R^2$ sortability is more subtle than varsortability as it cannot be just adjusted by renormalisation.\n(+) I like the overall presentation and motivation of the problem. The presented analysis for linear ANM and trees is insightful.", "weaknesses": "(--) It is not clear how big of an issue is $R^2$ sortability in real world settings. Section 4.2 illustrates it on Sachs dataset. But Sachs dataset is not a linear ANM, and hence has significant model misspecification. So I am unsure how concrete are the conclusions from 4.2. Having said that, I still find the overall contribution useful. For example, if it is the case that there is a $R^2$ sortability issue in a more relevant dataset, one could at least find out to what extent it exists based on the insights of this paper.\n(--) A theoretical analysis would have been more helpful. Right now, it is not clear if this is just a problem in linear ANM or can be generalised to nonlinear/nonparametric cases. More so, it might have made more clear whether it is actually the linear parameters which affect it the most in non tree situations (which would be the case in most practical settings).\n(--) Unlike the varsortability, which is clear shown that almost all the new causal discovery algorithms being developed exploit it, it is not necessarily the case here. For example, it is not clear whether new optimisation based/ neural network based methods exploit this to perform better. I also don['t think it is necessarily a bad thing to exploit it if it is the case that $R^2$ sortability exists in real world settings (it is yet unclear from this paper, see above point).", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your thoughtful review and for pointing out the impact of our findings on the causal structure learning community. We propose to make the following changes in response to your review and answer your questions individually below.\n\n---\n\n### Edit summary\n\n* **Edit 1** We will highlight the open nature of real-world $R^2$-sortability by adapting lines 19-21 (Abstract): \\\n_Our findings reveal high $R^2$-sortability as an assumption about the data generating process relevant to causal discovery and implicit in the choice of simulation parameters. It should be made explicit, as its prevalence in real-world data is an open question._\n  \n* **Edit 2** We will emphasize the impact of $R^2$-sortability on benchmarking practices 84 (Contribution): \\\n_Characterizing $R^2$-sortability and uncovering its role as a driver of causal discovery performance enables distinguishing between data with different levels of $R^2$-sortability when benchmarking or applying causal discovery algorithms to real-world data._\n\n* **Edit 3** We will emphasize the benefit a full theoretical characterization of $R^2$-sortability could provide in line 337 (Discussion): \\\n_A complete theoretical characterization of the conditions sufficient and/or necessary for extreme $R^2$-sortability could, if found, help decide where one may hope to exploit it in practice._\n\n---\n\n### Answers to your questions\n\n**Generating data with a specific $R^2$-sortability**\n\nBecause $R^2$-sortability arises in a complex interplay of simulation parameters there is no systematic way to achieve a target $R^2$-sortability when sampling ANM parameters iid that we know of. \nOne may use rejection sampling or update the properties of a given graph to achieve a target $R^2$-sortability with non-iid edge weights (cf. Sections 5 in Agrawal [2021], Squires [2022], as pointed out by reviewer zF5c).\n\n**$R^2$-sortability beyond ANMs and linear models?**\n\nWe believe that $R^2$-sortability can be an issue beyond ANMs and linear models, although its impact may vary depending on function class and parameters, and increasing functional complexity along the causal order could reduce or flip $R^2$-sortability.\nIn practice, estimating $R^2$ for nonlinear models requires the choice of a suitable regression method.\nPrompted by your question, we ran an experiment using our linear $R^2$-sortability as a first-order approximation on two different nonlinear Gaussian Process simulations used by Zheng [2020] and Reisach [2021].\nLinear $R^2$-sortabilities range from 0.10 to 0.85, indicating that extreme values are possible in non-linear ANMs and making this question an interesting topic for future work.\n\n|graph|SEM|Var-sortability|$R^2$-sortability|\n|---|---|---|---|\n|ER(20, 20)| Additive GP|0.71$\\pm$0.13|0.49$\\pm$0.18|\n||GP| 0.66$\\pm$0.14| 0.37$\\pm$0.13|\n| ER(20, 80)|Additive GP | 0.88$\\pm$0.06|0.41$\\pm$0.10|\n||GP| 0.61$\\pm$0.13| 0.10$\\pm$0.06|\n|SF(20, 20)| Additive GP| 0.83$\\pm$0.11|0.85$\\pm$0.10|\n||GP| 0.66$\\pm$0.16| 0.52$\\pm$0.17|\n| SF(20, 80) | Additive GP| 0.96$\\pm$0.02|0.75$\\pm$0.14|\n||GP|0.66$\\pm$0.11| 0.23$\\pm$0.11|\n\n----\n\n### Additional points\n\n**Real-world $R^2$-sortabilities**\n\nWe agree that real-world $R^2$-sortabilities are an important open question that cannot be answered by any single dataset (see line 336), and will communicate this more clearly (see **Edit 1**).\nReal-world $R^2$-sortabilities are not known and that is precisely why we think it is important to be aware that current simulations systematically tend to result in high $R^2$-sortabilities.\nAny in causal discovery (e.g. Gaussian vs non-Gaussian, equal vs unequal noise variances, etc.) may or may not hold on different real-world data, which is why they are treated separately in simulations.\nFor benchmarking, choosing ANM parameters that result in high $R^2$-sortability amounts to an assumption, which our contribution makes explicit.\nThis allows a distinction between data with different $R^2$-sortabilities in the same way we already distinguish between other assumptions to reflect different possible real-world scenarios.\nWe agree that this should be spelled out and will do so (see **Edit 2**).\n\n**Theoretical analysis**\n\n$R^2$ depends on a complex interplay of weights, noise parameters, and graph structure, which greatly complicates the development of a general theoretical characterization of $R^2$-sortability for generic and flexible parameter sampling schemes.\nWhile it is outside the scope (and page limit) of our contribution, we strongly share your view that more research on the conditions for the emergence of $R^2$-sortability would be beneficial, and hope that sharing our work with the NeurIPS community will inspire research leading to new theoretical insights. \nWe will state the potential benefit of such a result through **Edit 3**.\n\n**Meaning of $R^2$-sortability for existing algorithms**\n\nAs a property of the data, $R^2$-sortability affects the baseline performance that can be achieved and is thus relevant to the interpretation of the results obtained by causal discovery algorithms (we do not claim that other causal discovery algorithms exploit $R^2$-sortability).\nFor example, the same causal discovery performance may be less impressive if the $R^2$-sortability on a dataset is 0.95 rather than 0.5."}, {"review_id": "ItWJcRjocH", "reviewer": "Reviewer_zF5c", "summary": "The paper introduces the issue of \"$R^2$-sortability\" for synthetically generated data used in the evaluation of causal structure learning methods. $R^2$-sortability is a generalization of varsortability, which is invariant to re-scaling (e.g., standardizing) the simulated variables. They show that, using typical simulation parameters, sorting variables based on $R^2$-sortability can give good performance on causal discovery. Finally, they investigate how $R^2$-sortability is influenced by the data generation parameters.", "questions": "Please address the issues on clarity raised in the **Weaknesses** section.\n\n### Other suggestions\nThe authors may find it interesting that previous papers have designed a data generation process to control cause-explained variance, see eg. [1], Section 5 and [2], Section 5. Readers would likely find it helpful if the authors described such a data generation process to ameliorate the issue of $R^2$-sortability.\n\n[1] Agrawal, R., Squires, C., Prasad, N., & Uhler, C. (2021). The DeCAMFounder: Non-linear causal discovery in the presence of hidden variables.\n[2] Squires, C., Yun, A., Nichani, E., Agrawal, R., & Uhler, C. (2022). Causal structure discovery between clusters of nodes induced by latent factors.", "limitations": "The authors have adequately addressed the limitations of their work.", "rating": 7, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "**Clarity:** The paper clearly describes the motivation for introducing the concept of $R^2$-sortability, and the experimental results section clearly describes their experimental setup.\n\n**Significance:** There is significant value in papers which describe problems with current methods of evaluation. In particular, since evaluating causal structure learning methods relies heavily on synthetic data, it is important to characterize \"artifacts\" of synthetic data generation procedures, and make efforts to remove these artifacts. It is also helpful that the authors do some investigation into how $R^2$-sortability is influenced by the synthetic data generation parameters.", "weaknesses": "### Experimental Results\nIt is surprising that, in Section 4.1, the authors do not standardize the data so that *Var-SortnRegress* would perform poorly. It is not clear what message the authors are trying to convey by their experimental setup. I would expect that they would want to show that $R^2$-sortability is still an issue in simulated data, even when varsortability is not an issue. \n\n### Minor issues\n*Clarity*: Clarity could be improved in some places. \n1. Equation (3) seems to involve the number of paths between pairs of nodes $s$ and $t$, I think it would be more intuitive to write it this way and to also describe why we care about the number of paths rather than just the existence of a path.\n2. Equation (3) never decreases when adding an edge. It seems undesirable that the score would always (weakly) favor denser graphs. I see that there is a sparsity penalty in Algorithm 1 - why is it introduced there, instead of earlier?\n3. In Section 5, why do you switch from the product of squared edge weights to sum of the logs? It seems motivated by the appeal to the strong law of large numbers, but since you already lower bound the product of squared edge weights by the sum of the logs, then almost sure convergence also holds for the product. Right now, the logs seem ad hoc / unmotivated.\n4. When introducing the condition on $\\mathbb{E}(\\log |V|)$ in Equation (5), you should make it more clear what the conditions means for $P_W$. For example, $P_W = Unif([0.5, 1])$ would not satisfy the condition, would $P_W = Unif([0.5, 2])$ satisfy it?\n5. Please use $\\mathbb{E}$ for expectations, e.g. in Equation (5). This is typically preferred, but especially important when dealing with graphs where $E$ often denotes edges.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank you for your detailed review and for highlighting the significance of our contribution for the causal structure learning community. We propose to make the following edits in response to your review and answer your questions individually below.\n\n---\n\n### Edit summary\n\n* **Edit 1 (see Figure 1 in rebuttal PDF)** We have revised Figure 1 (and 3 & 4 likewise) to better show the strengths of our method. We will show results on standardized data to emphasize the scale-invariance of our method. For reference, we visually differentiate the baseline RandomRegress and include Var-SortnRegress on raw data as a dashed line. To simplify the interpretation of trends and uncertainty we will show moving averages using a window of width 0.1 (instead of binning by decile), and show error bars for the 95% confidence interval of the mean.\n\n* **Edit 2** We will add the following explanation to Equation (3) (line 168): \\\n_In effect, this is the fraction of directed paths of unique length between any two nodes that satisfy the sortability condition in the numerator._\n\n* **Edit 3 (see Figure 3 in rebuttal PDF)** We will add an appendix section showing an empirical comparison of sortabilities obtained using Equation (3) as is, compared to a version of Equation (3) that only considers the existence of a path. In the experiment we observe strong alignment (Kendall rank coefficients of $0.86$ for $R^2$-sortability and $0.84$ for var-sortability) between the two versions.\n\n* **Edit 4** We will highlight the novelty and usefulness as condition for the divergence of node variances of Equation (5) (line 274) and will motivate the use of the log (line 548): \\\n274: _We introduce the following sufficient condition for weight distributions to result in diverging node variances along causal chains of increasing length:_ \\\n548: _The application of the log in the final step allows us to lower-bound the product by a sum and employ the law of large numbers._\n\n* **Edit 5** We will include the references you suggested in the first paragraph of Section 5 (line 259): \\\n_$R^2$-sortability can be mitigated by introducing the assumption of a constant exogenous noise fraction (see for example Agrawal 2021, Squires 2022), which requires non-iid edge weight sampling in simulations._\n\n---\n\n### Answers to your questions\n\n**Experimental results**\n\nWe agree that the scale-invariance of $R^2$-sortability and -SortnRegress should be emphasized more and will make the changes described in **Edit 1**.\n\n**Minor issues (clarifications)**\n\n*1.* We will add an intuitive description of Equation (3) as outlined in **Edit 2**. Equation (3) builds on the definition by Reisach [2021] in a way that yields their original version of var-sortability as a special case. To address your question, we ran a comparison between the original sortability definition counting paths of different lengths compared to one that only considers the existence of a path and find no qualitative differences. We will include the experiment comparing the sortability variants as described in **Edit 3**.\n\n*2.* Sortability is a measure of a given data generating process and its causal structure. For example, for \n\n$A=N_A$, $\\quad B=10 A + N_B$,  $\\quad C=0.5 B + N_C \\quad$ with graph G1 $(A \\to B \\to C)$ and $\\text{Var}(N_A)=\\text{Var}(N_B)=\\text{Var}(N_C)=1$ we have $\\mathbf{v}_\\text{Var}((A,B,C), G1) = 2/3$.\n\nFor the data generating process \n\n$D=N_D$, $\\quad E=10 D + N_E$, $\\quad F=N_F \\quad$ with $\\text{Var}(N_D)=\\text{Var}(N_E)=\\text{Var}(N_F)=1$ and graph G2 $(D \\to E \\quad F)$ we have that $\\mathbf{v}_\\text{Var}((D,E,F), G2) = 1$.\n\nAs this example shows, a data generating process with denser graph can have lower sortability than a DGP with sparser graph (and vice-versa). \\\nAlgorithm 1 is used to estimate a causal graph from data and the sparsity constraint avoids returning a fully connected graph: after having sorted the variables we regress each node onto its predecessors in this order, and the sparsity constraint prunes superfluous edges.\n\n*3.* The application of the log allows us to make our statement about the divergence of node variances, and Equation 5 introduces a sufficient condition for divergence (the related conditions $\\mathbb{E}|V|>1$ or $\\mathbb{E}V^2>1$ are not equivalent to our condition as can be seen by the example of $\\text{Unif}(0.1, 2)$).\nIn particular, if Equation 5 is satisfied then this implies that for increasing chain lengths the sum of log-weights converges to $+\\infty$,\nwhich in turn implies that the product of squared weights converges to $+\\infty$.\nIf Equation 5 is not satisfied, another condition would be needed to establish convergence of the product to $+\\infty$.\nWe believe the motivation for applying the log and the role of equation (5) can be made clearer by implementing **Edit 4**. Thank you for pointing this out!\n\n*4.* We will add the values for $P_W = \\text{Unif}((-2, -0.5)\\cup(0.5, 2))$ ($\\approx 0.16$) in section 4.1, and for $P_W = \\text{Unif}((-0.5, 0.1)\\cup(0.1, 0.5))$ ($\\approx -1.29$) in Section A2.1, both with a pointer to Section 5. (Figure 2 shows further weight ranges.)\n\n*5.* Agreed. We will change the notation for the expectation from $E$ to $\\mathbb{E}$ throughout the manuscript.\n\n**Suggestions**\n\nThank you for pointing out the data generating mechanisms in these related works. We will include them as described in **Edit 5** and hope they can serve as a reference point to develop new and more realistic simulation schemes."}, {"review_id": "WBHzESvlUg", "reviewer": "Reviewer_UFwV", "summary": "The paper considers causal inference in the context of additive noise models (ANMs). In particular the paper points out a possible problem in simulation benchmarks for this setup: if the weights of the causal network are not chosen appropriately, the simulation may result in datasets in which identifying the causal ordering can be done with a trivial method which simply sorts the variables by the R^2 when the variable is predicted with other variables. Previously, a similar result has been published by ordering the variables by their variance, but this paper shows that the same outcome can be achieved by sorting the variables by the amount of variance explained, R^2, which is scale-invariant, and hence the problem can't be resolved by simple rescaling.\n", "questions": "Could you show a histogram of R^2 sortability values among simulated datasets on which Fig. 1 was based? Now the R^2 sortability is shown on the x-axis but it is not clear what are the relative amounts of different R^2 values in the simulations.", "limitations": "Yes.", "rating": 6, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "The clarity of the paper is good. The paper demonstrates possible deficiencies in existing benchmarks, which is important. Overall, I found this paper useful, but I also think the significance could be strengthened (see below).", "weaknesses": "The paper demonstrates that simulated datasets may have a large R^2 sortability. However, it does not go very far in clarifying how big a problem this is in practice, i.e., it is not really clear how R^2 sortable real-world datasets actually are (there's one example but that's still quite limited). If the real-world datasets actually often are R^2 sortable, then there is no problem. If, on the other hand, they are not R^2 sortable, it would have increased the significance of the paper to provide more concrete suggestions about how simulations should be conducted to in order for them to better correspond to real-world datasets. In real-world datasets the ground-truth may not be easily available, but at least it is easy to calculate the distribution of R^2 values for all variables in a dataset, for multiple datasets, and investigate how the weights in the causal graph simulator should be chosen for the R^2 distribution to be similar to real-world datasets.\n\nThe take-home message from the experiments is not always clear. In Fig. 1 the presented method R^2-SortnRegress seems always worse than the previously published Var-Sortnregress. Also, the conclusion in the caption: \"R^2Sortnregress performs well it R^2 sortability is high\" seems a bit tautological. It is clear that R^2 is scale-invariant unlike the variance based criterion, so why not simulate datasets that accordingly demonstrate its strength? Also, in the real-world data the results don't seem impressive, as the novel R^2 SortnRegress has the worst SID value. (What is the othe value SHD? I didn't find it defined.)\n\nThe theory considers a causal chain, where E(log|V|)>0, where V is simulated from the weight distribution, and shows that in this case the variance of the last node in the chain goes to infinity and the amount of variance explained converges to one when the length of the chain increases. It seems that if the expected value of |V| is smaller than 1, then the condition is not satisfied? Indeed, it is easy to imagine that the chain will diverge if we have coefficients that in general are larger than one, but this does not sound very realistic assumption in real-world data sets. It would have been interesting to investigate the distribution of estimated weights in some real-world datasets, to see how commonly this condition holds true in those. The paper says that if the condition did not hold, \"detecting statistical dependencies between variables connected by long paths\" would be difficult, but I would imagine this to be case with many real-world datasets, so this does not seem a proper reason for not considering such cases.\n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your constructive review and your suggestions on how to better present the strengths of our work.\nWe propose the following edits in response and individually answer your questions below.\n\n---\n\n### Edit summary\n\n* **Edit 1 (see Figure 2 rebuttal PDF)** We will include the $R^2$-sortability histograms for the settings shown in Figure 1 and 3 in Section A2. For ER graphs, the distribution of $R^2$-sortabilities is close to symmetric, while for SF graphs it has a strong left skew.\n\n* **Edit 2** We will highlight the open nature of real-world $R^2$-sortability (lines 19-21): \\\n_Our findings reveal high $R^2$-sortability as an assumption about the data generating process relevant to causal discovery and implicit in the choice of simulation parameters. It should be made explicit, as its prevalence in real-world data is an open question._\n\n* **Edit 3** We will emphasize the impact of $R^2$-sortability on benchmarking practices (line 84): \\\n_Characterizing $R^2$-sortability and uncovering its role as a driver of causal discovery performance enables distinguishing between data with different levels of $R^2$-sortability when benchmarking or applying causal discovery algorithms to real-world data._\n\n* **Edit 4 (see Figure 1 in rebuttal PDF)** We have revised Figure 1 (and 3 & 4 likewise) to better show the strengths of our method. We will show results on standardized data to emphasize the scale-invariance of our method. For reference, we visually differentiate the baseline RandomRegress and include Var-SortnRegress on raw data as a dashed line. To simplify the interpretation of trends and uncertainty we will show moving averages using a window of width 0.1 (instead of binning by decile), and show error bars for the 95% confidence interval of the mean.\n\n* **Edit 5** We will clarify the take-home message of the real-world experiment (line 254): \\\n_This indicates that, on real-world data, we may not expect to see consistently high $R^2$-sortabilities as much as we do in many simulation settings (see Appendix B.2).\nFor benchmarks to be representative of this factor, they should differentiate between settings with different levels of $R^2$-sortability._\n\n* **Edit 6** We will highlight the novelty and usefulness as condition for the divergence of node variances of Equation (5) (line 274): \\\n_We introduce the following sufficient condition for weight distributions to result in diverging node variances along causal chains of increasing length:_\n\n---\n\n### Answers to your questions\n\n__Histogram of $R^2$-sortabilities (see Figure 1 in rebuttal PDF)__\n\nThis is a great idea, thank you. We will include the histograms as described in **Edit 1**.\n\n---\n\n### Additional points\n\n**Real-world $R^2$-sortabilities**\n\nWe agree that real-world $R^2$-sortabilities are an important open question that cannot be answered by any single dataset (see line 336), and will communicate this more clearly (see **Edit 2**).\nReal-world $R^2$-sortabilities are not known and that is precisely why we think it is important to be aware that current simulations systematically tend to result in high $R^2$-sortabilities.\nAny in causal discovery (e.g. Gaussian vs non-Gaussian, equal vs unequal noise variances, etc.) may or may not hold on different real-world data, which is why they are treated separately in simulations.\nFor benchmarking, choosing ANM parameters that result in high $R^2$-sortability amounts to an assumption, which our contribution makes explicit.\nThis allows a distinction between data with different $R^2$-sortabilities in the same way we already distinguish between other assumptions to reflect different possible real-world scenarios.\nWe agree that this should be spelled out and will do so (see **Edit 3**).\n\nYou suggest calculating $R^2$ distributions across multiple datasets and investigating how to match simulations to those values, yet real-world data suiting ANM assumptions are scarce, hence the extensive use of simulations in causal discovery.\nSuch a study would have to find datasets, suitable function classes for $R^2$ estimation, and a strategy to adapt simulations (one could change the weights, graph structure, or noise variances) consistent with domain knowledge.\nWhile this is outside the scope of our work, we share your excitement for such a project and hope our work will inspire studies like these when shared with the community.\n\n**Take-home message of experiments**\n\nWe agree that the scale-invariance of $R^2$-sortability should be emphasized more and update accordingly (see **Edit 4**). Thank you for highlighting this. \\\nWe agree that none of the performances on the real-world dataset are impressive (line 250), which underlines the need for identifying potentially unrealistic benchmark properties such as high $R^2$-sortability. To clarify the take-home message, we will implement **Edit 5** and move the definition of the Structural Hamming Distance (SHD) to the main text.\n\n**Case of $E\\log|V|<0$**\n\nThank you for raising this.\nWe can better clarify the role of the condition $E\\log|V|>0$, which we prove is sufficient for the divergence of node variances to infinity (see **Edit 6**).\nIndeed, our condition is not satisfied if $E|V|<1$ (Jensen's inequality).\nConversely,\n$E|V|>1$ does not imply $E\\log|V|>0$ (e.g. consider $\\text{Unif}(0.1, 2)$),\nand is thus not sufficient for the divergence.\n\nWe agree that real-world processes may not have weights that result in diverging variances. \nThis is why we think it is important to point out that simulations often rely on parameters that result in high var- and $R^2$-sortability.\nUnfortunately, measuring weights in real-world ANM data (which are scarce) is not straightforward, since they are affected by the arbitrariness of the data scale. \\\nWe do consider the case of $E\\log|V|<0$ in our empirical analyses.\nFigure 2 shows that even $E\\log|V|<0$ result in high $R^2$-sortabilities for Scale-free graphs,\nand Section A2.1 shows a setting with small weights and $E\\log|V|<0$."}]

Author rebuttal:
Dear All,

We thank you for the review of our submission, valuable suggestions, and recognition of our work's significance and timeliness for the causal discovery community.
All reviews seem to agree on the good technical soundness (4x good)
and presentation (1x excellent, 3x good) of the overall contribution (3x good, 1x fair) our paper makes by characterizing a novel pattern implicit in the parameters of additive noise models that can strongly affect causal discovery.

We have revised our manuscript carefully and made improvements based on all reviewers' comments:

* _Updated Figures._ We have improved existing figures and created new figures to complement the presentation of our experimental results and better highlight some of the strengths of our work (scale-invariance of $R^2$-sortability and -SortnRegress, distribution of $R^2$-sortabilities, robustness to different ways of counting sortable node pairs).
* _Revisions._ We have rephrased or added individual sentences to clarify explanations and include two new references.

We summarize these edits below and provide detailed responses to all reviewers in our individual replies.

We are looking forward to discussing our work with you and hope our replies answer your questions and help in your assessment of our submission. 
We will continue to revise our manuscript if further suggestions arise during the author-reviewer discussion.

---

__Updated Figures__

- __see Figure 1 in attached PDF__ [zF5c, UFwV] (Figure 1, Section 4.1). \
We have revised Figure 1 (and 3 & 4 likewise) to better show the strengths of our method. We will show results on standardized data to emphasize the scale-invariance of our method. For reference, we visually differentiate the baseline RandomRegress and include Var-SortnRegress on raw data as a dashed line. To simplify the interpretation of trends and uncertainty we will show moving averages using a window of width 0.1 (instead of binning by decile), and show error bars for the 95% confidence interval of the mean.

- __see Figure 2 in the attached PDF__ [UFwV] (Appendix A2). \
We will include the $R^2$-sortability histograms for the settings shown in Figure 1 and 3 in Section A2. For ER graphs, the distribution of $R^2$-sortabilities is close to symmetric, while for SF graphs it has a strong left skew.

- __see Figure 3 in the attached PDF__ [zF5c] (new Appendix C). \
We will add an appendix section showing an empirical comparison of sortabilities obtained using Equation (3) as is, compared to a version of Equation (3) that only considers the existence of a path. In the experiment we observe strong alignment (Kendall rank coefficients of $0.86$ for $R^2$-sortability and $0.84$ for var-sortability) between the two versions.

- [j7gp] (Figure 5, Appendix A3). \
We will include DirectLiNGAM (Shimizu 2021) as an additional algorithm in our comparisons in Appendix A3, Figure 5 where we consider non-Gaussian noise, among other settings.
We will also add the SHD results corresponding to Figures 1 and 3 in A2 (the trends are the same as for SID).


__Revisions__

* We will highlight the open nature of real-world $R^2$-sortability [iqwe, UFwV] (lines 19-21, Abstract). \
_Our findings reveal high $R^2$-sortability as an assumption about the data generating process relevant to causal discovery and implicit in the choice of simulation parameters. It should be made explicit, as its prevalence in real-world data is an open question._

* We will emphasize the impact of $R^2$-sortability on benchmarking practices [iqwe, UFwV] (line 84, Contribution). \
_Characterizing $R^2$-sortability and uncovering its role as a driver of causal discovery performance enables distinguishing between data with different levels of $R^2$-sortability when benchmarking or applying causal discovery algorithms to real-world data._

* We will add the following explanation to Equation (3) [zF5c] (line 168, Section 3.1). \
_In effect, this is the fraction of directed paths of unique length between any two nodes that satisfy the sortability condition in the numerator._

* We will clarify the implications of our findings on real-world data for simulations [UFwV] (line 254, Section 4.2). \
_This indicates that, on real-world data, we may not expect to see consistently high $R^2$-sortabilities as much as we do in many simulation settings (see Appendix B.2).
For benchmarks to be representative of this factor, they should differentiate between settings with different levels of $R^2$-sortability._

* We will include the suggested references [zF5c] (line 259, Section 5). \
_$R^2$-sortability can be mitigated by introducing the assumption of a constant exogenous noise fraction (see for example Agrawal 2021, Squires 2022), which requires non-iid edge weight sampling in simulations._

* We will highlight the novelty and usefulness as condition for the divergence of node variances of Equation (5) (line 274) [zF5c, UFwV] (line 274, Section 5.1). \
_We introduce the following sufficient condition for weight distributions to result in diverging node variances along causal chains of increasing length:_

* We will emphasize the benefit a full theoretical characterization of $R^2$-sortability could provide [iqwe] (line 337, Section 6). \
_A complete theoretical characterization of the conditions sufficient and/or necessary for extreme $R^2$-sortability could, if found, help decide where one may hope to exploit it in practice._

Meta-review: This paper is concerned with finding causal order among the given variables according to the additive noise model. It demonstrates a possible R2 sortability problem in simulation benchmarks. The finding is interesting and has immediate implications for designing simulations to evaluate causal discovery methods.
Decision: Accept (poster)