Causal Feature Selection via Orthogonal Search
Authors: Ashkan Soleymani
Anant Raj
Stefan Bauer
Bernhard Schölkopf
Michel Besserve
Venue: Accepted by TMLR ()

Abstract: The problem of inferring the direct causal parents of a response variable among a large set of explanatory variables is of high practical importance in many disciplines. However, established approaches often scale at least exponentially with the number of explanatory variables, are difficult to extend to nonlinear relationships and are difficult to extend to cyclic data. Inspired by debiased machine learning methods, we study a one-vs.-the-rest feature selection approach to discover the direct causal parent of the response. We propose an algorithm that works for purely observational data while also offering theoretical guarantees, including the case of partially nonlinear relationships possibly under the presence of cycles. As it requires only one estimation for each variable, our approach is applicable even to large graphs. We demonstrate significant improvements compared to established approaches.

Published in Transactions on Machine Learning Research (08/2022)
Causal Feature Selection via Orthogonal Search
Ashkan Soleymani"
ashkanso mit edu
Department of Electrical Engineering and Computer Science Massachusetts Institute of Technology. Cambridge, US
Anant Raj
rnant.ruj@einia.ff
Inria, Ecole Normale Supérieure
PSL Research University, Paris, France
Stefan Bauer
baue @ kth.se
KTH, Stockholm, Sweden
Bernhard Schölkopf
bs @ tuebingen.mpg.de
Planck Institute for Intelligent Systems Tübingen, Germany
Michel Besserve
michel.besserveeeuehinngg
Max Planck Institute for Intelligent Systems Tübbingen, Germany
Reviewed on OpenReview:
Abstract
The problem inferring the direct causal parents of a response variable among a large set of ex- planatory variables is of high practical importance in many disciplines. However, established ap- proaches often scale least exponentially with the number of explanatory variables, difficult to extend nonlinear relationships and difficult to extend cyclic data. Inspired by Debiased machine learning methods, we study one-vs. -ss--tt--rstt feature selection approach to discover the direct causal parent of the response. We propose an algorithm that works for purely observational data while also offering theoretical guarantees, including the case of partially nonlinear relationships possibly under the presence cycles. As it requires only one estimation for each variable, our ap- proach is applicable even to large graphs. We demonstrate significant improvements compared to established approaches.
Introduction
Identifying causal relationships is profound and hard problem pervading experimental sciences such as biology (Sachs et al., 2005), medicine (Castro et al. (220). earth system sciences (Runge et al., (019), or robotics (Ahmed et al. 2020). While randomized controlled interventional studies are considered the gold standard, they are in many cases ruled out by financial or ethical concerns (Pearl [Sirtes al. 2000). In order to improve the understanding of system and help design relevant interventions, the subset of causes that have direct effect (direct causesldirect causal parents) often needs be identified based observations only. This paper assumes structural equation model (SEM) comprising (I) set of d covariates represented by random vector X Rdd whose values are determined by uniquely solvable set of d structural equations, possibly non- non-linear and possibly including cycles and confounding (2) a response variable Y E R, who is not a parent of any X and whose value is determined by a linear structural equation of the form,
Y = X0 + +U, with 0 E Rdd
(1)
Equal contributions

---

Published in Transactions on Machine Learning Research (08/2022)
where U is exogenous variable with zero mean, independent from any other exogenous variables of the SEM and (5.)) denotes the inner product. Such SEM is exemplified in Figure [ Uniquely solvability of SEMs amounts to not having self-cycles in the causal structure, but any other arbitrary non-linear cyclic structure between covariates is allowed (Bongers et al.) 2021), possibly including hidden confounders, as long as there is no hidden confounder for the response variable (this would violate the assumption of independence of U). Practically speaking, almost all causal discovery applications lie under the umbrella of simple SCMs (Bollen 1989. Sanchez-Romero al. 2019).
Besides, the assumption of not having self-cycles is usually assumed not-limiting in the literature (Lacerda et al. 2012 Rothenhäusler et al. 2015: Bongers et al. (016)).
In this paper, we investigate how to find the direct causes of Y among a high-dimensional vector of covariates X. From our formulation, a given entry of 0 should be non-zero if and only if the variable corresponding to that particular coefficient a direct causal parent (Peters et al. 2017), e.g., X1 and X2 in Figure || We restrict ourselves to the setting of linear direct causal effects Y (LDC, as specified in Equation[] and no feature descending from Y (NFD). LDC justified as an approximation when the effects of each causal feature are weak such that the possibly non- linear effects can be linearized; NFD is justified in some applications where we can exclude any influence of Y on covariate. This is, for example, the case when X are genetic factors, and Y is particular trait/phenotype. Our method, in particular, comes handy in this case due to the relatively complex non-lineat cyclic structure of these genetic factors in high- dighhdimensional regimes Yao et al. 2015
Meinshausen et al. 2016 Warrell & Gerstein 2020).
While applicable to full graph discovery rather than the simplified problem of finding causal parents, state-of-the- art methods for causal discovery often rely on strong assumptions or the availability of interventional data or have prohibitive computational costs explained in section I.] In addition and despite their strong assumptions, causal discovery methods may perform worse than simple regression baselines (Heinze-Deml al. 2018 [aazzigg 2019) Zheng et al. 2018). However, estimating 0 in high dimensional settings (i.e. # observations << d) using unregular- ized least squares regression, will lead to identifiability problems since there can be infinitely many possible choices for 0 recovered with equivalent prediction accuracy for regressing Y (Bühlmann & Van De Geer 2011). On the other hand, when using a regularized method such as Lasso, a critical issue is the bias induced by regularization (Javanmard & Montanari 2018). While regularized and unregularized Least Squares as well as debiased Lasso are applicable to the same problem, they have different statistical characteristies and we will explore their empirical performance in extensive experiments.
Double ML approaches (Chernozhukov et al. 2018a) have shown promising bias compensation results in the context of high dimen- sional observed confounding of a single variable. In the present pa- per, we use this approach to find direct causes among a large number of covariates. Our key contributions are:
Direct Causal Parents
X6
X5
X1
U
Y
We show that under the assumption that no feature of X is a child of Y, the Double ML (Chermozhukov et al.) 2018) principle can be applied in an iterative and parallel way to find the subset of direct causes with observational data.
X4
X2
0g
X
Y= 8, +0,x,
Non-Linear (with possible Cycles) Relationship
Our approach has a computational complexity requirement polynomial (fast) time in dimension d.
Our method provides asymptotic guarantees that the set can be recovered from observational data. Importantly, this re- sult neither requires linear interactions among the covari- ates, faithfulness, nor acyclic structure.
Graphical representation of Causal Feature Selection in our setting, for the case of two direct causal parents of Y, X2 and X2. Xa). such that Y = 0 + 02X₂ + U. U being an independent zero-mean noise. We propose an approach to find X, and X2 under assumptions discussed in the text. An example of this setup in the real- world is finding genes which directly cause phenotype.
Figure 1:
out of variables {X1,
Extensive experimental results demonstrate the state-of- the-art performance of our method. Our approach signif- icantly outperforms all other methods (even though under- lying data generation conditions favor them), especially in the case of non-lineat interactions between covariates, de- spite relying only on linear projection.
2

---

Published in Transactions on Machine Learning Research (08/2022)
1.1 Related work
The question of finding direct causal parents is also addressed in the
literature as mediation analysis (Baron & Kenny 1986 Hayes 2017; Shrout & Bolger 2002 Several principled
approaches have been proposed (relying. for instance, on Instrumental Variables (IVs)) (Angrist & Imbens 1995;
Angrist al. 1996 Bowden & Turkington| 1990) test for a single direct effect in the context of specific causal
graphs. Extensions of the IV-based approach generalized IVs-based approaches (Brito & Pearl 2012 Van der
Zander & Liskiewicz 2016) are the closest known result to discovering direct causal parents. However, no algorithm is provided in Brito & Pearl (2012) to identify the instrumental set. Subsequently, an algorithm is provided in Van der
Zander & Liskiewicz (2016) for discovering the instrumental set in the simple setting where all the interactions are
linear and the graph is acyclic. In contrast, our method allows non-linear cyclic interaction amongst the variables.
Several other works have also tried to address the problem of discovering causal features. The authors review work causal feature selection in Guyon & Aliferis (2000)) More recent papers on causal feature selection have appeared since (Cawley] 2008 Paul 2017: Yu et al., 2018), but none of those claims to recover all the direct causal parents asymptotically or non-asymptoticall, as we do in our case. There has been another line of works on inferring causal relationships from observational data based on conditional independences, such as the PC-algorithm, which can be used for more general causal inference purposes, at the expense of extra assumptions, such as faithfulness, allowing inference of the Markov equivalence class through testing iterative testing of d-separation statements (Mastakouri et al. 2019. Pearl 2009 Spirtes et al. 2000p. In contrast, under our L.DC, NFD and totally independent U assumptions, only small subset of d-separation relationships are relevant, and further, these assumptions imply that the observed relevant conditional independences and associated d-separation equivalent, so that full set of faithfulness assumptions about conditioning sets besides X-j are not necessary.
Another approach is to restrict the class of interactions among the covariates and the functional form of the signal-noise mixing (typically considered additive) or the distribution (e.g., non -on-aausianntt) to achieve identifiability (see (Hoyer ct al. 2009; Peters et al., 2014); this includes linear approaches like LiNGAM (Shimizu et al. 2006) and nonlinear generalizations with additive noise (Peters et al., 2011). For recent review of the empirical performance of structure learning algorithms and detailed description of causal discovery methods, we refer to Srrinze-Demil Recently, there have been several attempts at solving the problem of causal inference by exploiting the invariance of prediction under causal model given different experimental settings (Ghassami et al.) 2017 Peters et al., 2016). The computational cost to run both algorithms is exponential in the number of variables when aiming to discover the full causal graph.
al., 2018).
Our method mainly takes inspiration from Debiased/Double ML method (Chernozhukov et al., 2018a) which utilizes the concept of orthogonalization to overcome the bias introduced due to regularization. We will discuss this in detail in the next section. Considering a specific example, the Lasso suffers from the fact that the estimated coefficients are shrunk towards zero, which is undesirable (Tibshirani & Wasserman 2017). To overcome this limitation, a debiasing approach was proposed for the Lasso in several papers (Javanmard & Montanari] 2014 2018 /Zhang & Zhang 2014). However, unlike our approach, Debiased Lasso methods do not recover all the
coefficients of the parameter vector under the generic assumptions of the present work. To be more specific, (Javanmard & Montanari 2018) is built upon the Equation ① with the following differences to our setting: Noise is Gaussian; X has independent eero-mean Gaussian rows with covariance matrix Σ satisfying specific bounding conditions; sparsity conditions su = o(n/(log p) 2 and min(sq, sp) of /n/ log where so and sq are sparsity levels of the true coefficients 0 and the precision matrix Ω = E-1 of X respectively.
2 Methodology
Before describing the proposed method, we discuss our general strategy well as Double ML and Neyman orthogo- nality in the next sections, which will be helpful in building the theoretical framework for our method.
1p= d 1 is the covariate dimension.
3

---

Published in Transactions on Machine Learning Research (08/2022)
2.1 Reduction to a nonparametric estimation problem
According to Equation ① determining whether Xj is a parent of Y in our setting amounts to testing whether 0j + 0. Let x- = X \ Xj, this can be reduced to testing whether the following estimand vanishes:
Xj E [Y - B(Y x-;)) (X, B(X; xx-)))]
(2)
Indeed, U independent of X entails Y - E(Y X-)) = 0g (X, - E(X, x-,)) + U. This leads to
Xj = 0 (X, [X,, x x-)))²
OX (X, - B(x, X-1))]
(3)
Under mild assumptions, testing whether 0 * 0 thus reduces to testing whether X; * Equation shows that Xj constitutes non-parametric estimand, i.c. a model-free functional of the observed data distribution. Nonparametric estimation results (Robins et al. 2008 Van der Laan al. 2011 Chernozhukov et al. 2018a) make use of the efficient influence function of such estimand (see e.g. Hines et al. (2022)) to derive valid estimates and confidence bounds, while allowing the use of data adaptive estimation strategies, such as machine learning algorithms. The resulting strategies are known as target learning and debiased/double machine learning, and are suitable challenging settings such as ours when is high dimensional with possibly non- non-linear dependencies among components.
2.2 Double Machine Learning (Double ML)
Double ML constitutes one possible way to derive efficient nonparametric estimates. We introduce it with the partial linear regression setting introduced in Chernozhukov et al. (2018a) Example 1.1). Given a fixed set of policy vari- ables D and control variables X acting as common causes of and Y, we consider the partial regression model of Equation
Y = D8g + (X) + U, EU X, =
(4)
D mx) VV, E[V[X] 0,
where Y is the outcome variable, U, V are disturbances and go, mg Rdd R (possibly non-linear) measurable functions. An unbiased estimator of the causal effect parameter Og can be obtained via the orthogonalization approach as in Chermozbukov et al. (2018a), which is obtained via the use of the "Neyman Orthogonality Condition" described
below.
Neyman Orthogonality Condition: Let W denote the collection of all observed variables. The traditional estimator of 8g in Equation (4) can be simply obtained by finding the zero of the empirical average of a score function é such that ((W:0.)
DT (Y DØ x()).. However, the estimation of 6 is sensitive to the bias in the estimation orthogonalization approach to get an estimate for On that is more robust the bias in the estimation of nuisance parameter (mo, 9o). Assume for moment that the true nuisance parameter is no (which represents my and 9o in Equation (4) then the orthogonalized "score" function tb should satisfy the property that the Gateaux derivative operator with respect to vanishes when evaluated at the true parameter values:
of the function g. Neyman (Nemann 1979) proposed
0, (W; 0n vo) 7ol 0.
(5)
One way to build such a score, following Chernozhukov ated to maximum likelihood-like estimate. Let &(W; (f, m)) be the log likelihood function another smooth objective for which the true parameter is the unique maximizer. The true parameter then satisfies Edg((W; (fu, nu)) 0. suggesting to start with dgl(W; (0o, no)) (biased) score. In order to compensate the bias due to the nuisance parameters, we then subtract a linear function of the derivative of the likelihood with respect it, leading to the orthog- onalized score
2018a) leq. (2.7)), start from biased score associ-
* (:: 0, n)) O, (0, - (w: (0,)). 1))
where μ is determined by the constraint of Equation ( (see proof of Proposition 6 in appendix). The correspond- ing Orthogonalized Double/Debbaseed MI. estimator On solves constraint of vanishing empirical average of the
4

---

Published in Transactions on Machine Learning Research (08/2022)
orthogonalized score, based on n-iid samples (wi.e of the observed variables.
(W; f1)) = 0,
(6)
where fi is the estimator of To and &b satisfies condition in Equation 5 in Equation the orthogonalized score function vi is,
For the partially linear model discussed
(W:0, - DO -0x)(())( mXX))
with 77 (m, g). This leads an debiased estimator satisfying
8 D.D.
= - E (0X-))((.)))00 mo(()))
(8)
which relies on the "double" use of machine learning algorithm: once to learn (X;) and once to learn (
hence the name Double ML for such estimator. We can further relate this approach to the design an estimator of the non-parametric estimand of previous section.
Indeed by subtracting 0 E mo (D)
- mXx())) on both sides of eq. 8 get
1 0.- - 555 (D. mo (X,))) E(x. (Y, - (X,) (D. - -mo(X,))
Noticing that X| = [ 9o(X)
0 (X) + 9o()., the term 0oyioo (X,) G.(( in eq. (9 appears as
ML estimator of E X). such that we recognize the right hand side of Equation (9) a Double ML estimator of
E[[(Y (X [[]] which is special case of the non-parann non-parametric estimand x, defined in Equation 33.
for the setting Xj D and X = X-j- j. In practice, we directly learn an ML estimator of X| by predicting
using X, relying on the double robustness of the Xy estimands (Smucler et al. 2019). as described in section 2.5
From Double ML to Causal Discovery: The distinction between policy variables and confounding variables is not always known in advance. Fortunately, as described section 2.1] Double ML relies on estimating a conpaaamectrric estimand that does only depend on observational data and not on the causal model. This will allow us to exploit the same approach iteratively in the setting of causal discovery. To this end, we consider a set of variables X X2, Xa) which includes direct causal parents of the outcome variable Y well other variables. We also reiterate assumption that the relationship between the outcome variable and direct causal parents of the outcome variable is linear. The relationship among other variables can be cyclic and nonlinear. We now provide a general approach to scanning putative direct causes scaling "polynomiall"" in their number (see Computational Complexity paragraph in next section), based on the application of a statistical test and Double MI. estimators. We describe first the algorithm and then provide theoretical support for its performance.
2.3 Informal Search Algorithm Description
Pseudo-code for our proposed method (CORTH Features) is in Algorithm 1 The idea is to do a one-vs-rest split for each variable in turn and estimate the link between that particular variable and the outcome variable using Double ML. To do so, we decompose Equation [T] to single out a variable = X₁ as policy variable and take the remaining variables Z X-+ X as multidimensional control variables, and run Double ML. estimation assuming the partial regression model presented in Section 2.2 which now takes the form
Y D9 + (Z) +U, B 0,
(10)
D my Z) + +V. E [ = 0.
The step-wise description of our estimation algorithm goes as follows:
(a) Select of the variables X; to estimate its (hypothetical) linear causal effect 0 on Y.
5

---

Published in Transactions on Machine Learning Research (08/2022)
(b) Set all of the other variables X as the set of possible confounders.
(c) Use the Double MI. approach to estimate the parameter 0 i.e. the causal effect of X; on Y.
(d) If the variable X; is not causal parent, the distribution of the conditional covariance Xi (Proposition B is a Gaussian centered around zero. We use a simple normality test for X; to select or discard X; as one of the direct causal parents of
We iteratively repeat the procedure on each of the variables until completion. Pseudo-code for the entire procedure is given below in Algorithm [
Guaranties for this approach identify the true parents rely the assumptions stated in Section 2.5 Equations (1315). They notably allow for hidden confounders between covariates, as long as those are not direct causes of Y, not descendent of Y. On the contrary, if Y is an ancestor of any covariate, the search algorithm may fail in both directions (false positive and false negative).
Note that Equation (TO) is not necessarily a correct structural equation model to describe the true underlying causal structure. In general, for instance, when 1) actually causes Z, it is non-trivial show that the Double MI. estimation of parameter 0g will be unbiased (see Section 2.4
Algorithm Efficient Causal Orthogonal Structure Search (CORTH Features)
1: Input: response Y E RV, covariates X E IR Nxd, significance level number of partitions K.
2: Split N observations into K random partitions, 45 for k = ......& K, each having n N/ K samples.
3: for i = -,...., d do
4: for Subsample E [K] do
5: Dr x and Zx X
Fit m |/]
(Z) to DA, and fit 9 (Z to ylW
Daj m [Vb] (Zkj) for all j E Is
d( C)Ee 1/14 Dki) -1 1 Ejen,
vi (y
[1 MMM (Zkz))
Vij
9ij
[k] + 1
n LjEn (-Y, mij
(1 (Zks) [/6 (Zkg) Disgi (Zxi) m (Zkz (i
(2ks) + Y Y1AAADa)
10: lo!p + 1
n EEEEee (- v, m
m
[WA] (Zxz) Dijgij
(1 (Zkz) ) + -m Zzkz )os3 (11)
(ZRz) + +Y, (
Dij x!!h)²
11: end for
12: 0: K ΣLER 0/ Ri K EbEK x.t1 and o7 +- * Eoex(o!t1p
13: end for
14: for E [d] do
15: Gaussian normality test for Ei 2 N/ ((,, x+) with a significance level and select 7th feature if null-hypothese is
rejected.
16: end for
17: Return Decision Vector
Remarks on Algorithm
x vector which corresponds the samples chosen in the htth subsampling pro- x cedure, ||] (x|+) "..,,.".!!! x 1+1 141...... x x for any i E [d]. In general the subscript represents the estimation for the jth variable and super-script & represents the kah subsampling procedure. K represents the set obtained after sample splitting. m. m[/ are (possibly nonlinear) parametric functions fitted using (1st....,k k. 1th k 1th ..... K
subsamples.
Computational Complexity: For each subset randomly selected from the data, we fit two lasso estimators. Ac- celerated coordinate descent (Nesterov) 2012) can be applied to optimize the lasso objective. To achieve e error. O (Kmm log 1) number of iterations are required where Kima the maximum of the two condition number for both the problems and each iteration requires (nd) computation. Hence, the computational complexity of running our approach is only polynomial in d.

---

Published in Transactions on Machine Learning Research (08/2022)
2.4 Orthogonal Scores
Now we describe the execution of our algorithm for a simple graph with 3 nodes. Let us consider the following linear structural equation model as example of our general formulation:
Y : 0 X1 + 02 + Es. X2 == 412 +2, and X1 = 61.
(11)
Example 1. Consider the system of structural equation given in Equation (0)) If Et, E2 and Ea are independent uncorrelated noise terms with zero mean, Algorithm 7 will recover the coefficients 01 and f2.
A detailed proof is given in Appendix A.I A.I While the estimation of the parameter 81 is in line with the assumed partial regression model of Equation (TT). the estimation of 02 does not follow the same. However, it can be seen from the proof that 02 can also be estimated from the orthogonal score Equation (7)
We now show that this result holds for more general graph structure given in Figure 2 allowing for non non-linear cyclic interactions among features.
Proposition 2. Assume the structural causal model of Figure 2 with (possibly non-linear and confounded) assign- ments between elements of X = [X,, with X- = [2 2,1 parameterized by 7 = (71, 712).. Assume the unconfounded linear structural assignment Y = xx0 x 19++0 +U, with U zero mean random variable with finite variance of > O, independent of X. Then, the score
v (W; 0, B) = (Y Xa0 - x - rxx Xx+)),
(12)
with TXX-1 X_X - "1 folloos the Neyman orthogonality condition for the estimation of with nuisance parameters v (3, Y) which reads
E (v X&O X- X- BD) rxx xxxx,x 0.
Please refer to Appendix A.2 for the proof. Applying Equa- tion (6), this leads to the debiased estimator
Xx
Ø E. x X (X fxx X- x-x)
72
E, XiI x
712
which relies on ML estimates / and fxx_ Comparing the score in Equation (2T) with the score in Equation 7 there are two takeaways from Proposition 2 (i) the orthogonality condition remains invariant irrespective of the causal direction between X₁ and Z, and (ii) the second term in Section 2.4 re- places function m by the (unbiased) linear regression estimator for modelling all the relations; given that the relation between Z and linear, even if relationships between Z and Xş non non-linear (See Appendix B for concrete examples). Com- bining with the Double ML theoretical results (Chernozhukov al. 2018a), this suggests that regularized predictors based on Lasso or ridge regression are tools of choice for fitting func-
Zi
Z22
B
32
A generic example of identification of causal effect in the presence of causal and anti- causal interactions between the causal predictor and other putative parents, and possibly arbitrary cyclic and nonlinear assignments for all nodes except Y (see Proposition 2). We have X Z₁ U Zz.
Figure 2:
tions (m.g).
2.5 Statistical Test
We now provide a theoretically grounded statistical decision criterion for the direct causes after the model has been fitted. Consider (Y, X7, Y E R, X E Rdd, satisfying
== (0, + U. B(Y2) < co, B(U²2) < 00, R(U) E(U X) 0, and [(I co,
(13)
(14)
E
(X, B(X, |x,))) # 0, for all j E {1.....d,,
(15)
7

---

Published in Transactions on Machine Learning Research (08/2022)
where U is an exogenous variable and X-j represents all the variables except Xz. The assumptions made with the above formulation are standard in the orthogonal machine learning literature (Rotnitzky et al. 2019 Smucler et al.) 2019 Cheroozuukyy et al. 2018). They allow identifying causal parents based on estimates of conditional covariances Xj defined in Equation
Proposition 3. Let PAy - Ui E (L,...,,dd:: :0; + 03. Then under the conditions given in Equations for each j E {1.....pp
13 to
15
a) Xj = (X, - E(x, x_ x=)))² and j E PAy if and only if Xj + O.
b) We also have (with notations of Prop. Xj E [YY B(Y X x-1)) (XX_XxXx,X-)]]
The proof is given in appendix A.3 There are two main implications of the results provided in Proposition (i X is non-zero only for direct causal parents of the outcome variable, and X2 has double robustness property as shown in (Rotmitzky al., 2019; Smucler al. 2019; Chernozhukov et al., 2018). Having double robustness property means that while computing the empirical version of the X, which we denote as Xソ, one use regularized methods like ridge regression or Lasso to estimate the conditional expectation (function m). Afterward, one can perform statistical tests on top of it to decide between zero or non-zero tests. (ii) In line with the above orthogonal score results, we see that this quantity can be estimated using linear (unbiased) regression to fit the function m, although interactions between features may be non-linear. Next, we discuss the variance of our estimator so that a statistical test can be used to identify causal parents. For the sake of convenience, the case of partitions (K 2 is explained here.
Variance of Empirical Estimates of Xj= Suppose we have i..d. observations indicated by Dn {(X,, Y,), 10.,i=1.... n). nn. let T = 2, for = 2 let To = 1. For = 1.2, compute estimates of Tet (Y X-j) and Eet (X, X-j) using the data in sample T.. Following Smucler al. (2019), we can use estimates of E (Y X-y) and E (X, X-x) that are solutions of f1-regularized regression problems to obtain square root N convergence guaranties. We use Lasso as the estimator for conditional expectation in the experiments. Now, we compute the cross-fitted empirical estimates of Xj and associated empirical variances
Randomly split the data in two halves, say 1n1 and Dn2-- Take (1,..., . d). For k =
23 = Pakk
- YY (X,
X
X,
X
X
X-)) + x 2 (x, x-,) + YX, - x3) (17)
(X,
(16)
(7))² = Pakk
((-EE* (x, x
X, E
-xx)
where In denotes the empirical average over the half. Finally, let
1 2
& =
oj = +
(18)
Consistency of such estimators notably relies on sparsity assumptions for ground truth models in the asymptotic high- dimensional setting where covariate dimension p(= d - 1) and sparsity is allowed to vary with number of samples In the notations below, we drop this dependency and specialize to our case.
Definition& 4 (Approximate Linear-Sparse Sparse class (ALS)). Ground truth predictor c(X-;) belongs to the approximately sparse class whenever there exists 0" E TRP and a function ((x-s) satisfying
dZ) = 10" (()) + r(Z), where ||0*) lo < s and E[((z)*1 < (slog(p)/m).
Theorem 1 of (Smucler et al.
2019) provides general conditions under which (see also (Chernozhukov et al. 2018)).
when the estimators 12 (Y
X and t (X, X-;) are Lasso-type regularized linear regressions.
13 SS and j E {l,...,d} Assume: (i) For true parameter in (13). 1101o S with <lo(pp/n + 0; (ūi Txx_, X_ is in the ALS class with s slog(p)/n - O; (ii) All X have support bounded by the same constant; (iv) U is independent of X and has tail decaying at least as fast as an exponential random variable;
Proposition 5. Given eqs.
Extension arbitrary numher of data partitions (K 2 2) is straightforwand. Check Algorithm 1
8

---

Published in Transactions on Machine Learning Research (08/2022)
W) h has eigenvalues lower and upper bounded: (vi) - and - -X[x,XX-,lll are bounded with variance lower bounded.
Then the pair of estimators (Ri, 03) of eq. 18). using 1 regularization coefficient X 22 /log(p)/nLfo for both Lasso
estimates, satisfies XJ 3N X1
+00.
All stated bounds are uniform across n and strictly positive. Proof is provided in Appendix A.4 with the broader context on the doubly robust estimation framework. The provided conditions correspond to a concise special case of the general statements provided in (Smucler al. (2019). It is possible to relax some of these assumptions, notably to allow some misspecifications of the sparsity assumptions for the Lasso estimates. Overall, the mildness of these assumptions may explain the empirical success of this approach in the following section.
Under conditions of proposition the test that rejects Xj = 0 when 18y1 2 1.96 confidence level. The probability of rejecting the null when it is false is
o Vn
will have approximately 95%
(1R
oj VT
3j Vn
[ 1R 1.96
{8
X S lxll
1.96
+
In order to account for multiple testing, we use Bonferroni correction.
Conditional Independence Tests: Asymptotically, the conditional independence testing between Y and Xj given X-j is also a possible solution for proposed approach. Indeed, d-separation rules imply that true causes are condi- tionally dependent according to this test, while non-causes are conditionally independent (because X-j is not collider under our NFD assumption). However, conditional independence testing is challenging in hiigh-diiensiiona/nnnn-ieeee settings. Kernel-bace based conditional independence testing is computationally expensive (Zhang et al.) 2012). We used Xj in the paper because was already known from previous works (Smucler et al.) 2019 Chernozhukov et al. 2018b) that has double robustness property, which means one can use regularized methods like Lasso to estimate empirical conditional expectation from finite number of samples and the empirical estimator is still unbiased with controlled variance. Our work is related to the recent work of (Shah & Peters 2020), which proposes conditional independence test whose proofs rely heavily on (Chernozhukov et al. 201 8a). In this paper, we use for the first time such double
ML-based tests for the search problem.
3 Experiments
3.1 Experimental Setup
To showcase performance of our algorithm, we conducted two sets of experiments: i) Comparison with causal structure learning methods (Casual and Markov Blanket discovery) using data consisted of DAGs with high num- ber observations-to-nnniii
of variables ratio (n > d) which is applicable to causal structure learning methods. Markov Blanket discovery methods are included since under NFD, faithfulness, and no-hidden-coofooudee assump- tions, Markov Blanket of the target variable corresponds to the direct parents. Note that, faithfulness and no-hidden- confounders assumptions are not necessary for method. These experiments are discussed in details in Section B.1.1 ii) Comparison with inference by regression methods using data consisted of DAGs with high number of observations- to-number of variables ratio (n & d and n < d) to illustrate performance in high-dimensional regimes. This part is explained thoroughly in Section 3.1.2
3.1.1 Causal Structure Learning
For every combination of number of nodes (#nodes), connectivity (p,), noise level (a²), number of observations (n). and linerr probability (p,) (see C.II. 100 examples (DAGs) are generated and stored as esv files (altogether 72.000 DAGs are simulated, comprising a dataset of overall > 10GB). For each DAG, n samples are generated. We provide more details about the parameters (#nodes, Ps. Pn. and n)) and data generation process in Appendix I.1. C.I.I For future benchmarking, the generated files with the code will be made available later.
The baselines we compare our method against are categorized in two groups which are suitable for observational data: i) Causal Structure Learning methods: LINGAM (Shimizu al. 2006), order independent PC (Colombo)

---

Published in Transactions on Machine Learning Research (08/2022)
noise level
1,00
1.00
0.9
0.95
0.95
0.8
0.90
0.90
Accurney
Accuracy
0.85
Accuracy
0.85
3.7
0.60
0.00
0.75
0,70
0.70
0.65
0.4
0.65
0.60
nonlinearity
0.1
0.3
0.5
0.1
0.3
0.5
0.1
0.3
0.5
Spursity
Sparsity
Spansity
(a) Pr. 0, 0
0.1
(b)
0, 0
0,5
(c)
0, a
1.0
1.00
1.00
1.0
0.95
0.95
0.9
2.90
Accuracy
0.85
Assuracy
0.90
Accurac
0.8
0.85
0.89
0.7
0.80
0.75
0.6
0.70
0.75
0.5
0.65
0.70
0.1
0.3
0.5
0.1
0.3
0.5
0.1
0.3
0,5
Sparsity
Sparsity
Sparsity
(d) Po
0.5, o
0.1
(c) Pm = 0.5, o² 0.5
(f) Ps. = 0.5, o² 1.0
Figure 3: Overall performance for single random DAG with 100 simulations for each setting, having 20 nodes and 500 observations.
& Maathuis 2014). rankPC, MMHC (Tsamardinos et al., 2006). GES (Chickering 2003). y. rankGES, ARGES (adap-
tively restricted GES (Nandy et al. 2016)). rankARGES, FCI+ (Claassen et al. 2013), PCI (Shah & Peters 2020)
and Lasse [ (Tibshirani 1996). ii) Markov Blanket discovery methods: Grow-Shrink (GS (Margaritis & Thrun
1999)), Incremental Association Markov Blanket (IAMB (Tsamardinos et al. 2003b)) Max-Min Parents & Children (MMPC (Tsamardinos et al.) (003a)) FastlAMB (Yaramakala & Margaritis, 2005). and LAMB with FDR Correc- tion (Pena, 2008). The "CompareCausalNewwoks [ and "bnlearn: Bayesian Network Structure Learning, Parameter Learning and Inference [ Packages are used to run most of the baselines methods. We use 10-fold cross-validation to choose the parameters of all approaches. As direction of the possible causes in the defined setting is determined, the non-directional edges inferred by some baselines, c.g., PC are evaluated as direct causes of the target variable.
3.1.2 Inference by Regression
Similar to the previous section, for every combination of parameters, 50 examples are generated and stored, which means 15000 DAGs overall. Details are provided in Appendix [.1.2] We compare our algorithm to methods for
inference in regression models: Standard Regression, Lasso with exact post-selection inference (Lee et al.
2016).
Debiased Lasso @Jvanmard et al. 2015), Forward Stepwise Regression for active variables (Loftus & Taylor
2014
Tibshirani et al. 2016), Forward Stepwise Regression for all variables (Loftus & Taylor, 2014 Tibshirani et al.
2016)
LARS for active variables (Efron et al. 2004] Tibshirani et al. 2016), and LARS for all variables (Efron et al.
2004
Tibshirani et al. 2016). "selectivelnference: Tools for Post-Selection Inference" " R Package | leveraged to run most of these basclines. We used cross-validation to choose hyperparameters and confidence level for hypothesis testing
considered is 90%.
Regression Technique and method because the variance bound for Xj with Lasso type estimator of conditional expectation is provided in equa-
We use Lasso as the estimator of conditional expectation for our
NNone-zero Nonne-eenoceefficiinttts
coefficients are reported.
Netwarks /ites///ddxxee.. html
10

---

Published in Transactions on Machine Learning Research (08/2022)
noise level
Recsarert
Raperent
0 (a) Pn 0.3, 0² 0.3
(b) pn 0.3, 0² 0.5
Pm 0.3, 0² 1.0
Figure Distribution of the estimated 0 values for the true and false causal parents in 100 simulations of the graph with 20 nodes, 20000 observations and 0.3 connectivity. The vertical lines indicate the ground truth values for the linear coefficients corresponding to causal parents.
tion 17 Further, using more splits than 2 splits in the experiment relatively increases the performance of parameter estimation. See Figure for parameter estimations.
Evaluation: Recall, Fall-out, Critical Success Index, Accuracy, Fl Score, and Matthews correlation coefficient (Matthews 1975) are considered as metrics for the evaluation. These metrics are described in Appendix
3.2
Results
3.2.1
Causal Structure Learning
Results aggregated by the number of observations (corresponding to 24000 simulations per entry in the table) are illustrated in Table ! Our method performs better than the competing baselines in terms of accuracy and Fl score, especially for more connected structures, despite data being generated according DAG causal structures, which, dissimilar our method, an essential condition for them. To provide a visual comparison, we plot the accuracy of all methods w.r.t. the connectivity parameter (px) in Figure 3 for different values of pn. and o² on 1800 samples.
It can be observed that the accuracies of the competing baselines significantly drop with increasing noise level and nonlinearity, while our method is more robust to them. We also extensively compare all the metrics (Recall, Fall- out, Critical Success Index, Accuracy, F1 Score, and Matthews correlation coefficient) for all the methods Ap- pendix C.3.I According to these metrics, our approach performs better than baselines in most cases regardless of the set of parameters used for generating data. Our method shows in particular stability in performance w.r.t. the number of nodes (Table C.3), partially non- non-linear relationships (Table C.4. connectivity (Table(..5) number of observations (Table C.7, and noise level (Table C.6). We also show the plot of parameter estimation for direct causal parents vs. non-causal parents in Figure In the plots and tables, we denote approach CORTH Features.
3.2.2 Inference by Regression
Analogous to previous part, results are aggregated by nonlinear probability (corresponding to 3750 simulations per entry the table), number of observations (3000 simulations per entry the table), connectivity (5000 simulations per entry in the table) and beta distribution parameters are provided Tables Tables|C.8 to C.IT Based on these results, our method suggests more robustness w.r.t. the set of parameters used for generating data and relatively better performance compared to other methods.
3.3 Scaling Causal Inference to Large Graphs
Figure 5 shows the runtime the method in sees as a function of the graph's size. Notice that the runtime of our algorithm in the log-log plot roughly linear, supporting our above statement about the computational time being polynomial d. As used 5000 observations, additional overhead comes from cross-validation.
7 Please refer Appendix C.3.1 for thoruugh tables all parameters.

---

Published in Transactions on Machine Learning Research (08/2022)
Mithed
10
OORTH Features
Number of Observations
Method
100
500
1000
ACC
F1
ACC
F1
ACC
FI
GES
0.80
0.59
0.81
0.65
0.81
0.67
10a
rankGES
Runtime (sec)
ARGES
0.78
0.49
0.80
0.58
0.80
0.59
RankARGES
0.78
0.47
0.79
0.57
0.80
0.58
FCI+
0.84
0.67
0.86
0.75
0.87
0.78
LINGAM
0.84
0.65
0.91
0.74
0.94
0.88
10
rankPC
0.82
0.62
0.85
0.71
0.85
0.73
MMPC
0.77
0.37
0.53
MMHC
0.80
0.56
0.82
0.62
0.83
0.64
GS
0.79
0.43
0.84
0.59
0.86
0.62
10a
JAMB
0.74
Fast-IAMB
0.80
0.46
0.84
0.59
0.86
0.62
TAMB
0.37
PCT
0.83
0.59
0.91
0.85
0.93
0.89
5
20
50
100
200
500
1000
Lasso
0.87
0.81
0.89
0.85
0.89
CORTH Festures
0.88
0.78
0.93
0,91
0.94
0.92
Number of nodes
Table Performance across all the settings for dif- ferent number of observations (100, 500 and 1000). Each single entry in the table is averaged over 24000 simulations. Our method is almost state of the art in every case.
Figure 5: Runtime function of the number of variables for 10 simulations per number of nodes. In these simulations connectivity, number of observa- tions, nonlinaer prob., and noise level are set to 0.3, 5000, 0, and respectively.
3.4 Real-World Data
We also apply our algorithm to a recent COVID-19 Dataset (Einstein 2020) where the task is predict COVID-19 cases (confirmed using RT-PCR) amongst suspected ones. For existing and extensive analysis of the dataset with predictive methods, we refer to Schwab et al. (020). We apply our algorithm to discover the features which directly cause the diagnosed infection. We found that the following were the most common causes across different runs of our approach: Patient age quantile, Arterial Lactic Acid, Promyelocytes, and Base excess venous blood gas analysis. Lacking medical ground truth, we report these not as corroboration our approach but rather as a potential contribu- tion to causal discovery in this challenging problem. It is encouraging that some of these variables are consistent with other studies Schwab et (2020). Details on data preprocessing and more results are available Appendix D
4 Discussion
A recent empirical evaluation of different causal discovery methods highlighted the desirability of more efficient search algorithms (Heinze-Deml et al. 2018). the present work, we provide identifiability results for the set of direct causal parents, including the case of partially nonlinear eyclic models, as well as a highly efficient algorithm that scales well w.r.t. the number of variables and exhibits state-of-the-art performance across extensive experiments. Our approach builds the Double ML method Chemoozukov et al. 201 and properties of conditional covariance estimands, which can leverage Lasso predictors to obtain consistent estimators in high dimensional settings (Suucer et al. (2019). Theoretical properties of such approaches (Smucler et al. 2019 Section 2), such as rate double robustness, which accommodates non-strict (approximate) sparsity assumptions, as well as model double robustness, which allows some degree of model misspecification, may explain the empirically observed performance of our approach. Whilst not amounting full causal graph discovery, identification of causal parents is of major interest in real-worke applications, e.g., when assaying the causal influence of genes on the phenotype. A natural direction worth exploring is extend this approach for discovering direct causal parents in the case when nonlinear relationships exist between the output variable and its direct causal parents.
Acknowledgements
This work was supported by the German Federal Ministry of Education and Research (BMBF): Tübingen AI Center, FKZ: 011S 10338....B.....
12

---

Published in Transactions on Machine Learning Research (08/2022)
References
Ossama Ahmed, Frederik Träuble, Anirudh Goyal, Alexander Neitz, Yoshua Bengio, Bernhard Schölkopf, Manuel Wüthrich, and Stefan Bauer. Causalworld: A robotic manipulation benchmark for causal structure and transfer learning. arXiv preprint arXiv:2010.04296, 2020.
Joshua D Angrist and Guido W Imbens. Identification and estimation of local average treatment effects. Technical report, National Bureau of Economic Research, 1995.
Joshua D Angrist, Guido W Imbens, and Donald B Rubin. Identification of causal effects using instrumental variables. Journal of the American statistical Association, 91(434):444 455, 1996.
Reuben M Baron and David A Kenny. The moderator-mediator variable distinction in social psychological research: Conceptual, strategic, and statistical considerations. Journal of personality and social psychology, 51(6). 1986.
Kenneth Bollen. Structural equations with latent variables, volume 210. John Wiley & Sons, 1989.
Stephan Bongers, Jonas Peters, Bernhard Schölkopf, and Joris M Mooij. Theoretical aspects of cyclic structural causal models. arXiv preprint 1rXiv: 1.06221, 2016.
Stephan Bongers, Patrick Forré, Jonas Peters, and Joris M Mooij. Foundations of structural causal models with cycles and latent variables. The Annals of Statistics, 49(5):2885-29 2915, 2021.
Roger J Bowden and Darrell A Turkington. Instrumental variables, volume 8. Cambridge university press, 1990.
Carlos Brito and Judea Pearl. Generalized instrumental variables. arXiv preprint arXiv: 2012.
Peter Bühlmann and Sara Van De Geer. Statistics for high-dimensional data: methods, theory and applications.
Springer Science & Business Media, 2011.
Daniel Castro, Ian Walker, and Ben Glocker. Causality matters in medical imaging. Nature Communications, (1):1-10, 2020.
Gavin Cawley. Causal & non-causal feature selection for ridge regression. In Causation and Prediction Challenge. 2008.
Victor Chernozhukov, Denis Chetverikov, Mert Demirer, Esther Duflo, Christian Hansen, Whitney Newey, and James Robins. Double/debiased machine learning for treatment and structural parameters, 2018a.
Victor Chernozhukov, Whitney Newey, and James Robins. Double/de-biased machine learning using regularized riesz representers. arXiv preprint arXiv:1802.08667 2018b.
Victor Chernozhukov, Whitny K. Newey, and Rahul Singh. Learning L2 Continuous Regression Functionals via Regularized Riesz Representers. arXiv preprint 159.0248,,,, Sep 2018.
David Maxwell Chickering. Optimal structure identification with greedy search. J. Mach. Learn. Res., ((nul)): 507-554, March 2003. ISSN 1532-4435.
doi: 1011621153243032117771 URL
10.1162/153244505333
Tom Claassen, Joris M. Mooij, and Tom Heskes. Learning sparse causal models is not np-hard. In Proceedings of the Twenty-Ninth Conference on Uncertainty in Artificial Intelligence, UAr' 13, pp. 172-181, Arlington, Virginia, USA, 2013. AUAI Press.
Diego Colombo and Marloes H. Maathuis. Order-independent constraint-based causal structure learning. Jour- nal of Machine Learning Research, 15(116):3221-3963 2014. URL http:://mm.r..oo/pppprsss.
co Tombol 1 4a. html
Bradley Efron, Trevor Hastie, lain Johnstone, and Robert Tibshirani. Least angle regression. The Annals of statistics, 32(2).407-499, 2004.
Hospital Israelita Albert Einstein.
Diagnosis of COVID-19 and its clinical spectrum, 2020.
13

---

Published in Transactions on Machine Learning Research (08/2022)
Patrick Forré and Joris M Mooij. Markov properties for graphical models with cycles and latent variables. arXiv preprint arXiv: 1710008777, 2017.
AmirEmad Ghassami, Saber Salehkaleybar, Negar Kiyavash, and Kun Zhang. Learning causal structures using re-
gression invariance. In Advances in Neural Information Processing Systems, pp- 301 1-3021, 2017.
Isabelle Guyon and Constantin Aliferis. Causal feature selection. In Computational methods of feature selection. Chapman and Hall/CRC, 2007.
Andrew F Hayes. Introduction to mediation, moderation, and conditional process analysis: A regression-based ap- proach. Guilford publications, 2017.
Christina Heinze-Deml, Marloes H Maathuis, and Nicolai Meinshausen. Causal structure learning. Annual Review of Statisties and Its Application, 5:371-399 19111191122 2018.
Oliver Hines, Oliver Dukes, Karla Diaz-Ordaz, and Stijn Vansteelandt. Demystifying statistical learning based on efficient influence functions. The American Statistician, pp. 1-13, 2022.
Patrik O Hoyer, Dominik Janzing, Joris M Mooij, Jonas Peters, and Bernhard Schölkopf. Nonlinear causal discovery with additive noise models. In Advances in neural information processing systems, pp. 689-696, 2009.
Dominik Janzing. Causal regularization. In Advances in Neural Information Processing Systems, 2019.
A Javanmard et al. De-biasing the lasso: Optimal sample size for gaussian designs. arxiv, 2015.
Adel Javanmard and Andrea Montanari. Confidence intervals and hypothesis testing for high dith-eiimnsional regression. The Journal of Machine Learning Research, 15(1), 2014.
Adel Javanmard and Andreas Montanari. Debiasing the lasso: Optimal sample size for gaussian designs. The Annals of Statistics, 46(6A), 2018.
Gustavo Lacerda, Peter L. Spirtes, Joseph Ramsey, and Patrik O Hoyer. Discovering eyclic causal models by indepen- dent components analysis. arXiv preprint 1rX.v:1206.3273, 2012.
Jason D Lee, Dennis L Sun, Yuckai Sun, and Jonathan E Taylor. Exact post-selection inference, with application to the lasso. The Annals of Statistics, 44(3):907-9277 2016.
Joshua R Loftus and Jonathan E Taylor. A significance test for forward stepwise model selection. arXiv preprint arXiv:1405.399 2014.
Dimitris Margaritis and Sebastian Thrun Bayesian network induction via local neighborhoods. Advances in neural information processing systems, 12, 1999.
Atalanti Mastakouri, Bernhard Schölkopf, and Dominik Janzing. Selecting causal brain features with single condi- tional independence test per feature. In Advances in Neural Information Processing Systems 32, pp. 12532-12543, 2019.
Brian W Matthews. Comparison of the predicted and observed secondary structure of t4 phage lysozyme. Biochimica et Biophysica Acta Protenn Structure, 405(2))442-45)) 1975.
Nicolai Meinshausen, Alain Hauser, Joris M Mooij, Jonas Peters, Philip Versteeg, and Peter Bühlmann. Methods for causal inference from gene perturbation experiments and validation. Proceedings of the National Academy of Sciences, 113(2)):361--7663 2016.
Preetam Nandy, Alain Hauser, and Marloes Maathuis. High- High-dimensional consistency in acore-based and hybrid struc- ture learning. Annals of Statistics, 46, 03 2016. doi: 10.1214/17AOS11655
Yu Nesterov. Efficiency of coordinate descent methods on huge-scale optimization problems. SIAM Journal on Optimization, 22(2):341-36,, 2012.
Jerzy Neyman. C (a) tests and their Sankhyā: The Indian Journal of Statistics, Series A, 1979.
14

---

Published in Transactions on Machine Learning Research (08/2022)
Michael Paul. Feature selection causal inference: Experiments with text classification. In Proceedings of the 21st Conference on Computational Natural Language Learning (CoN)J.2017), pp. 163-172, 2017.
Judea Pearl. Causality. Cambridge university press, 2009.
Jose Pena. Learning gaussian graphical models of gene networks with false discovery rate control. In European con- ference on evolutionary computation, machine learning and data mining in bioinformatics, pp. 165-176. Springer. 2008.
Jonas Peters, Joris Mooij, Dominik Janzing, and Bernhard Schölkopf. Identifiability of causal graphs using functional models. In Proceedings of the 27th Annual Conference on Uncertainty in Artificial Intelligence ( UAI), pp. 589-598, 2011.
Jonas Peters, Joris M Mooij, Dominik Janzing, and Bernhard Schölkopf. Causal discovery with continuous additive noise models. The Journal of Machine Learning Research, 15(1):2009-2053, 2014.
Jonas Peters, Peter Bühlmann, and Nicolai Meinshausen. Causal inference by using invariant prediction: identification and confidence intervals. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 78(5):947- 1012, 2016.
Jonas Peters, Dominik Janzing, and Bernhard Schölkopf. Elements of causal inference: foundations and learning algorithms. MIT press, 2017.
James Robins, Lingling Li, Eric Tchetgen, Aad van der Vaart, et al. Higher order influence functions and minimax estimation of nonlinear functionals. Probability and statistics: essays in honor of David A. Freedman, 2:335-421, 2008.
Dominik Rothenhäusler, Christina Heinze, Jonas Peters, and Nicolai Meinshausen. Backshift: Learning causal cyclic graphs from unknown shift interventions. Advances in Neural Information Processing Systems, 28, 2015.
Andrea Rotnitzky, Ezequiel Smucler, and James M. Robins. Characterization of parameters with mixed bias property. arXiv preprint wrXivv:904..3/7 2019.
Jakob Runge, Sebastian Bathiany, Erik Bollt, Gustau Camps-Valls, Dim Coumou, Ethan Deyle, Clark Glymour, Mar- lene Kretschmer, Miguel D Mahecha, Jordi Muñoz-Marí, al. Inferring causation from time series in earth system sciences. Nature communications. 10(1):1-13, 2019.
Karen Sachs, Omar Perez, Dana Pe'er, Douglas A Lauffenburger, and Garry P Nolan. Causal protein-signaling net- works derived from multiparameter single-cell data. Science, (08(57)) 2005.
Ruben Sanchez-Romero, Joseph D Ramsey, Kun Zhang, Madelyn Glymour, Biwei Huang, and Clark Glymour. Estimating feedforward and feedback effective connections from fmri time series: Assessments of statistical meth- ods. Network Neuroscience, 3(2):274-306, 2019.
Patrick Schwab, August DuMont Schütte, Benedikt Dietz, and Stefan Bauer. predcovid-19) systematic study of clinical predictive models for coronavirus disease 2019. arXiv preprint arXiv:2005. 08302, 2020.
Rajen D Shah and Jonas Peters. The hardness of conditional independence testing and the generalised covariance measure. The Annals of Statistics, 48(3):1514-158,, 2020.
Shohei Shimizu, Patrik O Hoyer, Aapo Hyvärinen, and Antti Kerminen. A linear non- non-gaussian acyclic model for causal discovery. Journal of Machine Learning Research, 7(0c)):0003-2030, 2006.
Patrick E Shrout and Niall Bolger. Mediation in experimental and nonexperimental studies: new procedures and recommendations. Psychological methods, 7(4), 2002.
Ezequiel Smucler, Andrea Rotnitzky, and James M Robins. A unifying approach for doubly-robust & regularized estimation of causal contrasts. Technical report, 2019.
Peter Spirtes, Clark N Glymour, Richard Scheines, and David Heckerman. Causation, prediction, and search. MIT press, 2000.
15

---

Published in Transactions on Machine Learning Research (08/2022)
Robert Tibshirani. Regression shrinkage and others. Journal of the Royal Statistical Society: Series B (Methodologi- (8())26722 1996.
Ryan Tibshirani and Larry Wasserman. Sparsity, the lasso, and friends, 2017.
Ryan Tibshirani, Jonathan Taylor, Richard Lockhart, and Robert Tibshirani. Exact post-selection inference for sequential regression procedures. Journal of the American Statistical Association, 111(514);600-622)) 2016.
Ioannis Tsamardinos, Constantin F Aliferis, and Alexander Statnikov. Time and sample efficient discovery of markov blankets and direct causal relations. In Proceedings of the ninth ACM SIGKDD international conference on Knowl- edge discovery and data mining, pp. 673-678, 2003a.
Ioannis Tsamardinos, Constantin F Aliferis, Alexander Statnikov, and Er Statnikov. Algorithms for large scale markov blanket discovery. In FLAIRS conference, volume 2, pp. 376-380. St. Augustine, FL, 2003b.
Ioannis Tsamardinos, Laura Brown, and Constantin Aliferis. The max-min hill-climbing bayesian network structure learning algorithm. Machine Learning, 65:31-78, 10 2006. doi: (100000/11999000---660
Mark J Van der Laan, Sherri Rose, et al. Targeted learning: causal inference for observational and experimental data, volume 10. Springer, 2011.
Benito Van der Zander and Maciej Liskiewiez. On searching for generalized instrumental variables. In AISTATS, pp. 1214-1222, 2016.
Jonathan Warrell and Mark Gerstein. Cyclic and multilevel causation evolutionary processes. Biology & Philosophy, 35(5):1-36, 2020.
Shun Yao, Shinjae Yoo. and Dantong Yu. Prior knowledge driven granger causality analysis on gene regulatory network discovery. BMC bioinformatics, 16(1):1-18, 2015.
Sandeep Yaramakala and Dimitris Margaritis. Speculative markov blanket discovery for optimal feature selection. In Fifth IEEE International Conference on Data Mining (ICDM 005), pp. 4-pp. IEEE, 2005.
Kui Yu, Lin Liu, and Jiuyong Li. unified view of causal and non-causal feature selection. arXiv preprint arXiv: aeee:8002.08540 2018.
Cun-Hui Zhang and Stephanie S Zhang. Confidence intervals for low dimensional parameters in high dimensional
linear models. Journal of the Royal Statistical Society: Series (Statistical Methodology). 76(1), 2014.
Kun Zhang, Jonas Peters, Dominik Janzing, and Bernhard Schölkopf. Kernel-based conditional independence test and application in causal discovery. arXiv preprint arXiv:1202.3773. 2012.
Xun Zheng, Bryon Aragam, Pradeep K Ravikumar, and Eric P Xing. Dags with tears: Continuous optimization for structure learning. In Advances in Neural Information Processing Systems, pp. 9472-992999 9483, 2018.
16

---

Published in Transactions on Machine Learning Research (08/2022)
A Causal Discovery via Orthogonalization
A.1 Example 1
012
X1
X2
81
Bg
Figure A.1: An example with linear structural equations.
Proof of Example 7. Let us start from the easier case first (See Figure Let us first try to estimate the coefficient of interaction between X2 and Y but it is also very clear that the estimation of 02 will be unbiased as the given setting precisely match with the double machine learning setting. However, will see in this example that given the population, 0, can be approximated as well. Let us write down the structural equation model first:
Y = 02 X1 + 02X2 + 63
X₂ 01 X1 +
(19)
X1 E1
From the set of equations we have:
-1 Xi d1
-1 X2
Let also denote E[=i]] = of and E|s²] = o... Hence, E| = of. [[] = a12 and B[X3] = [12 +
E[E2X2] of + o².. Let first try to find the regression co-efficient fitting X₃ on Y.
Y 02 m-
Hence, 0₂ E] [[X3]
if 1 is independent of X2-
E[
X0,, + 02X2 + ea)]
of
8₂ E[X3]
02 01a12
(20)
E E[X²]
o2 + a
2'
Similarly, if we fit X2 X, then
X, = a X2 + 1/2
then 0
2[X,Xg E
However [2] can also be written as following:
X|X|X]]= = a - a a 22]
Hence,
(1
of
af of + aizo?
1
0
of + ai
o
Residual V = X, a Xz. Hence we can have
of
x(LX,)= = E =
=
Ele?]
R(e)
17

---

Published in Transactions on Machine Learning Research (08/2022)
We now calculate,
V(Y - 02X2)
- EE (Xi - a X2 - 02X2)
E (X1 - X2) 22) 01 + E3)
(02 10-0222 02)01201 + 0101 2 02) a}z0 0i 010120]
010702
0 + क
Last equation was written after a step of minor calculation. Since the estimator is
0 [E(VX)) "' VCY 02X2)] 01.
A.2 Influence of the interactions between parents
In this section, we use generic example shown in Figure which we show again in Figure A.2 to illustrate the role of interactions between the covariates the proposed causal discovery algorithm.
Xx
7/2
712
0
Zi
Zz
82
Y
Figure A.2: A generic example of identification of causal effect 0 in the presence of causal and anti-causal interac- tions between the causal predictor and other putative parents, and possibly arbitrary cyclic and nonlinear assignments for all nodes except Y (see Proposition We have X- Z1 U Zz.
The estimator discussed can simply be derived from the Neyman orthogonality condition. We now provide the below the proof for Proposition 2 For the sake of completeness, we also rewrite the statement of the proposition again.
Proposition (Restatement of Proposition 2). Assume the structural causal model of Fig, A.2 with (possibly non- linear and confounded) assignments between elements of X [X,, X with X_ [Z ZTTT, parameterized by 7 (Y,, 7/, 712). Assume the unconfounded linear structural assignment mean random variable with finite variance of 0, independent of X. Then, the score
X₂0 + X T U, with zero
(W; (Y -x00 XI,,8) X-1),
(21)
with Txx_1 = E[ T X - followss the Neyman orthogonality condition for the estimation of with nuisance parameters 17 (B,7) which reads
E (r - X&O - Xxxxxx Xxx] 0.
Proof of Proposition 2
We first use a likelihood based approach under an additional Gaussianity assumption before addressing the general setting.
18

---

Published in Transactions on Machine Learning Research (08/2022)
Special case Let assume U is Gaussian. Using the global Markov factorization for simple SCM Forré & Mooij
2017 Bongers et al. 202T).
P(W; 8,
due to linearity and gaussianity of the assignment of Y, obtain negative log likelihood of the form (up to additive constants)
(W;0,n)
2 - - X_x T - X - 3)) + J(X X X-ii7).
2
where f stands for the negative log likelihood of the second factor and 77 | Y is the nuisance parameter vector. Following the principle of the approach described in main text, use Chernozhukov et al. 2018a) [Eq. (2.7))
to define the Neyman orthogonalized score, leading to:
1 (Y -- Xa0 X- T
(W;0, m) Oel (0: - vd (W; (0,
o
+Y - X,0 -x- T -)
The quantity μ should be chosen to satisfy Neyman orthogonality of Equation which leads tc
Og (W; 0, n) 0, IE, (0, -,a,,- EQ, (0. n)) = 0
leading to the expression of given in Chernozhukov et al.
2018a) (eq. (2.8)):
J 9.7
with
[[0,,(WW,W = oy
2 2 [
0 02 E(0,S(X x mm))
Jnn - 0,
and
Jo.n H, E (&b n)) = o.
E X 0 0].
resulting in
E [x E [ 0 [ 0]]
Reintroducing μ in the expression of the correction term leads to
-x0 -x0- TB)x
8)X-
(O,,e(W((Omm))
x_[(X,,_xi1))
leads to
wa,e(W; (0. =
= (Y
X x X
which leads to the final expression of the orthogonalized score :h:
do (W; (0, n)) (.) (0, - og 1 - X TxX_,
Since is a positive constant of the problem, for this estimation problem, leading to the final expression
remove it while still having appropriate orthogonalized score
v WW 0, n) = (Y - X,0 - X T (X - Txx_, X-R)
"The necessary condition this stutement be uniquely solvability which equivalent having self-cycles in the causal structure. tte transpose 2,,TT should understood organizing columnwise the partial derivatives with respect each components of
19

---

Published in Transactions on Machine Learning Research (08/2022)
General case
We now drop the Gaussianity assumption. Interestingly, the least square objective
L(W: n) (Y - - X-+
T B) -
is still an appropriate (biased) score for the estimation problem, in which the nuisance parameters 7 do not play a role. Indeed, starting from the structural equation of Y
= X x
we get
Y X = B| 0 + E x
3,
which can be rewritten as
L(W (Oo. 77)) 0.
As consequence, the true parameter is the unique minimizer of the loss, identified by this vanishing paratial derivative.
We can thus apply the previous orthogonalization procedure in the same way as we did for the likelihood, defining
(W; 0, m) dwL(W; (0, (W; (0, 17)
and leading to the exact same orthogonalized score (up irrelevant sign change).
A.3 Proof of proposition
Proof. From Equation
13
B(Y x = B((0, X- + E(U (x-))
x--)
= -0-,, x_,) + X_,) = (0, - 1x-))
- -X,
- B(X, Xx-))).
Thus
Xj BIUU + +0, - B(X, x-)))) (X, - B(X, IX-;)] [X E(X, x-j)) 0,E [[X, B(X, X-))
= 1 [(X, TE(X, 1-,)
Since E (X, B(X, X > 0. E PAy if and only if + 0. proving a)-b). For c), we rely on the properties
of the Hilbert space of square integrable RV's L (?)), equiped with the scalare product (X, Y) = E[XY]. We rewrite
Xj E [Y - IE Xx1))) (X, X-,x,, [Y IE(Y X-,)) (rxx B(X, X-j))
Under our assumptions, B(Y is the orthogonal projection of Y on the subspace of G-measurable square integrable RV's L2( so Y [{Y|X_,) is orthogonal to any elements of L 9). Noticing that ( E(X, X-;)) is element of I. (22, the second right-hand side term of the above equation van- ishes and we get the result.
A.4 Background on the Bilinear Influence Function (BIF) class and proof of Proposition
We summarize the results in Smucler et al. (2019) that justify the normal convergence the Lasso-type estimates that we relying on in Section 2.5 to perform statistical hypothesis tests. This paper investigates the properties of f-- regularised machine learning estimators for particular family of non nor-maraittric estimands, called Bilinear Influence Function (BIF) functionals.
20

---

Published in Transactions on Machine Learning Research (08/2022)
The estimand is denoted X()) where denotes the model. The bias of the estimator due to the unperfect estimate f of can be quantified using the influence function framework, leading to the following Taylor expansion in the neighborhood of the true model
X = x(f) + P.X)
where Xg is the influence function, which is Pn is the empirical distribution of n. iid samples of the observation distribution.
random variable under distribution P. of the true model, and
Estimands belonging to the BIF class are characterized by an influence function of the following form, for observed random variables O including a vector Z,
()) Sasa Sasa (Z)b(Z) + m. (O,a) + m, + -((n)
where a and are L2 and ma and vi,, are linear in the second argument (Smucler et al. 2019 Definition 1),
Our estimand of Equation (3) belongs this BIF class as an expected conditional covariance (Smucler et al. 2019 Example 5), written in the notation of that paper as Lin.L, E and Lin. V.
B((Y - E(Y - [[D]Z)))
For this case, Hines al. (2022) (among others) provide the influence function:
x² = - ((( E() x(n)
which entails the BIF function parameters Sab 1, a(Z) EYY)Z), b(Z) E(D]Z), ma(0,a) Da. b) = Yb,S. = DY.
An estimation procedure of BIF functionals based on Lasso-type estimators of conditional expectations is described in Smucler et al. (2019 Section 3.1), and corresponds to ours described section 2.5 In short, estimates are chosen among a family models of the form (a ((f. (())) and (o0(,,, (())), with parameters A., c E {a, b). ,that solve the & regularized problem of the general form
0_ arg min Pn. [Qa(0, o, w)] A|10|11 BER
(22)
for e E b] with I,, n-sample empirical average, the regularization parameter and objective function
Qc(O: o; w) (Z ((.)))) + (0, me(w 6))
In our case, we can use linear models which entail an identity link function tp = id and anti-derivative xb = x²² /1. Moreover the weight can be set to = Smucler et al. (2019), leading to
Yo)) (0, (a(0;o:w) (10, ((0.)((2))2))))))))))))))))))))) (0,
(Y (0, (0,6(2))) Y²
and
(o () = (19, (()))())))² (o()) = (D (0, 2 -
which thus both reduce to least square linear regression objectives, and thus turn eq. 22 into the classical Lasso objective.
Asymptotic properties are established for the true a()) and/or b(2) belonging sequences of models associated to each sample size n with parameter dimension p and sparsity s, quantified as an upper bound on the number of non-zero coefficients of vector, i.e. its Eg norm [|-I|o- such that
8 log(p)/n 7 0.
In particular, (Smucler et al. 2019 Section 4) use the approximately generalized linear-sparse class, for some j. such that there exists 8" R2" and function r(Z) satisfying (dependence on n is dropped ease notation)
c(Z) = =((0", ((Z))) + +r())
where |10 S and E[r(Z)?) < (sslog(p)/n) In our main text Definition we specialize this class to the linear (Φ id) and
21

---

Published in Transactions on Machine Learning Research (08/2022)
Sketch of the proof of Proposition Statement (3) of Theorem in Smucler al. (2019) provides necessary condi- tions the estimators' error to converges to zero mean distribution with accurate empirical estimate of the the variance. These conditions are gathered in 3 subsets of assumptions called Lin.L, Lin.E and Lin. V (see conditions 1-3 in (Smucler et al., 2019, Section 5.1.1) together with statements of sufficient conditions to satisfy them. Based on these statements, we go through justifications for each assumptions.
Lin.L Our assumption (i) entails that this as Lin.L. trivially satisfied for E[y|X->] (see (Smucler et al. 2019 Section 4, Example 9)). Moreover, Lin.L.I is satisfied for [[[ explicitly by our assumption (ii).
For Lin.L., we use the statement of Smucler et al. (2019) that it can replaced by an assumption of "tails decays at least as fast as the tails an exponential random variable". This is achieved by the combined effect of our assumptions (iii) and (iv). For Lin.L.3-4, Lin.L.5, is trivial as in our case Sab =
the statement of [Smucler al. (2019) that assumptions (iii) and (v) are sufficient.
Lin.E For Ein.E.I we use the statement by Smucler al. (019)) that assumption (vi) is sufficient.
For Lin.E.2, a) results from the (strictly positíve) lower bound on the variances in our assumption (vi), combined with (iii-v) for the requested upper bounds.
Lin.V For Lin. V./ is directly stated in assumption (iii).
For Lin. V.2 is trivial in our case Sab 1.
For Lin. V3 is directly imposed by our assumption (iii,, as stated by Smucler et al. (2019) for their Example 9,
B
Examples
The result discussed in Proposition is not directly intuitive. In simple words, there are two takeaways from Propo- (i) the orthogonality condition remains invariant irrespective of the causal direction between Xg and Z, and (ii) the second term in Proposition 6 suggests use a linear estimator for modeling all the relations, given that the relation between Z and Y is linear.
sition 2
To generate more intuition, we provide few examples. Let us go back again to the three variable interaction assuming the following structural equation model:
Y = 0,X; 02X2 X2 f(X,) +
(23)
X1 :=
where is nonlinear function and c1, 82 and E3 mean Gaussian noises.
Consider the case when f(z) = 2 The goal is to estimate the parameter 0₁ which we call θ, We follow the standard double ML procedure assuming policy variable X, and control X2, although the ground truth causal dependency X, - X₂ in contradiction with such setting (see Equation (4)) The estimate of 0g following the double ML procedure, which we call O2
E[X,Y] E[X]] 0₂ + 0 X Similarly, we want to estimate
E[X
X, = xx2 + y) from which we get, a X,] Q and easy to see 0 = 0,
easy to see that [[ = E 0. Hence,
Consider now the more general case where is any nonlinear function. As in the previously discussed example, the goal is to estimate fj. We have 0₂
3[X,Y|
02 00 E[X]]
X[X Similarly, = 8x We
E[X]]
22

---

Published in Transactions on Machine Learning Research (08/2022)
substitute these estimates into the orthogonality condition in Proposition
6
E
(Y - X101 X262) -ax²) 0.
[[
E
Y 11 -xabs) [x1
X,0=
X₂ 0.
E[X3]
E
(xi(0) - -1) + (a2
02)X22 + E3
)
E| X2
0 = 0,.
From the above two examples, it is clear that even though the internal relations between the variables are nonlinear, all we need is an unbiased linear estimate to estimate the causal parameter.
C Data Generation and Evaluation Metric
C.1 Data Generation
C.1.1 Causal Structure Learning Data
For every combination of number of nodes (#nodes), connectivity (p.). noise level (o²), number of observation (n), and non-linear probability (pn) (look Table C.I). 100 examples (DAGss are generated and stored as files (altogether 72.000 DAGs are simulated, comprising dataset of overall > 10GB). For each DAG, number of samples are generated by sampling noise (e in Equation 25 with variance o² starting from root of the DAG, For future benchmarking, the generated files will be made available with the code later on.
We generate DAGs (Direct Acyclic Graphs) in multiple steps: i) a random permutation of nodes chosen as a topo- logical order of DAG. ii) Based on this order, directed edges are added this DAG from each node to its followers with certain probability Ps (connectivity). iii) For each observation, values are assigned to nodes according to the topological order of the DAG in such way that each node's 's value is determined by summing over transformations (linear nonlinear with certain nonlinear probability p..) of values of its direct causes with the addition of Gaussian distributed noise. The non- linear transformation used is tanh(bæ)) 10 with 0.5 and b = 1.5. If the set of parents for the node X' is denoted as PAx as before then value assignment for node X' is as follow:
X'=++
Σ 1e(p))0X + (1 (e-(pp)))) tanh(bXX, XEPAXI
(24)
where ε ~ N(O, o2) in which 2² represents noise level. 14 (X) an indicator functions which decides between linear or non-linear contribution of X in X'. We decide the value of (p())) by generating a binary randon number which is
The resulting values in the experiments are not concentrated around zero, and they are even up to 10ks for large graphs (~~ 50 nodes). With the nonlinearity feature of thur( for relatively large values taken into account, this good representer of nonlinear relationships.
nonlincar probability (p.)
noise level (o²)
connectivity (p.)
# nodes
# observ.
0.1
100
0.01
0.3
10
0.3
500
0.1
0.5
20
0.5
1.000
0.3
50
0.5
Table C.1: Experimental Setup: In the experiments we vary the connectivity parameter, the number of nodes in the graph, the non-linear probability, the number of observations and the noise level and generate 100 graphs for each setting.
23

---

Published in Transactions on Machine Learning Research (08/2022)
Density Plot for Beta Distribution with Different Parameters
Parameters
4
d 5,8
d 1, 3
2,8=2
PDF
2
0.0
0.2
0.4
0.8
1.0
Figure C.1: Beta distribution with different parameters.
with probablity Ps. and with probability 1 1 P-Pa- The value set to 2 for the small DAGs (number of nodes equal to 10) and 0.5 for large DAGs (number of nodes equal 20 or 50) due to the value exploitation that might happen in large graphs.
Pa.
We vary and investigate the effect of non-linear relationships, the number of nodes, number of observations, effect of connectivity and noise level while simulating the data. We summarize the factors in the data generation in Table C.I
C.1.2 Inference by Regression Data
Similar to Appendix C.I.I data generation follows the random topological order but with a slight difference, i.e., with the following value assignment,
X' = Σ u(px)0X + (1 te(pn)) tanh(bX), XEPAX
(25)
where E ~ B(a, B) in which a and are parameters of beta distribution. The reason for this is that most of the inference methods exploit normality tests and this way is possible challenge them. diverse of parameters are chosen for the beta distribution signifies this point (see Figure C.I). The set of parameters of DAGs varies according to Table C.2 For each setting, 50 examples are generated and stored (15000 DAGs overall) which will be available for future studies.
nonlinear probability (p.)
beta distribution parameter
connectivity (p.)
# nodes
# observ.
(a, 3)
0.1
100
0
10
(0.5, 0.5)
0.3
0.3
20
(5,
0.5
0.5
50
(1,
1
100
(2,
200
(2, 5)
Table C.2: Experimental Setup: In the experiments we vary the connectivity parameter, the non-linear probability, the number of observations and beta distribution parameters. We generate graphs for each setting.
Code: The code for the method and data generation is available in this GitHub repository
24

---

Published in Transactions on Machine Learning Research (08/2022)
C.2 Evaluation Metric
Correctly and incorrectly inferred direct causes are considered true and false. Let the total number of true positives, false positives, true negatives ,add false negatives denoted by TP, FP, TN, and FN, we evaluate our method using following metrics:
Recall (true positive rate):
TP TPR TP+ FN
Fall-out (false positive rate):
FP FPR
FP+
Critical Success Index (CST): also known Threat Score.
TP CSI TP + FN + FP
Accuracy:
ACC TP+TN P+N
F1 Score: harmonic mean of precision and sensitivity.
2TP 2T + FP + FN
Matthews correlation coefficient (MCC): metric for evaluating quality of binary classification introduced in (Matthews 1975).
TP TN FP FN MCC VITP+ (N)(T + FP)
In some rare cases, encountered rero-divided-by-zer and divided-by-ze be-zero cases for some of these metrics. In these situations, scores are reported and respectively while Fall-out is reported and 1.
C.3
Supplementary Tables for Performance in Inferring Direct Causes
Here, additional tables for the result of the experiments provided.
C.3.1 Compared to Causal Structure Learning Methods
In this section, supplementary tables supporting superior performance of CORTH Features compared to well- established Causal and Markov Blanket discovery methods are provided (See Tables C.3 is consistent w.r.t. connectivity (Table C.5). number of nodes (Table C.3) number of observations (Table C.7 nonlin-
C.7
This superiority
earity (Table C.4). and noise (Table C.6) using different evaluation metries.
C.3.2 Compared to Inference for Regression Methods
In this part, the superiority of our method in comparison to decent Inference for Regression methods, different settings of connectivity (Table beta distribution parameters (Table C.II) number of observations (Table C.9 and nonlinearity C..ble is provided.
25

---

Published in Transactions on Machine Learning Research (08/2022)
Table C.3: Performance across all the settings for different number of nodes. Each single entry in the table is averaged 18000 simulations. Our method is almost state of the art in every case.
Number of Nodes
20
50
Method
ACC
CSI
ACC
CSL
F1
ACC
CSI
F1
ACC
CSI
F1
GES
0,935
0.890
0.911
0.854
0.730
0.779
0.743
0.442
0.526
0.698
0.245
0.323
rankGES
0,923
0.857
0.883
0.846
0.700
0.753
0.740
0.428
0.514
0.697
0.237
0.316
ARGES
0,922
0.864
0.885
0.797
0.551
0.584
0.752
0.447
0.524
0.705
0.186
0.221
rankARGES
0,914
0.838
0.861
0.793
0.537
0.572
0.750
0.435
0.514
0.705
0.181
0.216
FCI+
0,963
0.918
0,932
0.873
0.744
0.808
0.830
0.602
0.703
0.766
0.368
0,486
LINGAM
0.991
0.978
0.982
0.953
0.865
0.889
0.891
0.712
0.778
0.750
0,318
0.385
PC
0,957
0.913
0,929
0.864
0.723
0.786
0.823
0.569
0.664
0.763
0,348
0.457
rankPC
0,946
0.891
0.854
0.701
0.768
0.813
0.541
0.638
0.754
0.324
0.431
MMPC
0.868
0,586
0.823
0.494
0.535
0.790
0.489
0,749
0.260
0,350
MMII
0,929
0.878
0.905
0.841
0.675
0.739
0.767
0,432
0.507
0.725
0.218
0,281
0,883
0.613
0.623
0.855
0.563
0.601
0.824
0.501
0,580
0,759
0.310
0,388
TAMB
0.850
0.571
0.585
0.791
0.508
0.561
0.806
0.500
0.587
0.768
0.337
0.424
FasIAMB
0.883
0.614
0.624
0.858
0.571
0.611
0.828
0.511
0.593
0.770
0.326
0.409
JAMB-FDR
0.869
0.584
0.593
0.831
0.494
0.526
0.825
0.484
0.558
0.770
0.322
0.406
PCI
0.984
0.965
0.972
0.922
0.844
0.875
0.888
0.734
0,782
0.773
0.414
0.491
Lasso
0.965
0.948
0.968
0.905
0.834
0.892
0.894
0.786
0.866
0.773
0.489
0.627
CORTH Features (Curs)
0.988
0.968
0.973
0.949
0.908
0.934
0.949
0.865
0.905
0.795
0.559
0.663
Number of Nodes
Method
TPR
FPR
MCC
TPR
FPR
MCC
TPR
FPR
MCC
TPR
FPR
MCC
GES
0.934
0.056
0.891
0.790
0.090
0.711
0.502
0.088
0.436
0.304
0.083
0.221
TanKGES
0.924
0.068
0.877
0.780
0.098
0.695
0.493
0.089
0.425
0.297
0.083
0.215
ARGES
0.903
0.046
0.906
0.590
0.041
0.841
0.500
0.073
0.557
0.220
0.020
0.794
rankARGES
0.897
0.054
0.896
0.584
0.044
0.832
0.495
0.075
0,549
0.216
0.020
0.789
FCI+
0,969
0.029
0,948
0.797
0.054
0.759
0.642
0.042
0.645
0.389
0.030
0.454
LINGAM
0,991
0.007
0.988
0.886
0.008
0.934
0.770
0.055
0.759
0.391
0.072
0.471
PC
0,950
0.024
0,941
0.759
0.041
0.759
0.600
0.032
0.650
0.363
0.021
0.468
rankPC
0,944
0.039
0,925
0.750
0.053
0.734
0,580
0.034
0.629
0.341
0.024
0.427
MMPC
0.588
0.006
0.965
0.498
0.011
0.852
0.417
0.006
0.684
0.261
0.003
0.528
MMHC
0.895
0.011
0.903
0.691
0.015
0.724
0.444
0.009
0.523
0.219
0.005
0.330
0.615
0.002
0.973
0.566
0.001
0.895
0.506
0.001
0.739
0.311
0.000
0.688
IAMB
0.573
0.003
0.960
0.511
0.002
0.848
0.505
0.001
0.711
0.338
0.001
0.660
FasilAMB
0.616
0.003
0.972
0.575
0.002
0.888
0,518
0.002
0.734
0.327
0.001
0.677
JAMB-FDR
0.585
0.001
0.975
0.494
0.001
0.909
0.485
0.001
0.766
0.322
0.001
0.661
PCI
0,992
0.017
0,981
0.875
0.028
0.890
0.754
0.016
0.839
0.430
0.030
0.638
Lasso
0,999
0.074
0.949
0.944
0.119
0.817
0.954
0.147
0.794
0.681
0.148
0.488
CORTH Features (Curs)
0.999
0.016
0.986
0.952
0.044
0.906
0.884
0.011
0.894
0.609
0.101
0.567
26

---

Published in Transactions on Machine Learning Research (08/2022)
Table C.4: Performance across all the settings for different number of nonlinear probabilities. Each single entry in the table is averaged over 18000 simulations. Our method almost state of the art in every case.
Nonlinear Probability
Method
ACC
CSI
F1
ACC
CSI
F1
ACC
CSI
F1
ACC
CSI
GES
0,803
0.583
0.646
0.806
0.566
0.622
0,811
0.577
0.632
0.810
0.581
0.641
rankGES
0.796
0.559
0.625
0.801
0.546
0.605
0,805
0.556
0.613
0.805
0.561
0.623
ARGES
0.781
0.476
0.515
0.786
0.486
0.525
0.792
0.506
0.546
0.818
0.581
0.628
rankARGES
0.778
0.461
0,503
0,782
0.474
0.515
0.788
0.490
0.531
0.814
0.564
0.615
FCI+
0,827
0.599
0.674
0,860
0.663
0.745
0,872
0.685
0.764
0.873
0.685
0.746
LINGAM
0.907
0.738
0.778
0,886
0.689
0.725
0,880
0.684
0.724
0.911
0.762
0.808
PC
0.818
0.574
0.641
0.854
0.641
0.720
0.864
0.665
0.7430
0.869
0.672
0.731
rankPC
0.813
0.560
0.630
0.841
0.614
0.694
0.848
0.627
0,704
0.864
0.656
0.720
MMPC
0.775
0.372
0.416
0.809
0.439
0.503
0.818
0.462
0.528
0.828
0.479
0.523
MMHC
0.797
0.516
0.578
0.815
0.549
0.610
0.823
0,566
0.625
0.826
0.571
0.620
GS
0,806
0.450
0.491
0.828
0.494
0.554
0.835
0,510
0,571
0.851
0.534
0.576
JAMB
0.762
0,389
0,440
0.799
0.463
0.538
0.809
0.488
0.565
0.830
0.520
0,576
FastAAM
0,807
0.454
0.497
0.835
0.503
0.566
0.842
0.522
0,587
0.855
0.543
0.587
JAMB-FDR
0.796
0.418
0.457
0.818
0.456
0.511
0.827
0.481
0.538
0.853
0.529
0.578
PCT
0.853
0.674
0.720
0.897
0.746
0.789
0.905
0.763
0.806
0.911
0.774
0.805
Lasso
0.847
0.694
0.773
0.891
0.776
0.853
0.902
0.797
0.869
0.896
0.790
0.857
CORTH Features (Ours)
0.871
0.768
0.824
0.934
0.830
0.873
0.943
0.851
0.891
0.933
0.852
0.887
Nonlincar Probability
Method
TPR
FPR
MCC
TPR
PPR
MCC
TPR
FPR
MCC
TPR
FPR
GFS
0.643
0.093
0.564
0.620
0.074
0.557
0.629
0.071
0.568
0.637
0.079
0.570
rankGES
0.633
0.100
0.550
0.612
0.080
0.546
0.620
0.076
0.557
0.628
0.083
0.559
ARGES
0,514
0.041
0.789
0.526
0.041
0.793
0.547
0.043
0.791
0.626
0.055
0.725
rankARGES
0.509
0.044
0,780
0.522
0.044
0.788
0.540
0.046
0.783
0.620
0.059
0.715
FCI+
0.638
0.045
0.637
0.704
0.037
0.708
0.728
0.035
0.731
0.728
0.037
0.730
LINGAM
0.775
0.025
0.832
0.723
0.028
0.759
0.722
0.034
0,741
0.819
0.053
0.822
PC
0.605
0.037
0.649
0.672
0.027
0.707
0.695
0.025
0.728
0.702
0.029
0.734
rankPC
0.597
0.043
0,626
0.656
0.040
0.680
0.668
0.036
0.695
0.692
0.031
0.714
MMPC
0.376
0.009
0.754
0.442
0.006
0.738
0.465
0.005
0.749
0.482
0.005
0.787
МMHC
0.528
0.017
0.581
0.561
0.008
0.623
0.578
0.007
0.636
0.582
0.008
0.639
GS
0.452
0.001
0.850
0.496
0.001
0.797
0.513
0.001
0.799
0.538
0.001
0.849
JAMB
0.847
0.003
0.773
0,891
0.002
0.853
0,902
0.001
0.869
0.896
0.001
0.857
FastLAMB
0.457
0.001
0.848
0.506
0.002
0.784
0.526
0.002
0.791
0.548
0,002
0.848
JAMB-FDR
0.418
0.001
0.837
0.457
0.819
0.481
0.001
0.822
0.530
0.001
0.833
0.712
0.057
0.792
0.765
0.848
0.781
0.010
0.845
0.794
0.013
0.864
Lasso
0.823
0.130
0.684
0.907
0.120
0.778
0.926
0.116
0.800
0.921
0.122
0.787
CORTH Features (Ours)
0.840
0.119
0.730
0.849
0.007
0.872
0.870
0.008
0.888
0.885
0.038
0.863
27

---

Published in Transactions on Machine Learning Research (08/2022)
Table C.5: Performance across all the settings for different connectivities. Each single entry in the table is averaged over 24000 simulations. Our method is almost state of the art in every case.
Connectivity
0.1
0.3
0.5
Method
ACC
CST
F1
MCC
ACC
CST
F1
MCC
ACC
CST
F1
MCC
GES
0.961
0.786
0.825
0.857
0.815
0.539
0.598
0.522
0.646
0.405
0,482
0.315
rankGES
0.954
0.746
0.790
0.840
0.809
0.522
0.584
0.511
0.642
0.398
0.475
0,308
ARGES
0.965
0.794
0.828
0,876
0.805
0.456
0.501
0.726
0.612
0.286
0.330
0.720
TankARGES
0.959
0.763
0.801
0,863
0.802
0.447
0.494
0.721
0.611
0.282
0.328
0.716
FCI+
0.974
0.819
0.853
0,910
0.866
0.631
0.714
0.674
0.734
0.524
0.629
0.521
LINGAM
0,966
0.763
0.796
0,889
0.896
0.710
0.753
0.761
0,827
0.682
0.727
0.715
PC
0,975
0.819
0.849
0.921
0.861
0.609
0.689
0.676
0.718
0.486
0.588
0.516
TankPC
0,971
0.797
0.831
0,912
0.852
0.587
0.670
0.653
0.701
0.458
0.560
0,470
MMPC
0,949
0,606
0.637
0.901
0.815
0.390
0,451
0.722
0.658
0.318
0.389
0.648
МMHC
0,978
0.834
0.867
0.901
0.830
0.497
0.561
0.574
0.639
0.321
0.397
0.385
GS
0.954
0.644
0.669
0.935
0.843
0.467
0.524
0.815
0.693
0,380
0.451
0.722
TAMR
0.969
0.692
0.745
0.864
0.841
0.463
0.522
0.807
0.692
0.377
0.452
0.709
FastlAMB
0.955
0.650
0.676
0.931
0.845
0.474
0.535
0.804
0.705
0,392
0.467
0.718
JAMB-FDR
0.950
0.608
0.626
0.961
0.832
0.436
0,492
0.816
0.689
0.369
0.446
0.707
PCI
0.986
0,902
0.920
0.954
0,906
0.716
0.759
0.838
0.783
0.600
0.661
0.720
Lasso
0.976
0.886
0.925
0.926
0.876
0.725
0.811
0.737
0.800
0.682
0.778
0.622
CCRTH Features (Ours)
0.988
0.915
0.934
0,959
0.926
0.813
0.858
0.833
0.847
0.747
0.814
0.724
Table C.6: Performance across all the settings for different noise levels. Each single entry in the table is averaged over 1 14400 simulations. Our method is almost state of the art in every case.
Noise Level
0.01
0.5
Method
CSI
MCC
ACC
CSI
F1
MCC
ACC
CSI
Fl
MCC
GES
0.804
0.579
0.639
0.559
0.808
0.571
0.629
0.562
0.818
0.586
0.644
0.589
rankGES
0.797
0.557
0.619
0.548
0.802
0.552
0.613
0.551
0.812
0.565
0.625
0.577
ARGES
0.810
0.572
0.625
0.653
0.789
0.496
0.534
0.814
0.774
0.434
0.460
0.897
rankARGES
0.804
0.549
0.605
0.643
0.786
0.483
0.523
0.806
0.774
0.433
0.459
0.895
FCI+
0.843
0.617
0.691
0.674
0.865
0.678
0.753
0.717
0.874
0.766
0.740
LINGAM
0.888
0.703
0.744
0.763
0.899
0.723
0.763
0.797
0.903
0.732
0.773
0.803
PC
0.837
0.595
0.664
0.683
0.859
0.659
0.731
0.716
0.870
0.686
0.752
0.745
rankPC
0.831
0.584
0.657
0.653
0.845
0.626
0.699
0.688
0.856
0.655
0.724
0.714
MMPC
0.796
0.405
0.456
0.762
0.812
0.453
0.510
0.756
0.825
0,480
0.533
0,780
MMHC
0.806
0.526
0.585
0.605
0.818
0.557
0.615
0.626
0.829
0.586
0.639
0.652
GS
0.820
0.468
0,518
0.824
0.836
0.513
0.566
0.819
0.846
0.538
0.586
0.833
JAMB
0.784
0.421
0.483
0.779
0.807
0.485
0.552
0.769
0.823
0.523
0.586
0.790
FasIAMB
0.821
0.469
0.520
0.819
0.842
0.526
0.582
0.814
0.852
0.548
0.600
0.828
JAMB-FDR
0.810
0.432
0.478
0.834
0.831
0.492
0.545
0.825
0.841
0.514
0.563
0.841
PCI
0.873
0.690
0.730
0.819
0.901
0.760
0.801
0.846
0.906
0.777
0,815
0.854
Lasso
0.868
0.728
0.807
0.725
0.891
0.780
0.852
0.779
0.898
0.794
0,861
0.793
CORTH Features (Ours)
0.899
0.789
0.839
0.795
0.929
0.842
0,883
0.858
0.934
0.854
0.891
0.866
28

---

Published in Transactions on Machine Learning Research (08/2022)
Table C.7: Performance across all the settings for different number of observations. Each single entry in the table is averaged over 24000 simulations. Our method is almost state of the art in every case.
Number Obscrvations
100
500
1000
Method
ACC
CSI
FI
MCC
ACC
CSI
F1
MCC
ACC
CSI
FI
MCC
GES
0.797
0.524
0.588
0.539
0.811
0.593
0.650
0.572
0.815
0.612
0.666
0.583
rankGES
0.788
0.495
0.561
0.522
0.806
0.576
0.636
0.564
0.810
0.595
0.652
0,573
ARGES
0.780
0.446
0.489
0.786
0.799
0.535
0.576
0.773
0.803
0.555
0.595
0.764
TankARGES
0.776
0.428
0.473
0.778
0.795
0.523
0.566
0,766
0.800
0.542
0.584
0.757
FCI+
0.837
0.589
0.671
0.652
0.865
0.684
0.755
0.720
0.871
0.702
0.771
0.732
LINGAM
0.840
0.578
0.650
0.678
0.908
0.719
0.743
0.825
0.941
0.858
0.883
0.862
PC
0.830
0.568
0.642
0.661
0.858
0.662
0.732
0.719
0.866
0.684
0.752
0.733
TankPC
0.821
0.544
0.617
0.632
0.849
0.639
0.711
0.696
0.855
0.660
0.733
0.707
MMPC
0.771
0.323
0.368
0.787
0.819
0.476
0.534
0.739
0.832
0,515
0.575
0.745
MMIC
0.800
0,495
0.557
0.579
0.820
0.570
0.625
0.633
0.826
0.587
0,642
0,647
GS
0.793
0.375
0.427
0.785
0.842
0.540
0.592
0.834
0.856
0.577
0.625
0.852
TAMB
0.745
0.316
0.390
0.705
0,815
0.510
0.574
0.794
0.835
0.556
0.614
0.821
FasiTAMB
0.803
0.401
0.461
0.770
0.844
0.541
0.593
0.833
0.857
0.574
0.623
0.850
JAMB-FDR
0.783
0.325
0.372
0.825
0.837
0.523
0.578
0.815
0.850
0.564
0.613
0.843
PCI
0.829
0.551
0.594
0.804
0.914
0.812
0.853
0.842
0.931
0.855
0.893
0.866
Lasso
0.870
0.729
0.812
0.732
0.889
0.778
0.848
0.773
0.893
0.786
0.854
0.780
CORTH Features (Ours)
0.883
0.710
0.780
0.754
0.935
0.874
0.906
0.874
0.942
0.891
0.920
0.887
Table C.8: Performance across all the settings for different nonlinear probabilities. Each single entry in the table is averaged over 3750 simulations.
Noulincar Probability
0
0.3
0.5
Method
TPR
CSI
F1
TPR
CSI
F1
TPR
CSI
F1
TPR
CSI
F1
Standard Regression
0.149
0.103
0.139
0.166
0.112
0.141
0.175
0.108
0.136
1.000
0.801
0.801
Lasso
0.237
0.116
0.176
0.285
0.126
0.202
0,360
0.165
0.263
1.000
0.046
0.046
Debiased Lasso
0.238
0.117
0.178
0.267
0.112
0.178
0,300
0.124
0.202
1.000
0.050
0.050
Forward Stepwise Reg_active
0.174
0.112
0.162
0.157
0.110
0.167
0.194
0.129
0.190
1.000
0.329
0.329
Forward Stepwise Reg_all
0.04
0.039
0.060
0.062
0.059
0.085
0.089
0.086
0.116
1.000
0.861
0.861
SARS_active
0.073
0.054
0.094
0.104
0.078
0.131
0.118
0.081
0.134
1.000
0.382
0.382
LARS_al all
0.017
0.016
0.028
0.030
0.029
0.048
0.039
0.037
0.057
1.000
0.866
0.866
CORTH Features (Ours)
0.481
0.287
0.407
0.436
0.258
0.366
0.364
0.220
0.313
1.000
0.610
0.610
Table C.9; Performance across all the settings for different number of observation. Each single entry in the table is averaged over 3000 simulations.
Number of Observations
10
20
50
100
200
Method
CSI
F1
CSI
F1
CSI
F1
CSI
F1
CSI
F1
Standard Regression
0.250
0.250
0.250
0.250
0.250
0.250
0.263
0.272
0.392
0.499
Lasso
0.075
0.117
0.100
0.155
0.127
0.192
0.131
0.196
0.132
0.198
Debiased Lasso
0.066
0.102
0.086
0.134
0.115
0.173
0.122
0.179
0.114
0.171
Forward Stepwise Reg active
0.193
0.199
0.161
0.177
0.103
0.134
0.071
0.111
0.322
0.439
Forward Stepwise Reg_al all
0.222
0.224
0.222
0.226
0.229
0.236
0.244
0.257
0.389
0.458
LARS active
0.193
0.200
0.171
0.191
0.143
0.177
0.090
0.128
0.160
0.242
LARS all
0.218
0.222
0.217
0.221
0.226
0.231
0.230
0.235
0.293
0.342
CORTH Features (Oura)
0.125
0.173
0.250
0.314
0.353
0.445
0.443
0.548
0.550
0.640
29

---

Published in Transactions on Machine Learning Research (08/2022)
Table C.10: Performance across all the settings for different connectivities. Each single entry the table averaged over 5000 simulations.
Connectivity
Method
TPR
F1
TPR
CSI
TPR
Standard Regression
0.395
0.267
0.290
0.372
0.292
0.313
0.350
0.284
0.310
Lasso
0.789
0.211
0.314
0.353
0.094
0.150
0.269
0.035
0.052
Debiased Lasso
0.787
0.211
0.314
0.296
0.054
0.087
0.270
0.037
0.055
Forward Stepwise Reg active
0.437
0.182
0.231
0.352
0.156
0.198
0.355
0.171
0.208
Forward Stepwise Reg_all
0.356
0.315
0.343
0.235
0.252
0.263
0.234
0.245
LARS active
0,360
0.149
0.192
0.312
0.145
0,183
0.340
0.169
0.196
LARS
0,289
0.246
0.264
0.267
0,233
0,246
0.259
0,231
0.239
CORTH Features (Ours)
0.494
0.367
0.417
0.540
0.274
0.355
0.677
0.390
0.499
Table C.11: Performance across all the settings for different parameters for the Beta distribution. Each single entry in the table averaged over 3000 simulations.
Beta Distribution Parameters (a, ()
(0.5, 0.5)
(1.3)
(2, 2)
(2, 5)
(5.
Method
CSI
CSI
F1
CSI
F1
CSI
F1
Standard Regression
0.282
0.305
0.280
0.304
0,280
0.303
0,280
0.303
0.283
0.306
Lasso
0.109
0.168
0.106
0.164
0.116
0.174
0.105
0.163
0.129
0.189
Debiased Lasso
0.096
0.148
0.092
0.143
0.152
0.093
0.144
0.173
Forward Stepwise Reg_active
0.169
0.211
0.168
0.211
0.214
0.215
Forward Stepwise Reg_all
0,261
0,279
0.257
0.275
0.265
0.284
0.278
0,284
LARS active
0.161
0.198
0.150
0.184
0.146
0.182
0.155
0.191
0.157
0.193
LARS_ai
0.235
0,248
0.236
0.249
0.241
0.254
0.238
0.251
0,234
0.247
CORTH Features (Gurs)
0.336
0.417
0.330
0.411
0.354
0.434
0.336
0.417
0.363
0.441
D Real-World Data Experiment-Covid19
D.1
Preprocessing
The preprocessing stage for this dataset is the same as (Schwab al., 2020) except that, for each target variable upsampling is used resolve data imbalance.
D.2
Results
The results obtained by leveraging CORTH Features is suprisingly consistent with (Schwab et al. 2020) which demon- strates the ability of this method in feature selection. The selected features are indicated in Tables D. to D.4
30

---

Published in Transactions on Machine Learning Research (08/2022)
Table D.1: Ranks of the features based on the times being predicted as direct causes of SARS-Cov-2 exam result out of 1000 different runs of our propsal approach. Not mentiond features were not predicted even once, note that preprocessed dataset has 331 fea- tures.
Table D.2: Ranks of the features based on the times being predicted as direct causes of Patient addmited to regular ward out of 1000 different runs of our propsal approach. Not mentiond features were not predicted even once, note that preprocessed dataset has 331 features.
Isneet aps quantie
Pinuer ago quaette
RC03 gas andyy's
bload realyii
Centa plasnyttonnoon
zi saazes blod seabyie
6.99
MISSING Lusis Debyingeras
0287
MISSING Sees pladee vileas:
6.992
Alasse
03543
MINAINN Lari tilgingra
2300
kigaIN Trimation somratted Far TNE
0391
Aerrieedd
695
Senm Claone
Mothrugees
6004
ji noaues Hol ga ssalyss
3001
6.50
Base exoos veasus Hood gn aealyss
0.341
Ladesgta
0.794
MISSING Ancead
0.33
CO2 snmaal Hoed zo unlviia
6499
Liiss Dasety
0334
Mayresier
0234
MISSING
0T13
6092
plonptokinne
5,00
BOU antrral babodd aaayyn
0046
0188
6000
07B
MERCSTN Magezdre
KINOOR Frniii
36
HD anunal Noed ga prolyde
0m1
35
Lice
0784
MISSING Adiced Fol
6X04
21
Rogkrann Seaggial VeaggialVVaaaaaaan Veasdeacas)
0002
6008
MISSINI: Magrecum
MISSING Feman
Repirrocy Syarydal fexcist
6014
MINSINU Alhamn
2010
Coueler phw.phokiaau.crc CPC
6010
500
Apoiiin
Nesteeptil
Infinrroa poutine
3007
28
Red Neod sell deeteetion with RDW
Consasicas HKUI Jeecial
0004
29
Unine Pretán eas
MISAINN Seekerz
0.022
Intteroa rape tre poatce
2a22 antral blod pas andyaa
Hb etorainn warse aloo yon Beatys
01822
Iifloros deterrd
0001
Iaflenza driecied
datsted
MISSING: Serom Ghooe
6001
Aup atuld antiane
Litee Sevrri sun
Table D.3: Ranks of the features based on the times being predicted direct causes of Patient addmited to semi-intensive unit out of 1000 different runs of our propsal approach. Not mentiond features were not predicted even once, note that preprocessed datasct has 331 features.
Table D.4: Ranks of the features based the times being predicted as direct causes of Patient addmited to intensive care unit out of 1000 different runs of our propsal approach. Not mentiond features were not predicted even once, note that preprocessed dataset has 331 features.
Bank Persue
Voae of heieg Psodicted is Daveet Uspas
loiy Poolctal as Dreei Chasc
Psncet quarilc
Patcet age quanto
C'metrine
MIsSiNG Meus plaidet volame
MISSING Lests Deby drogenac
gworelyds
Tresl veznee Nood gur analyde
8003
Magueiiuu
Alasis
Gunnu platang luan
Alazme timartnnne
Magnectrt:
cX2 anenal lood andoysa
HCT3 weven hloel gas avdycs
0.099
MISSING Losti Delgdrgeneee
Crnatrane
Ki Palcei Narra
0.796
POl veasus ga cnalpis
XISSING Artcea Firt
0,93
gas asdyyis
Hazz racmic vennue Nest ger zralyst
6374
pDx vancec Mose ga aralpux
MISSINI Anera
0707
Inkerakrraal nomeation
MISSINU Insmationa
INK
0396
Ro Rolooo gol dovilation wilh RDW
0300
Liire Dardly
0:14
29
Creaits phoptokinas CMK
0.391
6317
Baic espos vesses bood po stabois
0332
Theepter
616
3
Line Denily
6783
21
Repiraec Syacyial Vhas leacial
6061
Famiin
0174
21
MINSING Meas plasde velemee
6486
Metguoenoonn
0132
MINNING Pentin
675
Phanptor
sanus Mand ga mnajaic
Steppr proding
6013
MISSING Megredirs
5017
HINI dstetal
8016
Frversss nutha mg.l.
0006
ripit pooncc
2
Angalilo tranasuane
0003
25
Regiranyy Spnylial Yeas deacko)
3000
Convestias HKU1 deecie
6012
24
Köbnerdii Prieri Noree
0107
Asgantic
MaSSINNG Manzcotra
5004
6008
Antria Lacte Aca
MISSING Allunn
000
Ral bloe edi dosilatoo with RDW
lite Kerne 1istu Shore
Todtsphil
isteasiza doscied
0.05
Intheros diteae
0/01
deet
Inimmtioral nomaltadd min INK
0007
Cyyits Anerio
3
Ledoostes
0JMI
Smnjo A poliis
31


Reviews:
[{"review_id": "O2nLKSBEz8", "reviewer": "Reviewer_fVm5", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strengths:\n - Overall, the idea of (1) tailoring the solution to the special causal discovery problem of discovering the parents of one variable, and (2) using double ML to construct valid hypothesis tests is nice. This is certainly more efficient than attempting to learn a DAG for all variables, and could work a wide range of cases where nuisance functions needs to be estimated with regularization.\n - The method could scale to high-dimensional cases where regularization like LASSO is necessary to estimate E[Y | X], e.g., the case considered in [Farrell 2015](https://arxiv.org/abs/1309.4686) or [Belloni et al 2014](https://academic.oup.com/restud/article-abstract/81/2/608/1523757?redirectedFrom=fulltext).\n\nWeaknesses:\n - The presentation is far more complex than it needs to be, and does not clearly highlight the advantage of the proposed method. The paper is written as though the Double ML approach targets a different parameter from what other methods would target, but this is not true. The stated assumptions of the problem allow for an immediate reduction of the problem to testing whether the regression coefficient of each X is zero in the linear regression of Y on all X's simultaneously (or, equivalently, whether the partial correlation between Y and each X is zero, after projecting out the rest of the X's). See the [Frisch-Waugh-Lovell Theorem](https://en.wikipedia.org/wiki/Frisch%E2%80%93Waugh%E2%80%93Lovell_theorem) for this equivalence. The causal structure among the X's would not matter for any estimator of E[Y | X], since X is being conditioned on, and the structural noise terms are assumed to be independent of all X's. Rather, the strength of Double ML is providing (potentially) better statistical efficiency and a central limit theorem for hypothesis testing.\n - Some conditions are likely missing. In particular, one needs to have some notion of faithfulness in the X-Y relationship to perform structure discovery, even if the relationships among the X's are irrelevant. Similarly, in the assumptions for the statistical test, equation (9) does not eliminate the case of M-bias, where E[U | X_j] = 0 marginally, but E[U | X_j, X_k] \\neq 0, because conditioning on X_k opens a backdoor path. It seems to me that the condition in equation (5) is the one that's actually desired, where E[U | Z, D] = 0.\n - The methodology could be presented more clearly. It seems that the method actually centers around a set of partial covariance tests, not on estimation of parameters in the partial linear model, yet the methodology section begins with a description of estimating treatment effects in the partial linear model. The detour into the partial linear model and treatment effects seems irrelevant.\n - Too many details of the simulations are hidden in the appendix.\n - I am not sure that the baselines being compared to are appropriate. Some straightforward baselines based on hypothesis tests from ordinary linear regression or de-biased LASSO (Javanmard and Montanari) should be added.\n - The real data example does not add any insight, given the level of detail given in the paper.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "The authors propose a straightforward technique for discovering the causal parents of a target variable, under the assumptions that the structural equation for the target variable is linear in its parents, and all potential parents of the target variable are non-descendants. The technique focuses on each potential parent in turn, and constructs a hypothesis test of whether the conditional covariance (given all other variables) between the potential parent and target is zero, using the double ML framework. The double ML framework uses a Neyman orthogonal score function and cross-fitting to de-bias an estimate of the parameter of interest when nuisance functions are estimated with regularization. The authors construct many simulated DGP's and show that this strategy does a better job of identifying causal parents than more generic causal structure learning algorithms.", "claims_and_evidence_yn": "", "claims_explanation": "", "audience_yn": "", "audience_explanation": "", "requested_changes": " - Reframe the contributions to focus on debiasing and valid hypothesis tests, rather than identification under arbitrary causal structure in the X's. As stated under Weaknesses, if I understand correctly, the assumptions in the paper imply that ordinary linear regression of Y on all X's would identify the desired causal effects (the controlled direct effect of each X_j on Y), and would thus enable hypothesis testing about whether these effects are zero using more straightforward hypothesis tests, such as those that are built into linear regression routines. These approaches put no conditions on the causal structure of X, but only require that the covariance of X be non-degenerate. It seems that some of the examples attempt to show that the DML approach targets a distinct parameter from what one might target with more naive methods, but this is not the case. For example, the score equation given at the top of page 6 is exactly equivalent to a score equation that is solved by the standard linear regression coefficient on X_k in the multiple regression of Y on all X's, by the Frisch-Waugh-Lovell Theorem. If I am mistaken here, the authors need to provide a better justification of why the double ML procedure targets a unique estimand.\n\n - (Extension of last point) The more compelling contribution of double ML approaches is that it enables valid hypothesis testing in contexts where nuisance functions need to be estimated with regularization. In linear contexts, this could correspond to cases where the X's are high-dimensional. So it might make sense to re-focus the presentation on high-dimensional contexts.\n\n - Make sure that the conditions for validity are consistent and complete. As mentioned above, a faithfulness assumption seems to be missing, and assumptions about independence of the noise term are not consistent across the paper.\n\n - Focus the presentation on the particular hypothesis testing task, rather than the partial linear model.\n\n - Provide more details of the simulation in the main text. For example, it is very difficult to interpret what the various simulation parameters are from the main text alone.\n \n - Include baselines from standard linear regression and de-biased LASSO (Javanmard and Montanari).\n\n - Either eliminate the real data example, or provide a more thorough description of the problem and associated variables in the main text.", "broader_impact_concerns": "The real data example concerns COVID diagnosis, which could be a sensitive topic. The paper caveats the inconclusiveness of the analysis relatively well, but if the authors expand on this example, they should remain careful about the strength of claims made here.", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "YeRhyYrGmt", "reviewer": "Reviewer_9rrS", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strength: \n+ In general, the paper is well written. The idea and the method are clearly introduced. For example, I really like the example in Sec. 2.3, which makes the point much easier to be understood.\n+ The rigorous analysis and theoretical extension of the double ML are appreciated.\n+ The statistical test is appreciated. This is indeed necessary for causal discovery or causal feature selection, which is of lack in current causal methods with ML techniques.\n\nWeakness:\n\nThe theoretical analysis and the derivation in Sec. 2.3 and Sec. 2.4 requires more effort to be self-contained and for the clarification and the transparency. The proof and the theoretical results seem right to me, but due to the mentioned weakness, I cannot always follow the proofs and the derivations. It is likely that I may miss or misunderstand some parts, but please clarify the corresponding points if possible to fit a large range of the (potential) readers of TMLR and optimize their benefits of reading the paper. In the following “Requested changes”, I show the detailed comments. ", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "The authors aim at the causal feature selection by leveraging the connection with double ML. \n\n1. Assumptions: the model generally assumes that the variable of interest has no children and that the direct causal effects on the variable of interest are linear while the other relationships can be nonlinear and/or contain cycles. \n2. Under the assumptions, the authors extend the conclusion in double ML so that it can be applied in more situations that are necessary for the causal feature selection, as in Prop. 2. And then based on the extension, proposed an efficient algorithm as in Alg. 1. \n3. Furthermore, a statistical test is provided for determining the direct causal effect, as in Prop. 3.\n4. Based on the synthetic data generated according to Eqn (18), the proposed method has state-of-art results. The method is also applied to the COVID-19 dataset.", "claims_and_evidence_yn": "", "claims_explanation": "", "audience_yn": "", "audience_explanation": "", "requested_changes": "**Sec. 2.1** \n\nIt seems that the authors assume that the audience has certain knowledge about Neyman's orthogonality condition, i.e., Eqn (3) and (4). However, based on the paper, I don’t know and would like to ask for clarification on several points:\n\n1. how to have the estimator such that $\\hat{\\eta}$ satisfies Eqn (3).\n    * based on the algorithm and the proof for Example 1, it doesn’t seem to require an extra procedure, but only fitting nonlinear regressions.\n2. what is the expression of the estimator $\\theta$ and how to get it from Eqn. (4).\n    * based on the algorithm and the proof for example 1, it seems to have an expression. So I would like to know what is it and how to get it, especially for the audience without sufficient knowledge about Neyman's orthogonality condition and the related works.\n\n**Sec. 2.2**\n\nIt would be much better to see the robustness of the proposed algorithm, especially when the assumptions don’t hold. It can be shown via experiments or simply by providing a discussion about them. For example, \n\n1. What if in Z of Eqn. (10), there is a covariate that is not a confounder or a chain. It can be the independent noise or only a parent of Y.  \n2. What if Y is the cause of some variables.\n3. What if there are unknown confounders. \n\nCan the method still work and what would be the consequence?\n\n**Proof of Prop. 2**\n\n- To be self-contained, it would be much more convenient to include Eqn 2.7 and 2.8 in (Chernozhukov et al., 2018a).\n- I am not sure how to get the Neyman orthogonal score? It would be better to show the procedures of how to get $\\mu$ and $J$. Due to a lack of derivation steps, I didn’t check the correctness of the derivation. Since the derivation is basically plugging in terms, it doesn’t seem to be wrong. But I highly suggest making each step clear for the derivation.\n- It would be better to make it clear what is the unbiased estimator expression used in the algorithm. And how to get the expression (maybe the orthogonality condition implies so, but please make it clear).\n- The current proof of Prop. 2 is based on linear functions with $\\theta$, $\\beta$ and $\\gamma$. Note that the algorithm and the claims of the paper are based on nonlinear m and g, while the proof of Prop. 2 is based on linear m and g.\n    - I am fine with the focus of work is to prove the linear case, but I would suggest making the claims clear, i.e., what has been proved and what hasn't. And why so. In fact, I didn’t see a reason why the nonlinear cases don't hold or are more difficult? I may miss some parts about the nonlinear case; otherwise, I would suggest including the proof of the nonlinear case to make it transparent.\n- In the proof of Prop. 2, it assumes the gaussianity of the assignments of Y. It seems that there is no discussion or introduction of this in the main paper. If so, please make the claim and the assumption clear.\n\n**Proof of Prop. 3**\n\nFor the proof of Prop. 3, I can understand and follow the proof of a), but I feel it is less straightforward to see how b) is proved and to understand why c) is proved.\n\n**Variance of Empirical Estimates $\\chi$**\n\nI would suggest including the important conclusions from the other papers to make them self-contained and easier to read. Moreover, I didn’t see the connection that why after deriving Eqn (11) and the given that Eq (12) is LASSO-type estimators, one can get the asymptotic expression.\n\n**Experiments**\n\nIt requires a justification of the data generation process, especially, Eqn. (18).\n\n-----\n**Minor comments**\n* From page 24 to page 83, they are basically the experiments plots under different configurations without discussion and analysis. It doesn’t seem to be necessary to attach all the experiments plots in this way. I would suggest summarizing them a bit and to show what are the take-away messages. If there are no more things to discuss or analyze, maybe consider keeping the typical ones and removing the others.\n* In Figure 1, there are X1 to X6 involved while the caption shows {X1, ..., X11}.\n* In Algorithm 1, Line 7, is it supposed to be $\\hat{V}_{ij}^{[k]}$?\n* In the paragraph after Example 1, it seems that $\\theta_2$ fits the double ML model because X1 is a confounder, which can be estimated directly. While the text seems to show that $\\theta_1$ is straightforward to be estimated.\n* I didn't understand the paragraph above Sec. 2.4. \n\n", "broader_impact_concerns": "Looks good. I didn't have any concern on the ethical implications of the work.", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "CpzcSp86v5", "reviewer": "Reviewer_Xh2v", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strengths:\n\nSelecting causal features is an important but challenging problem. The authors extend the method of Double Machine Learning to learn causal parents of the target variable of interest is reasonable.\n\nThe analysis and the algorithm are sound and presented in a logical way.\n \n\nWeaknesses:\n\nThe assumptions of linear direct causal effects of the response variable and no hidden confounder for the response variable seem a bit strict.\n\nLack of some discussion in related work and experiments,  e.g., Markov blanket-based methods.\n", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "In this paper, the authors consider the problem of selecting causal features from an observational dataset. In detail, the authors provide a new method to learn the causal parents of a response variable from the explanatory variables. Besides, the authors show that all parents of the target variable of interest can be identified under certain assumptions. The experimental results demonstrate the effectiveness of their approach. ", "claims_and_evidence_yn": "", "claims_explanation": "", "audience_yn": "", "audience_explanation": "", "requested_changes": "Related work: As far as I know, finding the Markov blanket of the response of interest is another line for causal feature selection. Unfortunately, I can not find any descriptions in this paper. One should discuss the connections and differences between their method and  Markov blanket-based methods. In fact, in the paper's setting (there is no hidden confounder for the response variable and all explanatory variables are not the descendants of the response variable), one can identify all parents using the existing Markov blanket-based methods. \n\nThe authors may need to compare their method with Markov blanket-based methods in the Experiments section.  \n\n﻿\"Extensions of the IV-based approach ….are the closest known result to discovering direct causal parents\" in Section 1.1: Is this sentence reasonable? Can you explain how?\n\nIn the last paragraph on Page 3: no definition for W.\n\nIn the first equation on Page 4: no definitions for i and n.\n\nIt would be helpful to give the Conclusion Section.", "broader_impact_concerns": "N/A", "additional_comments": "", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Meta-review: This paper studies feature selection for a linear regression SEM when the covariate is in high dimensions. The authors show that under the assumption that no feature in X is a child of Y, they can derive a double ML method (Chernozhukov et al. 2018) to recover the direct parents of Y in X. Computational complexity and asymptotic normality analyses are provided. Upon revision, a variety of experiments have been conducted to demonstrate the improvement of feature selection over other feature selection methods.

Reviewers all recommended "accept/leaning accept" to this submission, so an agreement is reached.

I would ask the authors to again consider the suggestions from reviewer fVm5. Please clarify your assumptions and address fVm5's comment on faithfulness, and add a discussion on comparing your approach with debiased LASSO (not just adding experimental comparisons only). This will help clarify your approach better.
Decision: Accept (poster)