ANOCE: Analysis of Causal Effects with Multiple Mediators via Constrained Structural Learning
Authors: Hengrui Cai
Rui Song
Wenbin Lu
Venue: ICLR 2021 Poster (2021)
Keywords: Causal network
Constrained optimization
COVID-19
Individual mediation effects
Structure learning

Abstract: In the era of causal revolution, identifying the causal effect of an exposure on the outcome of interest is an important problem in many areas, such as epidemics, medicine, genetics, and economics. Under a general causal graph, the exposure may have a direct effect on the outcome and also an indirect effect regulated by a set of mediators. An analysis of causal effects that interprets the causal mechanism contributed through mediators is hence challenging but on demand. To the best of our knowledge, there are no feasible algorithms that give an exact decomposition of the indirect effect on the level of individual mediators, due to common interaction among mediators in the complex graph. In this paper, we establish a new statistical framework to comprehensively characterize causal effects with multiple mediators, namely, ANalysis Of Causal Effects (ANOCE), with a newly introduced definition of the mediator effect, under the linear structure equation model. We further propose a constrained causal structure learning method by incorporating a novel identification constraint that specifies the temporal causal relationship of variables. The proposed algorithm is applied to investigate the causal effects of 2020 Hubei lockdowns on reducing the spread of the coronavirus in Chinese major cities out of Hubei. 

Published as a conference paper at ICLR 2021
ANOCE: ANALYSIS OF CAUSAL EFFECTS WITH MULTIPLE MEDIATORS VIA CONSTRAINED STRUCTURAL LEARNING
Hengrui Cai
Rui Song
North Carolina State University Raleigh, USA hca i 5@ncsu.edu
North Carolina State University
Raleigh, USA
rseng@ncsu. edu
Wenbin Lu
North Carolina State University Raleigh, USA wlu4@ncsu.ent
ABSTRACT
In the era of causal revolution, identifying the causal effect of an exposure the outcome of interest is an important problem in many areas, such as epidemics, medicine, genetics, and economics. Under a general causal graph, the exposure may have a direct effect on the outcome and also an indirect effect regulated by a set of mediators. An analysis of causal effects that interprets the causal mechanism contributed through mediators is hence challenging but on demand. To the best of our knowledge, there are no feasible algorithms that give an exact decomposition of the indirect effect on the level of individual mediators, due to common interaction among mediators in the complex graph. this paper, we establish new statistical framework to comprehensively characterize causal effeets with multiple mediators, namely, ANalysis Of Causal Effects (ANOCE), with a newly introduced definition of the mediator effect, under the linear structure equation model. We further propose a constrained causal structure learning method by incorporating a novel identification constraint that specifies the temporal causal relationship of variables. The proposed algorithm is applied to investigate the causal effects of 2020 Hubei lockdowns on reducing the spread of the coronavirus in Chinese major cities out of Hubei.
INTRODUCTION
In the era of causal revolution, identifying the causal effect of an exposure on the outcome of interest is an important problem in many areas, such as epidemics (Hernán 2004). medicine (Hernán et al.) 2000), education (Card) 1999). and economics (Panizza & Presbitero 2014). Under general causal graph, the exposure may have direct effect the outcome and also indirect effect regulated by set of mediators (or intermediate variables). For instance, during the outbreak of Coronavirus disease 2019 (COVID-19), the Chinese government has taken extreme measures to stop the virus spreading such as locking Wuhan down on Jan 23rd, 2020, followed by 12 other cities in Hubei, known as the "2020 Hubei lockdowns". This approach (viewed the exposure), directly blocked infected people leaving from Hubei; and also stimulated various quarantine measures taken by cities outside of Hubei (as the mediators), which further decreased the migration countrywide in China, and thus indirectly control the spread of COVID-19. Quantifying the causal effects of 2020 Hubei lockdowns on reducing the COVID-19 spread regulated by different cities outside Hubei is challenging but of great interest for the current COVID-19 crisis. An analysis of causal effects that interprets the causal mechanism contributed via individual mediators is thus very important.
Many recent efforts have been made on studying causal effects that are regulated by mediators. Chakrabortty et al. (2018) specified the individual mediation effect a sparse high-dimensional causal graphical model. However, the sum of marginal individual mediation effect is not equal to the effect of all mediators considered jointly (i.e. the indirect effect) due to the common interaction among mediators (VanderWeele & Vansteelandt 2014). Here, 'interaction' means that there exists at

---

Published as a conference paper at ICLR 2021
M
M₁
Y
M,
M
M,
(a) causal graph with parallel mediators.
(b) causal graph with interacted mediators.
Figure 1: Illustration of causal graphs with different types of mediators, where A is the exposure,
(M,,"., ". Mp} are mediators, and Y is the outcome of interest.
least one mediator that is regulated by other mediator(s) (see Figure [b for illustration), in contrast to the simple 'parallel' case (shown in Figure Lak. Vansteelandd & Daniel (2017) considered an exact decomposition of the indirect effect with a two-mediator setting based on the conditional densities of mediators, while there was no feasible algorithm provided to solve their proposed expressions yet. Therefore, a new framework with a computational friendly algorithm that gives an exact decomposition of the indirect effect on the level of individual mediators is desired under the complex causal network.
To estimate the underlying causal network, structure learning algorithms of the directed acyclic graph (DAG) are widely used. Popular methods such as the PC algorithm (Spirtes et al., (2000) that uses conditional independence tests to examine the existence of edges between each pair of variables, require strong assumptions and thus have no guarantee in the finite sample regime. Recently, Zheng et al. (2018) opened up another class of causal discovery methods by directly formulating pure optimization problem over real metrics with a novel characteristic of the acyclicity. Yu et al. (2019) further extended Zheng al. (2018)'s work with deep generative model, and showed better performance on the structure learning with weaker assumptions on the noise. See more follow-up works in Lachappele al. (2019) and Zhu & Chen (2019). However, the current cutting-edge methods neglect the temporal causal relationship among variables, and thus cannot appropriately represent the causal network with pre-specified exposure and outcome.
In this paper, we consider establishing new statistical framework to comprehensively characterize causal effects with multiple mediators, namely, ANalysis Of Causal Effects (ANOCE), under the linear structure equation model (LSEM). Specifically, individual mediators, the natural direct effect and the natural indirect effect for a mediator, denoted as DM and IM, respectively. Our proposed DM can be interpreted the direct effect of particular mediator on the outcome that is not regulated by other mediators, while the IM is the indirect effect of the mediator controlled by its descendant mediators. We prove that the DM is valid in the sense that it exactly decomposes the indirect effect of the exposure on the outcome, followed by an ANOCE table to explain different sources of causal effects. To bridge the cutting-egge graphical learning approaches with the temporal causal relationship of variables, we extend the variational auto-encoder (VAE) framework in Yu et al. (2019) with a novel identification constraint that specifies the topological order of the exposure and the outcome. The proposed constrained VAE algorithm is AVAE'.
propose two causal effects the level of
then used to estimate causal effects defined in our ANOCE table, named
Our contributions can be summarized in the following three aspects: 1). Conceptually, we define different sources of causal effects through mediators with a newly introduced definition of direct and indirect mediator effects, and give an exact decomposition of the indirect effect on the level of individual mediators, under the linear structure equation model.
2). Methodologically, we incorporate the background knowledge (the temporal causal relationship among variables) when using an optimization approach to the causal discovery. Such prior knowledge can be generalized for any measured variable and the possible set of their parents. Our proposed constrained structural learning can be casily extended to other score-based algorithms.
3). Practically, extensive simulations conducted to demonstrate the empirical validity of the
2

---

Published as a conference paper at ICLR 2021
proposed algorithm and its competitive performance among existing causal discovery algorithms. Our method is applied to investigate the causal effects of 2020 Hubei lockdowns on reducing the COVID-19 spread in China, by quantifying the individual effect for each city.
2
BACKGROUND AND RELATED WORK
2.1 GRAPH TERMINOLOGGY
Consider a graph g == (X, (X,,) F)) with a node set X and an edge set E. There is at most one edge between any pair of nodes. If there is edge between X, and X,, then X, and Xj are adjacent. A node X; is said to be parent of X2 if there is a directed edge from X; to Xj. Let the set of all parents of node Xj in 9 PAx, (9). A path from X, to X; in G is a sequence of distinct vertices, J C V such that ag X, and d₂ = Xj. A directed path from X; to X, is a
= { 01, path between X; and X, where all edges are directed toward X,. A directed cycle formed by the directed path from X; to X, together with the directed edge X, to X4, A directed graph that does not contain directed cycles is called a directed acyelic graph (DAG). directed graph acyelic and only if it has a topological ordering. Suppose DAG G = (X, F.) that characterizes the causal represents a random vector and
relationship among [X] = d nodes, where X = [X X2. ... , an edge X; + X, means that X, is direct cause of Xj. Let B { } Iisid,siss be d x d matrix, where bi.5 is the weight of the edge X, + X, E E, and bi.j = otherwise. Then, we say that G (X, B) is weighted DAG with the node set X and the weighted adjacency matrix B (the edge set E is nested in B).
2.2 RELATED WORK
Our work connects to the literature of the causal graphical model. Pearl al. (2009) provided comprehensive review of recent advances in the analysis of causes and counterfactuals using 'do- operator' by graphical methods. Later, Maathuis al. (2009) started use an unknown DAG without hidden variables to estimate the causal effects from the high-dimensional observational data. Nandy al.] (2017) extended the work of Maathuis et al. (2009) with the linear structure equation model. followed by the individual mediation effect defined in Chakrabortty et al. 2018). All of these models rely on the PC algorithm to search the Markov equivalence class of the partial DAG, and usually require strong assumptions due to the computational limit. Our ANOCE is established under the
same causal structure of Chakrabortty et al. (2018) but without sparsity and normality assumptions.
Wide literature causal discovery can be summarized in three classes. The first type focuses on local conditional independence tests to find a causal skeleton and then determine the orientation of edges, such the well- PC algorithm (Spirtes et al. 2000; Kalisch & Bühlmann, 2007). However, testing the conditional independence of continuous variables is not easy (Shah & Peters 2018). The second class specifies properly functional causal models with additional assumptions on data distribution, including the ICA-LiNGAM (Shimizu et al.) 2006) and the causal additive model (CAM) (Bühlmann et al. 22144). The last class, the score-based method, includes the greedy cquivalence search (GES) (Chickering 2002) and the fast GES (IGES) (Ramsey et al. (2017) that use for example Bayesian scores in searching space of causal models. Recently, [Zheng et al. (2018) opened up another track of score- score-based methods by constructing an optimization with an acyclicity constraint under the LSEM, i.e. the NOTEARS, A follow-up work using a VAE parameterized by graph neural network that generalizes LSEM was proposed in Yu et al. (2019) with a more computational friendly constraint, namely DAG-GNN. Also see ZZuu & Chen (2019) and Lachapelle ct al. (019) for other cutting-edge structural learning methods.
The improvement of our ANOCE-CVAE over the state-of-the-arts is as follows. We consider a new constrained structural learning. by incorporating the background knowledge (the temporal causal relationship among variables) into the score- based algorithms. We formulated such prior information as the identification constraint and add it as the penalty term in the objective function for the causal discovery. In this paper, we typically extend the DAG-GNN for an illustration. Note that the proposed constraint is not limited to the DAG-GNN and can be easily extended to other score-based algorithms.
3

---

Published as a conference paper at ICLR 2021
3 ANALYSIS OF CAUSAL EFFECTS
3.1 STATISTICAL FRAMEWORK AND ASSUMPTIONS
Let A be the exposure/treatment. M [M M2,... M₂]" be mediators with dimension p. and Y be the outcome of interest. Suppose there exists a weighted DAG g (X, B) that characterizes the causal relationship among X [A, M Y]", where the dimension of X is = p 2. Let (A = .M m) be the potential outcome that would be observed after receiving treatment a and setting mediators as m, and M (AA = a) be the potential mediators that would be observed after receiving treatment a. As standard in the causal inference (Rosenbaum & Rubin 1983), we assume that there is no unmcasured confounder: (A1) the effect of the treatment on the outcome Y is unconfounded, i... (- = a, M m) 1 A, Va, m; (A2) the effect of the treatment A on the mediators M is unconfounded, i.e. M* = a) A, Va; (A3) the effect of the mediators M on the outcome Y is unconfounded given the treatment A, i.e. Y"(A In addition, as standard in the graphical causal discovery, we also make the Markov condition, the faithfulness condition, causal sufficiency assumption, and the linear structural equation model (LSEM) such that X = M. Y] characterized by the pair (G, c) is generated by
a, M m) T. M| Ya, m.
== X +,
where e is a random vector of jointly independent error variables.
Denote all directed paths in G that start with the exposure A and end with the outcome Y as set (T)Y(()) If there exists at least one directed path x² E {I (9)} such that the length of n" is larger than 2, we say there is an interaction among mediators, as shown in Figure otherwise, we call mediators are 'parallel' as shown in Figure la In this paper, we consider all possible causal structures with multiple mediators under assumptions (AI-A3).
We next give the total effect (TE), the natural direct effect that is not mediated by mediators (DE), and the natural indirect effect that is regulated by mediators (T)) defined Pearl et al. (2009).
Definition 3.1 Pearl al. 2009
TE = 1, [{Y] a)},
DE
IE = E{Y (d - - - a)))
where do(A= a) is a mathematical operator to simulate physical interventions that hold A constant as a while keeping the rest of the model unchanged, which corresponds to remove edges into and replace A by the constant a in G. Here, mn(a) is the value of M if setting do(A = a), and m m(++1) is the value of M if setting A + 1). Refer Pearl et al. (2009) for more details of 'do-operator'
Note that in the assumed linear model, the slope of the line is the same everywhere; for convenience and simplicity, we use a and a + 1 to present the change of the treatment in the definition.
3.2 NATURAL DIRECT AND INDIRECT EFFECT FOR INDIVIDUAL MEDIATORS
We first give the definition of the natural direct effect for individual mediator (DM). Definition 3.2
Natural direct effect for Mi:
DM, = [[E(M, = = + 1)} E{M, do(A = =a] a)}
(2)
[E(YYdo(AA- a, M, m 1,9, = o F[((do(A==))]]
where mi (a) is the value of M; when setting do A = a), Ω = M \ M; is the set of mediators except M,. and o. is the value of Q4 when setting do(A= a).
Remark 3.1 From Definition 3.2 the natural direct effect for M, is the product of the total effect of the treatment A on the mediator M; and the direct effect of the mediator Mg on the outcome Y. The
4

---

Published as a conference paper at ICLR 2021
second multiplier is in line with the classical meaning of "natural' in the causal inference literature PPaarl et al. 2009). Thus, the DM can be interpreted as the causal effect through a particular mediator from the treatment on the outcome that is not regulated by its descendent mediators.
The natural indirect effect for an individual mediator (IM) can be defined similarly.
Definition 3.3 Natural indirect effect for M,:
IM; [M 1)) E{M, -a))] =a)]]
E[Y]do(A== a, m; (a) + 1}} E(Y = m( + +, = o0
Remark 3.2 The second multiplier in the T TM; captures the indirect effect of a particular mediator on the outcome regulated by descendent mediators. We show the individual mediation effect (r)) in Chakrabortty et al. (2018) can be decomposed into the DM and the 1M in Section F in the appendix when the LSEM assumption holds, i.e., πソ- = DM; + MM, fo for i-th mediator.
Next, give explicit expressions of defined causal effects under the LSEM. Specifically, write the linear structural model [ under assumptions (A1-A3) as
A =
A M Y
0 pxf
0
A
T
+
OF
B
B M
M
+
6Mp
(3)
Y
L
Y
where y is scalar, B, and 0g px1 are p x 1 vectors, Ba is p x matrix, and c [CA, T [T. Here, presents the weight of the edge A + Y. the i-th element of a corresponds to the weight of the edge A -> M,, and the i-th element of /3 is the weight of the edge M; + Y. Note that by assumptions (AI-A3) we have the exposure A has no parents and the outcome Y has no descendants, equivalently, the first and the last column of B are all zeros (i.e., the first column and the last row B are all zeros). Notice the exposure can be presented by its own noise, i.e., A = <A, since A has no parents, any exposure (with arbitrary noise distribution) will satisfy the LSEM assumption.
Next, we obtain expressions of causal effects under the LSEM in the following theorem. The proof can be found in Section G. of the appendix.
Theorem 3.1 Under assumptions (AI-A3) and Model 1 we have:
1). the natural direct effect is DE 7 3 (Ip B ax where Ip is a x identity matrix; 3). the total effect of A on Y is TE = + (Ip -B -1;);
2). the natural indirect effect is IE = =3
4). the natural direct effect of M, on Y is DM; = Bir - where Br is the i-th element of B corresponding to the weight of M, + Y. and [(p - Bbr) -1 Jass is the i-th element of
())
(4p -BM)- aa as the total effect of A on M, i.e. E{Mi a - E do(A = a)}-
One may refer to section in the appendix for the invertibility of I, BT-- Also, a toy example is provided in section E to illsstrate how to manually compute the causal effects defined above. Note that there is no explicit expression of the IM due to the complex interaction among mediators, while we provide its theoretical form in Section G.2 with its numerical form in Section B
Remark 3.3
Based on the result 2) and 4) in Theorem 3.1 the IE can be presented as an additive form of DMs, as shown in Theorem 3.2. Thus, the proposed natural direct effect of individual mediators is valid in the sense that it exactly decomposes the indirect effect of the exposure on the outcome.
Theorem 3.2 Under assumptions (Al-A3) and Model [ the IE can be decomposed through DMs:
IE
DM.
3.3 ANALYSIS OF CAUSAL EFFECTS TABLE
Based on the result TE = DE + IE in Pearl et al. (2009) and Theorem 3.2 we summarize the defined causal effects and their relationship in Table for the analysis of causal effects (ANOCE).
5

---

Published as a conference paper at ICLR 2021
Firstly, the causal effect of A on Y has two sources, the direct effect from A and the indirect effect via p mediators M (M,.... Mp). Next, the direct source has the degree of freedom (d.f.) as 1, while the indirect source has f. as p from p mediators. Note the true d.f. of the indirect effect may be smaller than p. since A may not be regulated by all mediators. Then, the causal effect for the direct source is the DE and for the indirect source is the / F., where the F can be further decomposed into DMS and each component corresponds to the natural direct effect for a specific mediator. The last row in the table shows that the DE and the E compose the total effect ' E with d. as p 1.
Table 1: Table of Analysis of Causal Effects (ANOCE Table).
Source
Degree of freedom
Causal effects
Direct effect from A Indirect effect via M
DE IE
P
M
1
DM DM2
M2
1
:
Mp
1
DM,
Total
+ p
TE
4 CONSTRAINED STRUCTURAL LEARNING FOR ANOCE
We next estimate the weighted adjacency matrix B with our causal framework under the LSEM to calculate causal effects. To better capture the sampling distribution faithful to the DAG, we consider a deep generative model that generalizes the LSEM instead of using a regression that heavily relies on assumptions of noise (see more discussion in Section 2.2 Specifically, the LSEM | can be
rewritten as (Ip+2 B ) c, where 1p+2 is (p+2) (p + 2) identity matrix. Inversely, we have X (Ip+2 B )- e. Following the VAE architecture in Yu et al 2011)), treat the random error c as the independent latent variables to generate X, by two multilayer perceptrons as the encoder and the decoder, with weights denoted 8. We adopt their acyclicity constraint on B as,
hi (B)= tr| (Ip+2 I. Byb+
(p + 2) 0,
where tr(-) is the trace of a matrix, / is a hyperparameter that depends on an estimation of the largest eigenvalue of B, and denotes for the dlement-wise square.
Next, to incorporate the background knowledge of the temporal causal relationship among variables, we propose an identification constraint that indicates the topological order of the exposure and the outcome. As mentioned in Equation 3 under assumptions (A1-A3) the exposure A has no parents, ie. AA (g) Ø, and the outcome Y has no descendants, i.e. Y & (9). Or equivalently, we have the first column and the last row of B should equal to zero. Therefore, the matrix B must satisfy
p-2
p++2
(B) = Σ
Σ j=2
10 0,
(5)
where brg is the element of the matrix B in i-th and j-th column. The above constraint forces the topological order of the exposure as 1 while the outcome as p + 2, under which the DAG is searched within a restricted regime. The prior knowledge in in can can be generalized for any measured variable
and on the possible set of their parents, by connecting the topological order to the weighted matrix B.
Following Yu et al. (2019). the objective function is the evidence lower bound with two constraints:
min B,
p+2 E DxL (xu f( 1
(log.((
(6)
s.t. (B) = and bB)
where (-1-- is the Kullback- divergence, p(e) is the prior distribution of e, g(elX.) is the reconstructed empirical posterior distribution of e, and p( ((Xxle) is the likelihood function. Then, we have the loss function based on the augmented Lagrangian as
L B, 0, >1, A2) (), + X1hy (B) X2h2 (B) + c(,,(B)) + dB))

---

Published as a conference paper at ICLR 2021
(a) True. (b) v2 50
(c) 500
(d) True. (e) 50
(1) 500
(g) True. (h) n
50
ii 500
Honna 1990
(i) True.
(k) n (1) 1 500
(m) Y
(m) 72 500
Figure 2: The averaged estimated B under the ANOCE-CVAE: (a-c) correspond to Scenario 1
(d-f) correspond to Scenario 2; (g-i) correspond Scenario (i) is the shared true DAG for Scenario 4 and 4"; (k-l) are results under Scenario 4; (m-n) are results under Scenario 4".
where d and A2 are Lagrange multipliers, and and are penalty terms. To minimize the loss in and satisfy both h, (B) 0 and fı2 (B) 0, simultancously update Å and X2 and increase and to infinity, by modifying the basic technique in Yu et al. (209)). Here, the minimization can be solved using blackbox stochastic optimization such as 'Adam' in Kingma & Ba (2014). Denote the estimated matrix B from the above constrained structural learning. Under Theorem 3.1 we can estimate causal effects in the ANOCE table based on the learned B. We name the above algorithm as ANOCE-CVAE, with detail pseudocode provided in Section B
We incorporate the temporal causal relationship among variables when using an optimization approach to the causal discovery. Such constrained structural learning is not limited to the VAE framework and can be extended to any score-based algorithms. For instance, one can add constraint into the objective function in Zheng et al. (2018) or the reward in Zhu & Chen (019))
Remark 4.1
5 EXPERIMENTS
We conduct extensive simulation studies to investigate the proposed method on learning causal effects with multiple mediators, followed by comparison to the popular structural learning algorithms. The dataset and the code publicly available at https: /ttt::////ttttt com/ anw.e- AVACE
5.1 SIMULATION STUDIES
Scenarios are generated as follows. Scenario 3, we fix the dimension of M as p = 10 while increasing the complexity of the true graph to examine the sensitivity of our algorithm to sparsity. Specifically, Scenario is the simplest causal graph with only one edge (A - Y) shown in Figure 2a and Scenario 2 has fully connected graph with independent mediators (corresponding the parallel Scenario 3, we consider interacted mediators such allow p 30 with interacted mediators to examine the stability of our method under the high-dimensional setting. Here, the true DAGS in Scenarios and 4 are generated from the Erdös-Revyi (ER) model with expected degree as 2. Note that we consider fully identifiable models in the experiments so that it is meaningful to evaluate causal effects from the estimated graph. The synthetic datasets {A, M,Y} Y) are generated from Model | with Gaussian errors in Scenario 1-4. We also set A E {-1,1} } in Scenario 4 to show that our algorithm is capable to handle both discrete and continuous exposure, denoted as Scenario 4". The sample size is chosen from {50, 500} to be consistent with the scale of our real data. See more details of the data generation Section C. and the implementation in Section C.2 in the appendix.
case, i.e. BM = 0pxp) illustrated in Figure 2d
that BM + 0 pxxp demonstrated in Figure 2g For Scenario
The averaged estimated matrix T over 100 replications under the proposed ANOCE-CVAE is illustrated in Figure The numerical results are summarized in Table 2 (for Scenario 1 to 3) and

---

Published as a conference paper at ICLR 2021
(a) True.
(b) ANOCE
(c) PC
(d) ICA
(c) NOTEARS
(f) DAG-GNN
Figure 3: The averaged estimated T for Scenario 3 (ER2) with the graph threshold as 0.3.
Table 3 (for Scenario 4 and 4*) in the appendix, including the bias of the estimated TE, DE, IE, DM and IM for each mediator with their standard error. It can be observed that our proposed method could correctly identify most of the edges in the causal graph when n = 500 in almost all cases. Based on Table 2 aannd] 3 the estimated causal effects close the true values the sample size increases, indicating the good performance of our proposed method on identifying the causal effects regardless of the sparsity, the distribution of the exposure, and the dimension of mediators.
5.2 COMPARISON
We next compare our approach against the PC (Spirtes et al. 2000), the ICA-LiNGAM (Shimizu ct al. 2006), the NOTEARS (Zheng et al. 2018). and the DAG-GNN (Yu et al. 2019). Random graphs are generated from both the ER and the Scale-Free (SF) networks with the expected degree as 1, 2, and 4, denoted as Cases ERI, FR2, ER4, SF1, SF2, and SF4, respectively. To be consistent with Section 5.1 refer Scenario 3 (generated by the ER with the degree as 2) ER2, and set d = 12 (i.e. = 10) with Gaussian errors for other five cases, under n 500. Details of the data generation and the implementation of each method reported in Section C. and C.2 Here, we use graph threshold as 0.3 (commonly used in other methods) and 0.4 to prune the noise edges for fair comparison. The averaged estimated matrix B over 100 replications under different methods is shown in Figure with a graph threshold as 0.3 for Scenario (i.c., Case ER2) as an illustration. See other cases in Figures 5 to 15 in the appendix. All the numerical results of six cases are reported in Tables 4 and 5 in the appendix, including the false discovery rate, the true positive rate, and the structural Hamming distance. It is shown our algorithm performs the best among five methods in most cases, followed by the NOTEARS and the DAG-GNN. The comparison studies not only support the choice of the extension on the score-based algorithm (by comparing the results of the NOTEARS and the DAG-GNN with other methods), but also validate the improvement of our method over the DAG-GNN by introducing the background knowledge the causal discovery.
REAL DATA ANALYSIS: COVID- 19 OUTBREAK
From early Jan 2020 to late Feb 2020 (the Spring Festival period), COVID-19 spread to every province-level division of China, exacerbated by the Chinese new year migration and human to human transmission. The Chinese Government locked Hubei down on Jan 24th, which directly blocked infected people leaving from Hubei, and also indirectly control the spread of COVID-19. Here, for simplicity, we attribute the causal effects of the measures taken by cities outside Hubei to the original and main action of interest, i.e. Hubei lockdowns. Thus, the COVID-19 example satisfies the considered causality framework for studying the causal mediator effects.
We collect the data from the National Health Commission (NHC) of China and Baidu Qianxi for analysis. Specifically, let the exposure A as if Hubei is on lockdown, 0 for unlocked (before and on Jan 23rd), and for locked (on and after Jan 24th). We select 30 candidate cities outside Hubei that contain most potential infected people, as mediators M. The daily migration scale index (MSI) of each city is used as the value of each mediator, which is the migration magnitude of large groups of people from one geographical area to another (Chen et al. 2020) and is comparable among cities. Lastly, we use the daily increase rate of confirmed cases out of Hubei characterize the severity of the virus spreading with a one-weck delay (due to the diagnose and incubation period of COVID-19 LLuer et al. 2020 Y, Confirmed cases out of Hubeis Canlirmed cases oui of Hubei Here, the time starts Contirmed cases Hubei,+7
from Jan 12th to Feb 20th, 2020, since Jan 19th, 2020 the earliest date with an available number of
8

---

Published as a conference paper at ICLR 2021
15 Cargit
(a) The estimated weighted matrix for
(b) The estimated (in)direct effects for selected cities.
Figure 4: The analysis of the causal effects of 2020 Hubei lockdowns on reducing the COVID-19
spread in China regulated by Chinese major cities outside Hubei.
confirmed cases out of Hubei (to compute Y-=1 on Jan 12th), and after Feb 20th, 2020, the pandemic under control outside Hubei with the evidence of the work resumption in China. The final dataset yields total of 38 records. More details of data collection can be found in Section of the appendix.
The proposed algorithm is applied to the COVID-19 data with 100 replications by setting different random seeds in the neural network. The estimated weighted adjacency matrix is shown in Figure 4a with the detailed ANOCE table reported in Table 6 in the appendix. The total effect of 2020 Hubei lockdowns on the daily increase rate of confirmed cases outside Hubei as -0,497, where the direct effect -0,078 and the indirect effect is -0.419, In other words, by locking Hubei down, China successfully reduced 49.7% of the daily new cases outside Hubei; 84% of which is the indirect effect of Hubei, and the rest 16% owes
contributed via the reduced migration of cities (the mediators) to the direct effect of Hubei lockdowns since infected people were constrained in Hubei after the lockdown. Thus, the lockdown is effective in reducing the COVID-19 spread China.
The total indirect effect of the lockdown can be further broken down by cities' direct effects (DMs, corresponds to the intensity of transmission within a particular city). We compare cities' DMs with their associated indirect effects (,Ms, describes the secondary migration from a particular city to other places) in Figure 4b where a positive effect means spreading the virus while negative means control. Note that the selected 30 cities ordered by their cumulative MSI during the data period. 1). From Figure [ the majority of cities have a negative DM (colored in blue), which implies the infection within cities outside Hubei have been effectively controlled under the lockdown. 2). There are more cities with a positive IMs (red), which is in line with the intuition that the secondary migration among cities may exacerbate the pandemic. 3). The positive effects (red) are more likely located at the first 20 nodes, which corresponds to the cities with large MSI, while the last 10 cities with relatively small MSI are almost all blue. This accords with the migration peak among big cities during the Spring Festival period that aggravated the spread of the virus.
7 CONCLUSION
We conclude our paper with the following discussions. First, the proposed DM can be extended beyond the LSEM assumption.
generalized definition of the DM from a graphical perspective is given in Section F.2 of the appendix without the LSEM. Second, due to possibly unmeasured confounders our real data, such as cities' features and periodic effect, we may consider extending model with topological order that contains confounders for wider utility, such as forcing the topological order of & confounders as I to & followed by the exposure as + k. Third, our proposed identification constraint can be generalized to other background knowledge.
ACKNOWLEDGMENTS
The authors are grateful to the anonymous reviewers for valuable comments and suggestions. Rui Song's 's research is partially supported by grant from the National Science Foundation

---

Published as a conference paper at ICLR 2021
REFERENCES
Chen Avin, Ilya Shpitser, and Judea Pearl. Identifiability of path-specific effects. 2005.
Albert-LászIó Barabási and Réka Albert. Emergence of scaling in random networks. science, 286 (5439):509-512, 1999.
Peter Bühlmann, Jonas Peters, Jan Ernest, et al. Cam: Causal additive models, high-dimensional order search and penalized regression. The Annals of Statistics, 42(6))2526-2556, 2014.
David Card. The causal effect of education on earnings. In Handbook of labor economics, volume 3, pp. 1801-1863, Elsevier, 1999.
Abhishek Chakrabortty, Preetam Nandy, and Hongzhe Li. Inference for individual mediation effects and interventional effects in sparse high-dimensional causal graphical models. arXiv preprint arXiv: /809. 10652, 2018.
Huijie Chen, Ye Chen, Baijun Sun, Ping Wang, Lihai Wen, Zhiyong Lian, Ying Lu, Ying Qi, Shuo Zhao, Linlin Zhang, et al. Correlation between the migration scale index and the number of new confirmed novel coronavirus pneumonia cases in china. 2020.
David Maxwell Chickering. Optimal structure identification with greedy search. Journal of machine learning research, 3(Nov):507-5544 2002.
Miguel Angel Hernán. A definition of causal effect for epidemiological research. Journal of Epidemiology Community Health, 58(4):265-271, 2004.
Miguel Ángel Hernnn, Babette Brumback, and James M Robins. Marginal structural models to estimate the causal effect of zidovudine on the survival hiv-positive men. Epidemiology, pp. 561-570, 2000.
Markus Kalisch and Peter Bühlmann. Estimating high-dimensional directed acyclic graphs with the pe-algorithm. Journal of Machine Learning Research, 61Mar)66136636 2007.
Diederik P Kingma and Jimmy Ba. Adam: method for stochastic optimization. arXiv preprint arXiv:1412.6980, 2014.
Sébastien Lachapelle, Philippe Brouillard, Tristan Deleu, and Simon Lacoste-Julien. Gradient-based 2019.
neural dag learning. arXiv preprint
Stephen A Lauer, Kyra H Grantz, Qifang Bi, Forrest K Jones, Qulu Zheng, Hannah R Meredith, Andrew S Azman, Nicholas G Reich, and Justin Lessler The incubation period of coronavirus disease 2019 (covid-19) from publicly reported confirmed cases: estimation and application. Annals of internal medicine, 2020.
Marloes H Maathuis, Markus Kalisch, Peter Bühlmann, al. Estimating high-dimensione interven- tion effects from observational data. The Annals of Statistics, 37(6)):3133-11 3164,1 2009.
Preetam Nandy, Marloes H Maathuis, Thomas S Richardson, et al. Estimating the effect of joint interventions from observational data in sparse high nigh-dimensional settings. The Annals of Statistics, 45(2):647-674, 2017.
Ugo Panizza and Andrea F Presbitero. Public debt and economic growth: is there a causal effect? Journal of Macroeconomics, 41:21-41, 2014.
Adam Paszke, Sam Gross, Soumith Chintala, Gregory Chanan, Edward Yang, Zachary Detooo, Zeming Lin, Alban Desmaison, Luca Antiga, and Adam Lerer. Automatic differentiation in pytorch. 2017.
Judea Pearl et al. Causal inference in statistics: An overview. Statistics surveys, 3:96-146, 2009.
Joseph Ramsey, Madelyn Glymour, Ruben Sanchez-Romen and Clark Glymour. A million variables and more: the fast greedy equivalence search algorithm for learning high -ish-iinnssoonl graphical causal models, with an application to functional magnetic resonance images. International journal of data science and analytics, 321):121- 2017.
10

---

Published as a conference paper at ICLR 2021
Paul R Rosenbaum and Donald B Rubin. The central role of the propensity score in observational studies for causal effects. Biometrika, 70(1))41-55, 1983.
Rajen D Shah and Jonas Peters. The hardness of conditional independence testing and the generalised covariance measure. arXiv preprint arXiv:1804.00220. 2018.
Shohei Shimizu, Patrik Hoyer, Aapo Hyvärinen, and Antti Kerminen. linear non-gaussian acyclic model for causal discovery. Journal of Machine Learning Research, ((0ct):2003-2030 2006.
Pater Spirtes, Clark Glymour, Richard Scheines, Stuart Kauffman, Valerio Aimale, and Frank Wimberly. Constructing bayesian network models of gene expression networks from microarray data. 2000.
Tyler VanderWeele and Stijn Vanstcelandt. Mediation analysis with multiple mediators. Epidemio- logic methods, 2(1):95-115, 2014.
Stijn Vansteelandt and Rhian M Daniel. Interventional effects for mediation analysis with multiple mediators. Epidemiology (Cambridge, Mass.), 28(2):258, 2017.
Sewall Wright. Correlation and causation. Journal of agriculuurll research, 20(7):557-585, 1921.
Yue Yu, Jie Chen, Tian Gao, and Mo Yu. Dag-gnn: Dag structure learning with graph neural networks. arXiv preprint 1rXiv:190..1000 2019.
Xun Zheng, Bryon Aragam, Pradeep K Ravikumar, and Eric P Xing. Dags with no tears: Continuous optimization for structure learning. In Advances in Neural Information Processing Systems, pp. 9483, 2018
Fei Zhou, Ting Yu, Ronghui Du, Guohui Fan, Ying Liu, Zhibo Liu, Jie Xiang, Yeming Wang, Bin Song, Xiaoying Gu, et al. Clinical course and risk factors for mortality of adult inpatients with covid-19 in wuhan, china: retrospective cohort study. The Lancet, 2020.
Shengyu Zhu and Zhitang Chen. Causal discovery with reinforcement learning. arXiv preprint 14X1471906.04477, 2019.

---

Published as a conference paper at ICLR 2021
A ADDITIONAL GRAPH TERMINOLOGY
Given the node set, the weighted DAG can be uniquely determined by its weighted adjacency matrix, i.e. there is one-to-one transformation between & and B. Suppose the graph nodes X in G are sorted in its topological order (corresponding to elementary transformation of the matrix), then the matrix B is strictly upper triangular with the diagonal elements as O. Therefore, for identity matrix with the same dimension as B, I - B is invertible since all its diagonal elements are (positive).
B ALGORIHTM: ANOCE-CVAE
The first part of the ANOCE-CVAE algorithm is on learning causal DAG from the observational data in the constrained space, by minimizing the loss function in Equation 7 using blackbox stochastic optimization solvers. Here, to minimize the loss in 7 and satisfy both h1 (B) = 0 and h2 (B) O. we simultaneously update a and X2 and increase and d to infinity, by modifying the basic technique in Yu et al. (2019), corresponding to Part Il. and Part One. in Algorithm [ The second part is to estimate causal effects in the ANOCE table from the learned causal structure, based on the results from Theorem 3.1 Here, we numerically calculate the natural indirect effect for mediator / M
based on Corollary F. in step IV of the second part.
Algorithm 1 Analysis of Causal Effects via Constrained VAE (ANOCE-CVAE)
Global: Dataset X {A, M,, }. sample size n, dimension of mediators p. max iteration K, number of epoch H, original learning rate ro. tolerance of constraint to zero 8, parameter update bound U, tuning parameters p and w, and penalty terms c and d: Local: mean and standard variance of flg and oe, mean and standard variance of X ax and ox. weights in multilayer perceptrons of encoder and decoder 0 {W(1) W(2), W3)), w(4),, Lagrange multipliers a and X2 penalty terms e and d. (p + 2) x (p + 2) matrix B, Loss function L, old and new values of the first constraint A and h. old and new values of the second constraint hahd and hn and learning rate r; Output: estimated matrix B, total effect TE, natural direct and indirect effect DE and IE, and natural direct and indirect effect for mediator DM and IM.
Part One: Generate matrix B via Constrained Variational Auto-Encoder;
I. Initialization: A1 f 0; )2 + 0; c < 1: d t 1 +- ro; B 0l) h1th +- 0o; high <- 50;
II. For step k, l 1.... K:
A. While e x < U
a). For epoch i, 1,..., H::
1. Build Encoder (e), +- BTMMIP.X W(), W
2. Build Decoder (xx, t MLP(([p42 B - (33)), W(1); Calculate values of constraints h:pww + -i((B) and hgo t (B2(
and the loss function L t Le, W(), W(2), w()), w(t1, Ar. d2);
4. Use backward to update parameters (B, W(1), W(2), W3), w(4)):
5. Update learning rate r;
b). If h:pw > >ph and hg-w > phg c c c x wi d t- d x wi
Elseif ht > phgli and new phold, c ct x w;
Elseif hiw < phgl and hnew phgld: d x- d w;
Else: Break;
B. hegld <- h.pew; high hgww; > + A1 x hppw X22 + X2 x hp
C. If hyow < 8 and keww < 8: Break;
III. Output B < B;
12

---

Published as a conference paper at ICLR 2021
Algorithm 2 ANOCE-CVAE (cont.)
Part Two: Estimate causal effects in ANOCE based on matrix B:
I. According to Equation
Get y as the direct effect DD;
B. Get ā as the effect of on M, 3, and the inside matrix Bas; l that represents the causal effect of A on M; III. Get { that represents the total natural indirect effect IE; For each mediator M,, 1,... pi Define the natural direct effect for M, as DM a[d][]]
II. Get (4p ( BM-
IV. Get the natural indirect effect for mediator: For each mediator M,, 1,... p.
A. Delete M, from the matrix B and get Bi;
B. Repeat step II. with reduced matrix B; and get 3 and <:;
C. Calculate the effect difference the total mediation effect B 3 Tm
D. Define the natural indirect effect for M, as IM[i { - 3 {) DM[i]:
Define the total effect TE= 7 + 3.
C ADDITIONAL SIMULATION STUDIES
In this section, give more details simulation studies to investigate the finite sample performance of the proposed method for learning causal effects with multiple mediators, in comparison to the popular causal discovery methods, including the PC, the ICA-LINGAM, the NOTEARS, and the DAG-GNN. The computing infrastructure used is virtual machine in the compute engine of Google Cloud Platform with 8 processor cores and 32GB memory. The average runtime for each result around 1 to 2 hours.
C.1 DATA GENERATION
We first generate a random DAG from the Frdös-Reóyi (ER) or the Scale-Free (SF) network (Barabási & Albert 1999) with expected node degree. Then, remove all in-edges (from precedent nodes) of the first node as and remove all out-edges (from descendent nodes) the last node as Y. and randomly thus, the remaining nodes are the mediators M. Edges in DAGs for all scenarios assigned with weights (w E 1} with equal probability) to obtain the weighted adjacency matrix B. Specifically, the DAGs Scenarios and are generated from the Erdös-Reñyi (ER) model with an expected degree as 2, where we set number of nodes d = 12 (i.e. p 10) in Scenario and d 32 (i.e. p = 30) in Scenario 4. Note that we consider fully identifiable models Section so that it is meaningful to evaluate causal effects from the estimated graph. In Section 5.2 we repeat the above generation procedure with d the Scale-Free (SF) networks with the expected degree as 1, 2, and 4, denoted as Cases ER1, ER2, ER4, SF1, SF2, and SF4, respectively. Here, be consistent with Section [.1] we refer Scenario (generated by the ER with the degree 2) Case ER2.
to generate the true graph from both the ER and
The synthetic datasets (A, M, Y} are generated from Model 1 where the error variables in e = S
[eA, (a. <y] E binary exposure in Scenario 4' Here, we add a baseline of 1.0 on the outcome Y. Note that the Gaussian exposure in Scenario and the binary exposure in Scenario 4' have the same mean and noise and thus their results are comparable.
independently follow normal distribution with mean 0 and noise 0.5 except for the
C.2 IMPLEMENTATION DETAILS
We detail the implementation for the proposed ANOCE-CVAL and comparison partners as follows:
ANOCE-CVAE: The ANOCE-CVAL is implemented based on PyTorch (Paszke et al. (017). using Adam (Kingma & Ba 2014) minimize the loss function in Equation] 7 We set the batch
13

---

Published as a conference paper at ICLR 2021
size as 25 for = 50 and 100 for 500 with hidden nodes as p², the initial learning rate as 0.003 with an update rule as r + r/(log(c) + log(d) + 0.01} where e and are penalty terms for two constraints, and the parameter update bound as U 10 for all settings. Following the recommendation of Yu et al. (2019). we find that their tuned parameters work well in our settings, and we adopt the Huber-norm regularization of B for better convergence. Here, the variational posterior and the likelihood are parameterized as Gaussian with unit noise to approximate the underlying true model. The code is publicly available an anonymous repository at https: //tp://gihttbt com/ anwaannccvaaae Anoce--vae/ANCCE-CVA. ANOCE-CVCE
0.25 and w = 10 also
PC (Spirtes al., 2000 We set the Fisher-z test for the PC algorithm with the p-value as 0.01 for all settings. The implementation is available through the py-causal package at IJ written in highly optimized Java codes. Also pwwpkkkcccddypccaaaaaaaaa 1/blallllloo ob/ /eeee 1ompent...t//// ipp.
https: / thub com wttpc////wwcccccccccooooeooe see http:////.....aaaa//-0//00///2/0
example/py-ccuse
ICA-LINGAM (Shimizu et al. 2006): The ICA-LiNGAM assumes linear non-Gaussian additive model to recover the weighted adjacency matrix. We implement the ICA-LiNGAM with default
hyper-parameters through the lingam package for all settings. See their repository at https
7/github. com/ wow.cdtt//minnnn
NOTEARS (Zhngg et al. 2018): The NOTEARS estimates the weighted adjacency matrix by formulating the optimization with an acyclicity constraint. The implementation is available at their repository at https: / gigggitthhh thub. com/ xunzheng noteaze We set the loss function as the least square error with the L1 regularization. We find the NOTEARS is sensitive to the choice of the Li regularization in our settings. For a fair comparison, we set the Lı penalty parameter 0.03 (instead of the default 0. for all settings, which achieves an overall good performance in most cases. Note the author modified their acyclicity constraint in their codes to be the one used in Yu al. (2019) (i.e. Equation 4.. We also use the same acyclicity constraint for NOTEARS, DAGGNN, and our method for fair comparison. Other ayper-parameten are set as default in their repository.
DAG-GNN (Yu et al. 2019): The DAG-GNN incorporates the variational auto-encoder into causal discovery with a modified smooth characterization on acyclicity in the evidence lower bound as the loss function. Codes are available at their repository at https:: / htt/gtttbu..///// com/ fishmoon / GNN based on PyTorch (Paszke et al., (017)). We set the same hyper- parameters used in our ANOCE-CVAE for fair comparison. Specifically, we use Adam (Kingma & Ba 2014) to minimize the loss function, and set the batch size as 25 for n = 50 and 100 for 500 with hidden nodes as p². The initial learning rate is set 0.003 with update rule as r + r (log(ca) + 0.01) where og is penalty term for the acyclicity constraint. The rest settings are the same as the default in their codes.
In the comparison studies (see Section 5.2 and C.4 we use a uniform graph threshold as 0.3 (commonly used in current literature) for all algorithms to prune the noise edges for fair comparison. In addition, we also provide the results under the graph threshold as 0.4 for additional comparison.
C.3 ADDITIONAL RESULTS OF ANOCE-CVAE
In this section, we provide additional simulation results for the ANOCE-CVAE. Following Section [.I the numerical results are summarized in Table (for Scenario 3) and Table 3 (for Scenario and 4"), including the bias of the estimated TE, DE, IE, DM and IM for each mediator with their standard error, over 100 replications. Note that due to limited space, we save the numerical results of
the 7M in Table 3
From the results in Table 2, it is clear that the estimated TE, DE, IF, DM and TM for each mediator are close to the true values as the sample size increases in Scenario I to 3, which indicates the good performance of our proposed method on identifying the causal effects regardless of the sparsity. With the expected node degree increasing, one can observe a slightly larger bias and standard error of the estimated causal effects as expected, as shown in Table 2 Based on Table 3 the results of Scenario 4 and Scenario 4" merely identical under different sample sizes, indicating our proposed method can handle either discrete or continuous exposure. In addition, by comparing the results of Scenario 3 and 4 where fix the expected node degree as 2, one observe a slightly larger bias
14

---

Published as a conference paper at ICLR 2021
of the estimated causal effects but of similar small scale, the dimension of mediators p increases, which implies the stability of our method under the high- nigh-dimensional setting.
Table 2: The bias and standard error of the estimated causal effects for Scenario 1 to 3 (p = 10).
Scen. (simplest)
Scen. 2 (parallet)
Scen. 3 (interacted)
Results
n =
50
n = 500
n = 50
n = 500
n = 50
n = 500
TE
0.19 (0.02)
0.01 (0.00)
0.66 (0.05)
0.03 (0.02)
0.15 (0.05)
0.05 (0.03)
DE
0.23 (0.02)
0.01 (0.00)
0.78 (0.02)
0.12 (0.06)
0.78 (0.03)
0.15 (0.06)
IE
0.04 (0.01) 0.00 (0.00)
1.44 (0.04)
0.08 (0.06)
0.63 (0.04)
0.20 (0.04)
DM1
0.01 (0.00) 0.00 (0.00) 0.69 (0.03) 0.02 (0.01) 0.48 (0.04)
0.06 (0.02)
DM2
0.00 (0.00) 0.00 (0.00) 0.43 (0.03) 0.01 (0.01) 0.00 (0.00) 0.00 (0.00)
DM3
0.00 (0.00) 0.00 (0.00) 0.65 (0.02) 0.01 (0.02) 0.00 (0.00) 0.00 (0.00)
DM4
0.00 (0.00) 0.00 (0.00) 0.65 (0.03) 0.02 (0.01) 0.20 (0.02)
0.06 (0.01)
DMs
0.01 (0.00) 0.00 (0.00) 0.68 (0.03)
0.02 (0.01)
0.44 (0.02)
0.10 (0.01)
DM6
0.01 (0.00) 0.00 (0.00)
0.67 (0.02) 0.02 (0.02) 0.27 (0.04)
0.03 (0.01)
DM7
0.00 (0.00) 0.00 (0.00)
0.62 (0.03) 0.03 (0.01) 0.67 (0.03)
0.16 (0.03)
DMs
0.00
0.00 (0.00) 0.66 (0.03) 0.01 (0.01) 0.59 (0.03)
0.14 (0.03)
DMg
0.00 (0.00)
0.00 (0.00) 0.65 (0.02)
0.01 (0.01) 0.00 (0.00) 0.00 (0.00)
DM1o
0.01 (0.00) 0.00 (0.00) 0.66 (0.03) 0.01 (0.01) 0.10 (0.02) 0.05 (0.01)
IM1
0.00 (0.00) 0.00 (0.00) 0.08 (0.01) 0.08 (0.01) 0.76 (0.03)
0.16 (0.04)
TM2
0.00 (0.00) 0.00 (0.00) 0.09 (0.01) 0.01 (0.01) 0.00 (0.00)
0.00 (0.00)
TM3
0.00 (0.00) 0.00 (0.00)
0.03 (0.01) 0.03 (0.01)
0.00 (0.00) 0.00 (0.00)
1MA
0.00 (0.00) 0.00 (0.00)
0.03 (0.01) 0.05 (0.01)
0.16 (0.01) 0.03 (0.01)
IMG
0.00 (0.00) 0.00 (0.00)
0.06 (0.01) 0.12 (0.01)
0.02 (0.01)
0.03 (0.00)
IM6
0.00 (0.00) 0.00 (0.00) 0.03 (0.01) 0.04 (0.01)
0.66 (0.02) 0.09 (0.02)
1M7
0.00 (0.00) 0.00 (0.00)
0.09 (0.01)
0.00 (0.01)
0.05 (0.02) 0.04 (0.01)
IMs
0.00 (0.00)
0.00 (0.00)
0.03 (0.01)
0.02 (0.01)
0.11(0.02) 0.06 (0.01)
IMg
0.00 (0.00)
0.00 (0.00)
0.03 (0.01)
0.01 (0.01)
0.00 (0.00) 0.00 (0.00)
IM
0.00 (0.00)
0.00 (0.00)
0.07 (0.01) 0.08 (0.01) 0.03 (0.01) 0.00 (0.00)
15

---

Published as a conference paper at ICLR 2021
Table 3: The bias and standard error of the estimated causal effects for Scenario 4 and 4" (p 30).
Scen. 4 (continuous)
Scen. 4" (binary)
Results
71 50
500
50
500
TE
0.29 (0.04)
0.20 (0.02)
0.09 (0.03)
0.05 (0.01)
DE
0.90 (0.01)
0.28 (0.02)
0.85 (0.01)
0.25 (0.02)
IE
1.19 (0.04)
0.48 (0.03)
0.94 (0.03)
0.30 (0.02)
DM1
0.12 (0.02)
0.08 (0.01)
0.18 (0.02)
0.07 (0.01)
DM2
0.61 (0.02)
0.02 (0.01)
0.45 (0.02)
0.02 (0.01)
DM3
0.00 (0.01)
0.03 (0.01) 0.01 (0.01)
0.02 (0.01)
DM4
0.31 (0.03)
0.25 (0.02) 0.09 (0.02)
0.07 (0.01)
DM;
0.01 (0.01) 0.07 (0.01) 0.05 (0.03)
0.10 (0.02)
DM,
0.00 (0.00)
0.00 (0.00) 0.00 (0.00)
0.00 (0.00)
DM7
0.37 (0.02) 0.19 (0.02) 0.60 (0.04)
0.18 (0.02)
DMs
0.02 (0.01) 0.00 (0.01) 0.05 (0.01)
0.02 (0.02)
DMg
0.00 (0.00) 0.00 (0.00) 0.00 (0.00)
0.00 (0.00)
DM10
0.02 (0.01) 0.09 (0.01) 0.14 (0.01)
0.15 (0.01)
DM1
0.00 (0.01)
0.00 (0.00) 0.01 (0.01)
0.00 (0.00)
DM12
0.02 (0.00)
0.01 (0.01)
0.06 (0.01)
0.01 (0.01)
DMaa
0.25 (0.02)
0.11 (0.02)
0.51 (0.04)
0.16 (0.03)
DM4
0.02 (0.01)
0.10 (0.01)
0.09 (0.02)
0.16 (0.01)
DM15
0.00 (0.00)
0.00 (0.00)
0.00 (0.00)
0.00 (0.00)
DM1
0.09 (0.01)
0.10 (0.02)
0.18 (0.03)
0.13 (0.02)
DM17
0.03 (0.01) 0.07 (0.02) 0.04 (0.03)
0.06 (0.03)
DM1s
0.01 (0.00) 0.00 (0.00) 0.00 (0.00)
0.00 (0.00)
DM19
0.02 (0.00) 0.01 (0.00) 0.00 (0.00) 0.00 (0.00)
DM20
0.01 (0.00) 0.00 (0.00) 0.08 (0.01)
0.01 (0.01)
DM21
0.00 (0.00) 0.00 (0.00) 0.00 (0.00) 0.00 (0.00)
DM22
0.00 (0.00) 0.00 (0.00) 0.00 (0.00) 0.00 (0.00)
DM23
0.00 (0.00)
0.00 (0.00) 0.00 (0.00) 0.00 (0.00)
DM24
0.00 (0.00) 0.00 (0.00) 0.00 (0.00) 0.00 (0.00)
DM25
0.00 (0.00)
0.00 (0.00) 0.00 (0.00) 0.00 (0.00)
DM26
0.00 (0.00)
0.00 (0.00) 0.00 (0.00) 0.00 (0.00)
DM27
0.00 (0.00)
0.00 (0.00) 0.00 (0.00)
0.00 (0.00)
DM28
0.01 (0.01)
0.02 (0.01)
0.01 (0.01)
0.01 (0.00)
DM
0.00 (0.00)
0.00 (0.00)
0.00 (0.00)
0.00 (0.00)
DM30
0.00 (0.00)
0.00 (0.00) 0.00 (0.00)
0.00 (0.00)
16

---

Published as a conference paper at ICLR 2021
C.4 ADDITIONAL COMPARISON STUDIES
(a) True.
(b) ANOCE
(c) PC
(d) ICA
(c) NOTEARS
(f) DAG-GNN
Figure 5: The averaged estimated B for Case ERI under different methods with threshold 0.3.
(a) True.
(b) ANOCE
(c) PC
(d) ICA
(c) NOTEARS
DAG-GNN
Figure 6: The averaged estimated B for Case ER4 under different methods with threshold 0.3.
(a) True.
(b) ANOCE
(c) PC
(d) ICA
(c) NOTEARS
(1 DAG-GNN
Figure 7: The averaged estimated BT for Case SF1 under different methods with threshold 0.3.
(a) True.
(b) ANOCE
(c) PC
(d) ICA
(e) NOTEARS
(1 DAG-GNN
Figure 8: The averaged estimated B for Case SF2 under different methods with threshold 0.4.
(a) True.
(b) ANOCE
(c) PC
(d) ICA
(c) NOTEARS
(1 DAG-GNN
Figure 9: The averaged estimated B for Case SF4 under different methods with threshold 0.4.
This section provides more results on comparison studies against the existing methods. The data generation and implementation details are provided in Section C.1 and C.2 Following Section 5.2
17

---

Published as a conference paper at ICLR 2021
Table The empirical comparison results under different causal discovery methods with the graph threshold as 0.3.
Methods
Case
ERI
ER2
ER4
SF1
SF2
SF4
FDR
0.00 (0.14)
0.00 (0.09)
0.14 (0.03)
0.21 (0.06)
0.17 (0.05)
0.13 (0.06)
ANOCE-CVAE
TPR
0.50 (0.00)
1.00 (0.08)
0.93 (0.05)
1.00 (0.03)
0.96 (0.05)
0.81 (0.06)
SHD
5.00 (1.01)
0.00 (3.57)
8.00 (2.44)
3.00 (1.49)
6.00 (2.15)
10.00 (2.65)
FDR
0.00 (0.10)
0.50 (0.04)
0.23 (0.05)
0.00 (0.01)
0.29 (0.04)
0.27 (0.05)
PC
TPR
0.40 (0.01)
0.26 (0.01)
0.41 (0.04)
1.00 (0.00)
0.46 (0.02)
0.34 (0.02)
SHD
6.00 (0.55)
19.00 (0.51)
29 29.00 (2.01)
0.00 (0.17)
19.00 (1.19)
25,00 (1.18)
FDR
0.00 (0.18)
0.08 (0.16)
0.15 (0.12)
0.00 (0.14)
0.00 (0.18)
0.00 (0.14)
ICA-LiNGAM
TPR
0.40 (0.12)
0.52 (0.14)
0.41 (0.10)
0.64 (0.17)
0.50 (0.17)
0.47 (0.13)
SHID
6.00 (1.50)
12.00 (4.08)
26.00 (4.67)
4,00 (2.98)
13.00 (6.29)
17.00 (5.36)
FDR
0.00 (0.02)
0.00 (0.06)
0.04 (0.06)
0.00 (0.00)
0.00 (0.07)
0.04 (0.03)
NOTEARS
TPR
0.50 (0.00)
0.78 (0.09)
0.63 (0.08)
1.00 (0.00)
0.58 (0.08)
0.72 (0.07)
SHD
5,00 (0.10)
5.00 (2.29)
15.00 (3.90)
0.00 (0.00)
11.00 (2.82)
9.00 (2.59)
FDR
0.29 (0.07)
0.15 (0.06)
0.13 (0.04)
0.29 (0.07)
0.13 (0.06)
0.11 (0.05)
DAG-GNN
TPR
0.50 (0.00)
0.74 (0.06)
0.80 (0.06)
0.91 (0.06)
0.77 (0.07)
0.75 (0.05)
SHD
7.00 (0.49)
9.00 (2.22)
12.00 (3.08) 5.00 (1.66)
9.00 (2.23)
10.00 (1.71)
and C.2. we use a graph thresholds as 0.3 (commonly used in current literature) or 0.4 (for additional comparison) all algorithms prune the noise edges for fair comparison.
The estimated graphs (after pruning) are evaluated by three metrics: the false discovery rate (FDR), the true positive rate (TPR), and the structural Hamming distance (SHD, the smallest number of edge additions, deletions, and reversals to convert the estimated graph into the true DAG). Here, the SHD takes into account both false positives and negatives and lower SHD indicates better estimate of
the causal graph. The FDR, TPR, and SHD of the averaged estimated matrix B with their standard deviation over 100 replications reported in Table for the graph threshold as 0.3 and in Table for the graph threshold as 0.4, under different methods for all six cases with sample size n = 500.
18

---

Published as a conference paper at ICLR 2021
Table The empirical comparison results under different causal discovery methods with the graph threshold as 0.4.
Methods
Case
ER4
SF1
SF2
SF4
FDR
0.00 (0.13)
0.00 (0.08)
0.07 (0.04)
0.08 (0.05)
0.08 (0.04)
0.07 (0.05)
ANOCE-CVAE
TPR
0.50 (0.00)
1.00 (0.10)
0.90 (0.05)
1.00 (0.03)
0.92 (0.05)
0.78 (0.07)
SHD
5.00 (0.92)
0.00 (3.52)
6.00 (2.91)
1.00 (0.90)
4.00 (1.89)
9.00 (2.62)
FDR
0.00 (0.10)
0.40 (0.04)
0.20 (0.05)
0.00 (0.01)
0.31 (0.04)
0.27 (0.05)
PC
TPR
0.40 (0.01)
0.26 (0.01)
0.39 (0.04)
1.00 (0.00)
0.42 (0.02)
0.34 (0.02)
SHD
6.00 (0.55)
19.00 (0.51)
29 29.00 (2.01)
0.00 (0.17)
20.00 (1.19)
25,00 (1.18)
FDR
0.00 (0.30)
0.00 (0.15)
0.08 (0.12)
0.00 (0.15)
0.00 (0.19)
0.00 (0.15)
ICA-LiNGAM
TPR
0.20 14)
0.39 (0.14)
0.29 (0.09)
0.45 (0.16)
0.27 (0.17)
0.22 (0.13)
SHID
8.00 (1.72)
14.00 (3.90)
30.00 (4.28)
6.00 (2.51)
19.00 (5.48)
25.00 (4.99)
FDR
0.00 (0.00)
0.00 (0.04)
0.04 (0.05)
0.00 (0.00)
0.00 (0.06)
0.05 (0.03)
NOTEARS
TPR
0.50 (0.00)
0.65 (0.08)
0.59 (0.07)
1.00 (0.00)
0.58 (0.09)
0.56 (0.08)
SHD
5.00 (0.00)
8.00 (1.86)
17.00 (3.54)
0.00 (0.00)
11.00 (2.83)
14.00 (2.72)
FDR
0.29 (0.09)
0.15 (0.05)
0.11 (0.05)
0.17 (0.05)
0.06 (0.06)
0.09 (0.03)
DAG-GNN
TPR
0.50 (0.00)
0.74 (0.07)
0.78 (0.07)
0.91 (0.07)
0.65 (0.08)
0.66 (0.06)
SHD
7.00 (0.67)
9.00 (1.64)
12.00 (3.36)
3.00 (1.08)
10.00 (2.32)
(1.71)
(a) True
(b) ANOCE
(c) PC
(d) ICA
(c) NOTEARS
( DAG-GNN
Figure 10: The averaged estimated B for Case ERI under different methods with threshold 0.4.
(a) True.
(b) ANOCE
(c) PC
(d) ICA
(e) NOTEARS
DAG-GNN
Figure 11: The averaged estimated B for Case ER2 under different methods with threshold 0.4.
(a) True.
(b) ANOCE
(c) PC
(d) ICA
(e) NOTEARS
(O DAG-GNN
Figure 12: The averaged estimated B for Case ER4 under different methods with threshold 0.4.
Besides Figure 3 for Case ER2 with graph threshold as 0.3 shown in the main text, we also illustrate the averaged estimated matrix B. over 100 replications under different methods for Case ERI, ER4,
19

---

Published as a conference paper at ICLR 2021
(a) True.
(b) ANOCE
(c) PC
(d) ICA
(c) NOTEARS
(1) DAG-GNN
Figure 13: The averaged estimated B for Case SFI under different methods with threshold 0.4.
(a) Truc.
(b) ANOCE
(c) PC
(d) ICA
(c) NOTEARS
(f) DAG-GNN
Figure 14: The averaged estimated B for Case SF2 under different methods with threshold 0.4.
2.
-
(a) True.
(b) ANOCE
(c) PC
(d) ICA
(c) NOTEARS
( DAG-GNN
Figure 15: The averaged estimated B for Case SF4 under different methods with threshold 0.4.
SFI, SF2, and SF4 with graph threshold as 0.3 in Figure 5 ERI. ER2, ER4, SF1, SF2, and SF4 with a graph threshold as 0.4 in Figure 10 respectively, under
6
9 respectively, and for Case
[2
I3
14 15
500.
From Table and 5 it is clear that our algorithm performs the best among the five methods in most cases, followed by the other two score-bas based methods, i.e. the NOTEARS and the DAG-GNN. While the traditional methods (the PC and the ICA- ICA-LiNGAM) perform the worst with large SHD and small TPR. This finding supports the choice of the extension on the score-based method. Moreover, by comparing our performance with the DAG-GNN, one can observe substantial gain in terms of the SHD and the TPR in most cases with comparable FDR. This validates the improvement of our method over the DAG-GNN by introducing the background knowledge in the causal discovery. Another supports are illustrated in Figure 3
and 15 for different settings, where the averaged estimated matrix B under the ANOCE-CVAE, is approximately the same as the ground true graph B when n 500. However, the PC and the ICA (CA-LiNGAM can hardly recognize the true causal pattern. In addition. methods have a slightly better performance in terms of FDR and SHD while a slightly worse performance terms of TPR under the graph threshold 0.4, in comparison to the results under the graph threshold as 0.3.
g
12 13
D ADDITIONAL REAL DATA RESULTS
In this section, we provide additional real data analysis on the COVID-19.
D.1 DATA COLLECTION
To better characterize the causality the virus spreading under the Hubei lockdowns in China, assume: Hubei the centre of the COVID-19 outbreak in China (Zhuu et al. 2020); 2)
20

---

Published as a conference paper at ICLR 2021
Hubel Lockdown
Other Cities Migration
Virus
Spread
Figure 16: A direct acyclic graph illustrating the temporal causal relationship of the spreading network of the coronavirus outbreak under 2020 Hubei lockdowns China.
Marth tot. 2000
(a) The cumulative confirmed cases of selected citics.
Chunge
Sharghal
acropan
Chompinn
Changet
uumiing
Nanchung
shron
subou
Turtn
Anding
-1.0
-0.5
0.5
0.5
1.0
(b) The partial heat map of correlation matrix between cities and the treatment A the outcome Y.
Figure 17: Additional figures for the real data analysis of the COVID-19 spread.
the decreased migration outside Hubei was largely stimulated by the lockdown; individual who departed to one destination would not return to the original departure due to travel restrictions in temporal causal relationship to describe
China. Under the above assumptions, it is reasonable the spread of COVID-19 under 2020 Hubei lockdowns as in Figure I6
Next, give more details how to present components in Figure 16 with appropriate variables. First, we set the exposure A as if Hubei on lockdown, O for unlocked (before and on Jan 23rd) and I for locked (on and after Jan 24th). To select the candidate cities that contain most potential infected people, we rank cities outside Hubei by their received Wuhan (identified over 60% cases in China reported by the NHC) migration between Jan 1st. 2020 and Jan 22nd, 2020 (before the lockdown), and choose the top 30 cities (account for 69. 69.17% of total Wuhan migration) as mediators M, to control the noise. We use the daily migration scale index (MSI) of each city as the value of each mediator, defined in Baidu Qianxi to describe the migration magnitude. By noticing the following facts: 1) it took usually 2 days to diagnose the COVID- 19; 2) the estimated median incubation period is 5 days (Lauer et al. 2020); the outcome of interest Y is defined as the daily increase rate of confirmed cases out of Hubei with a one-week (7=2+5) delay. We delete extreme data points (two outbreaks at jails), and the final dataset yields total of 38 records.
21

---

Published as a conference paper at ICLR 2021
D.2 INITIAL DATA ANALYSIS
Figure 1a demonstrates the location of selected cities on the Chinese map with the color representing its cumulative confirmed cases by March Ist, 2020. It can be observed that the selected mediators either located around Hubei province or are big central cities such Beijing that has large population and migration scale. In addition, we provide the partial heat map of the correlation matrix between cities and the treatment A or the outcome Y as illustrated in Figure 17b Here, one could observe that cities' MSIs are highly positive correlated with the lockdown while are highly negative correlated with the daily increase rate of confirmed cases.
D.3 ADDITIONAL RESULTS
Table 6 lists all the numerical facts of the selected cities, including the population (million), the cumulative migration scale index (MSI) during the data period (Jan 12th to Feb 20th, 2020), the ratio of received Wuhan migration between Jan 1st to Jan 22nd, 2020 (before the lockdown), the cumulative confirmed cases by March stt, 2020, and cities' direct (DM) and indirect effect (1 M). Note the selected 30 cities in Table are sorted according to their cumulative MSI, and its order is used the order of the mediators in this paper. It can be seen from Table Table/othha the population, the cumulative MSI, the ratio of received Wuhan migration, and the confirmed cases are highly correlated among selected cities as expected. Note that we list these factors to assist the interpretation of the results, while none of them is used to estimate the DAG of interest. Therefore, similar values of factors don't 't necessarily imply similar causal effects.
Besides the general causal pattern of cities' estimated DMs and IMs stated in the main text, we provide more interpretation on the level of the individual city. Here, we compare three groups of cities that are of similar scale (population and migration) geographic position, to specify First, we compare the results between Beijing and Shanghai. where both cities have comparable population scale and are the center of politic economic in China (the correlation between the MSI of Beijing and the outcome yields the same value as of Shanghai as 0.50). It can be seen from Table that Beijing and Shanghai have similar positive indirect effects as 0.247 and 0.235, respectively, while Shanghai has a slightly higher direct effect on controlling the virus as -0.069, which is possibly due the smaller MSI of Shanghai. Second, we compare three cities located at Guangdong province (see Figure [ and 18b). including Shenzhen, Guangzhou. and Dongguan. All three cities have positive effects, among which Guangzhou yields the highest direct effect as 0.847, followed by Dongguan. The effect size of these cities agrees with their correlation coefficients with the outcome, where the correlation between Guangzhou and Y is 0.61, followed by Dongguan as 0.58 and Shenzhen as 0.50. The last comparison is among cities in southeastern China (see Figure T8a and T8b), including Suzhou, Hangzhou, and Wenzhou, all of which have a negative direct effect on the virus control and a positive indirect effect for the virus spread. Here, Wenzhou achieves the largest absolute value of the negative direct effect as -0.650, which conforms to its strict local shelter in home order after the Hubei lockdowns, where the correlation between Wenzhou and the outcome achieves the highest value as 0.83.
results.
We summarize all different sources of causal effects in Table the ANOCE table of 2020 Hubei lockdowns on reducing the COVID-19 spread in China. Note that due to cities' different levels of control measures on the coronavirus outside Hubei well as other possible confounders. there are some inconsistency between cities' DMs and their cumulative MSIs We leave the extension with confounders for further investigation. Further interpretation of the real data analysis requires domain experts.
One can refer Figure 18a and 18b for cities' DMs and IMS on the Chinese map. To check the reasonability of our results, we plot the estimated weighted matrix in Figure 4a where the first node (indexed 0) represents the Hubei lockdowns, the last node (indexed 31) is the daily increase rate out of Hubei, and the middle 30 nodes (indexed 1-30) correspond to 30 selected cities in Figure [b) From Figure 44a we can observe: 1). The color of the first column is almost all blue, indicating that locking Hubei down can reduce the migration of selected cities; 2). An approximate red upper triangular among first 20 nodes implies a migration trend with positive effects among central cities with large MSI, i.e. relatively smaller cities tend to have positive effects other cities with relatively larger MSI; 3). An approximate red lower triangular among last 10 nodes indicates a weaker migration trend with opposite direction among non-central cities with small MSI; 4). There is
22

---

Published as a conference paper at ICLR 2021
Table Analysis causal effects of 2020 Hubei lockdowns on reducing the COVID-19 spread in China regulated by Chinese major cities outside Hubei.
Source
d.f.
Causal effects
Direct effect from Hubei lockdowns
-0.078
Inirect effect via Chinese major cities outside Hubei
30
1E
0.419
Ratio Wuhan
Confirmed Cases
Cumulative MSI
Index
City name
Population
DM
TM
Shenzhen
10.36
3.02%
418
275.42
0.021
0.002
Guangzhou
12.70
3.22%
346
270.38
0.847
0.176
Beijing
23.00
5.50%
414
260.02
-0.003
0.247
Chengdu
14.05
2.57%
143
233.62
-0.201
-0.004
5
Shanghai
23.02
4.09%
338
232.80
-0.069
0.235
6
Dongguan
8.22
1.14%
99
211.99
0.257
0.003
Suzhou
10.47
1.42%
87
171.46
-0.481
0.553
Xian
10.00
1.87%
120
136.56
0.015
-0.089
Hangzhou
8.70
1.77%
169
135.73
-0.335
0.144
10
Zhengzhou
10.14
3.04%
157
132.58
-0.019
0.309
Chongging
30.17
4.97%
576
128.98
0,313
0.050
12
Changsha
7.04
5.21%
242
126.78
0.399
0.222
13
Nanjing
8.00
1.78%
93
103.87
-0.210
-0.100
14
Kunming
6,43
1.24%
53
95.35
0.186
-0.075
15
Tianjin
12.94
1.02%
136
90.50
-0.092
-0.372
16
Hefei
5.70
2.11%
174
87.28
-0.182
0.249
17
Nanning
6.66
1.01%
55
76.63
0.1
0.152
Wenzhou
9.12
1.02%
504
60.38
-0.650
0.200
19
Nanchang
5.04
2.15%
230
46.95
-0.109
0.095
20
Zhoukou
8.95
1.50%
76
42.01
0.045
-0.043
21
Fuyang
7.60
1.25%
155
39.26
0.031
-0.023
22
Shangqiu
7.36
1.12%
91
35.48
0.030
-0.057
23
Yueyang
5.48
2.31%
156
29.56
0.006
-0.050
Zhumadian
7.23
2.34%
139
28.86
-0.070
-0.075
Changde
5.72
1.05%
82
28.49
0.002
-0.003
Nanyang
10.26
1.91%
156
27.71
-0.039
-0.106
27
Yichun
5,42
0.90%
106
24.91
-0.051
-0.085
28
Xinyang
6.11
5.00%
274
24.78
-0.080
-0.027
29
Anging
5,31
1.57%
83
23.61
-0.022
-0.012
30
Jiujiang
4.73
2.06%
118
22.29
-0.072
-0.047
Selected cities
295.95
69.17%
S790
3204.24
30
-0.419
Total
1 p = 31 TE
0.497
23

---

Published as a conference paper at ICLR 2021
Sprees, Contel
(a) The estimated direct effect of cities ( (DMs).
(b) The estimated indirect effect of cities (LMs).
Mai gernad Path
28 10
The estimated DAG for the coronavirus data.
(d) The spreading network among citics that received most of the Wuhan migration during the data period (2.57% ~ 5.50%, presented by the color of nodes).
Figure 18: Additional analysis the causal effects 2020 Hubei lockdowns on reducing the
COVID-19 spread in China regulated by Chinese major cities outside Hubei.
also an almost all blue rectangle in the right top of the estimated matrix, showing non-central cities tend to have negative effects on central cities. Overall speaking, all the above finding accords with the migration trend during the Spring Festival period and intensive mutual communications among central cities in China, though there also some noisy causal directions opposite with the main trend in each area due to possible confounders and identifiability issue in the linear Gaussian model.
In addition, we provide the estimated DAG for the coronavirus data in Figure ISclas as the complete spreading network among major cities outside Hubei, where the exposure A, mediators M (cities' corresponding index can be found in Table G). and the outcome Y are colored red, blue and green, respectively. It can be observed that the in-degree is larger than the out-degree for nodes with small index, while an opposite rule is applied for nodes with large index. This finding is consistent with migration trend identified in our main text. Lastly, we give the spreading network among cities that received most Wuhan migration during the data period, including Beijing. Shanghai, Guangzhou, Shenzhen, Chengdu, Chongqing, Zhengzhou, Changsha, and Xinyang, plus Wuhan, in Figure [&d] to illustrate the partial interaction trend among cities. Each node refers to city with the color of the node presenting the percentage of received Wuhan migration, ranging from 2.57% 5.50%.
24

---

Published as a conference paper at ICLR 2021
M
M,
0.5
M,
0.5
0.4
M
0.7
Figure 19: weighted causal DAG g
E TOY EXAMPLE
Example E.1 Here, we use toy example of weighted DAG glos under the LSEM given in Figure 19to better demonstrate our definitions.
In g p mediators are included, where M1 0.2A + E,, M2 0.2M, + 0.5M, + e2. M3 & 0.4A +5, E3, M4 + 0.5M3 + Y = 0.5A + 1.1M2 + 0.7M3. There are 4 directed path from
to Y:
1) "A + Y' with length 1;
2) 'A + M3 + Y' with length 2:
"A > M1 M2 + Y with length 3;
4) + M3 M4 - M2 Y' with length
Since there exists 3 directed path x² E AT (ghou)} such that the length of x² is larger than 2, we have the mediators in Gloy are interacted.
From the weighted DAG Gloy, we have the direct effect of A Y is 2oy 0.5, a = (0.5,0,0.4,)) 3o0y = 00 1.1,0.7,00 OJT, and
0
0
0
toy
0.2
0.5
BM
0
0
0
0
0.5
Note one may recover Gloy from the weight matrix Bloy as long as the order of the vertices in Bloy given. Then, we have
0
(Lp Bi y
0.2
0.25
0.5
0
0.5
Thus, the indirect effect of A on Y
Bloy T
(Tp
Log
a toy
M
0.2
1
0.25
0.5
10, 0.7,0
[0.5, 0, 0]
(E.1)
0
0
0.5
= [0, 1.1.0...7,,000001. 0.2, 0.2
=0 0 + 0.22 + 0.28 + 0.5.
From the Equation E. and Theorem 3.1 we have the DM of M2 and M3 is 0.22 and 0.28, respec- rively, while other DM are O. Note that there is no explicit expression of the natural indirect effect through the directed path from A Y (IM) due to the complex interaction among mediators, while we provide its theoretical form in Equation [.66 based on the path method in Wright (192]) and Nandy et al. (017g) Specifically, for each directed path from A to Y, have:
25

---

Published as a conference paper at ICLR 2021
1) the direct effect through 'A (- 0.5;
2) the effect of path -A M3 + 0.7
3) the effect of path 'A -A- .M + x. 0.2
4) the effect of path 'A + M3 -
5) the indirect effect of A on Y - 0 x 0.5 x 1.1 0.11;
0.5;
6) and the total effect of A on Y (TE) is
7) the indirect effect for M1 (1 corresponds the effect of path 'A Mi + M2 Y', thus is 0.11;
8) the indirect effect for M2 (12)) zero since there is no path first goes through M2 followedd by other mediators;
9) the indirect effect for M3 (lMa) corresponds to the effect of path + M3 - M4 - M2 - Y
thus is 0.11;
10) the indirect effect for M4 (1M4) corresponds to the effect of path 'A - M3 - M1 + M2 - Y',
thus is 0.11.
We can calculate the last edge-specific effect directly from its definition. Since there is no M1 + Y and M4 -> Y gloy, we have LE LE = 0. By deleting M2 - in Gloy, the total effect reduced by 0.11 + 0.11 0.22, so LE 0.22; similarly, after deleting M3 + Y, the total effect decreases 0.28, thus LE = 0.28. One may notice that the last edge-specific effects are equal the DMs, and the additive of the last edge-specific effects is exact the last step of calculation Bloy (Lp - B toy -2 in Equation
F CONNECTION TO LITERATURE
We establish the connection between our proposed method the literature from three different angles. First, we show that the individual mediation effect defined in Chaaraborrtyy et al. (2018) can be decomposed into our defined DM and /M when LSEM assumption holds. Next, we give an equivalent definition of the DM through type of special edge (last edge) in the causal graph. Lastly, we prove that the proposed DM is consistent with the interventional effect via particular mediator defined in Vansteelandt & Daniel (2017) under the LSEM.
F.1
FROM INDIVIDUAL MEDIATION VIEWPOINT
Chakrabortty et al. 2018) defined the individual mediation effect under the LSEM as follows.
Definition F.1
Chakrabortty et al.
2018) Individual mediation effect for M₁:
17 EEM do a 1)} -)
(R.)
E(Y m. + 1)) (YY mm))]
In the following theorem, we show that the summation of n is strictly larger than the IE if the mediators are not parallel. The proof is given Section G.2
Theorem F.1 If there exists least one directed path n' E (()) such that the length of larger than 2, and the element in B is nonnegative, then
Σn>I EN>IE. IE.
(F2)
Remark F.1 From the above theorem, it is clear that the mediator effect defined in Chakrabortty lt al. (2018) is not appropriate for interpreting the decomposition of the indirect effect, when there exists interaction among mediators (a common situation described in the introduction). Here, keep the condition that the element in B is nonnegative, as the multiple count mediation effects in Chakrabortty al. (2018) may cancel out in some cases and their summation would equal the I by chance.
Inspired by the proof of Theorem F.I the mediator effect can be decomposed into two parts, the natural direct and indirect effect for i-th mediator, shown in the following corollary.
26

---

Published as a conference paper at ICLR 2021
Corollary F.1 Under assumptions (AI-A3) and Model we have
Tバ DM, + IM,
(F.3
Remark F.2 Corollary F.1 logether with the definition of th in Chakrabority et al. (2018) provides feasible way to numerically calculate the natural indirect effect IM, Specifically, by deleting the mediator M, in the causal graph, the reduced treatment effect corresponds no then IM, Tl DM,, where the explicit expression of the DM, is provided in Theorem 3.1 See more implementation details in Section
F.2 FROM GRAPHICAL PERSPECTIVE
Next, we give the definition of the edge-specific effect following Avin et al. (2005). Suppose directed edge of interest as X; + X, in weighted DAG G. Define a new weighted DAG 9 by deleting the directed edge X; + Xj in G, i.e. Fh = (X; Xj).
Definition F.2 Avin al. 2005) Edge-specific effect:
ET(X;, X,) TEg "
(R.4)
where TÈg means the total effect in graph G.
We next give equivalent definition of our proposed DM from graphical perspective. Let the edge G that starts with i-th mediator and ends with node Y, i.e. M₁ + Y, as the i-th last edge. Denote the graph G deleting the i-th last edge (M; + Y) G. We define the ith last edge-specific effect as
Definition F.3 Last adge-specifie effect for M,:
LE = { 0,
TEg T TEg
if there exists edge Mj - Y in G;
(F.5)
otherwise.
By Theorem 3.1 we have (Ip - -B) -1 ax is the causal effect of A on M. Let (T (J)-BB) with its i-th element G = {(Ip B B Next, we show that the i-th last edge-specific effect can be presented as BiGi under the LSEM in the following theorem, where 3; is the i-th element of vector (3 and corresponds to the weight the edge M; + Y. The proof can be found in Section
G
Theorem F.2 Under assumptions (AI-A3) and Model we have
LE, BiG..
(F.6)
Based Theorem 3.1 the natural direct effect of M; Y can be expressed as DM; BiG.. Thus, with the result of Theorem F2 it is easy to show the following corollary.
Corollary F.2 Under assumptions (Al-A3) and Model 7 the natural direct effect of Mi is equal to the i-th last
effect:
LF,, DM, = B.G
(F.7)
Remark F.3 Here, both definitions describe the direct impact of one mediator M, the outcome. The natural direct effect of particular mediator M; can be understood as the influence when removing the direct edge between Mt and Y. Thus, we have the equivalence between two definitions.
Then, we can decompose the total natural indirect effect into p last edge-specific effects or p DMS as the following additive form, based on Theorem 3.2 and Corollary F2
Corollary F3 Under assumptions (Al-A3) and Model
we have
IE IE= Σ
BiGe =
Σ DM, Σ LEg-
(F8)
1=1
1=1
i=1
27

---

Published as a conference paper at ICLR 2021
In fact, based on the uniqueness of each last edge, the natural indirect effect can be decomposed into p last edge-specific effect regardless the LSEM setting through the graphical perspective. We give the following intuitive conclusion. The proof can be found in Section G.4
Theorem F.3 The IE can be decomposed through LEs
1E LE
(F99
1=1
Remark F.4 One can view the last edge-specific effect as the generalized definition of the natural direct effect for mediator without the LSEM assumption.
F.3 FROM INTERVENTIONAL EFFECT LEVEL
Finally, we show the consistency of our defined DM to the interventional effect via a particular mediator defined in Vanstcelandt & Daniel (2017) under the LSEM.
Definition F.4 via M, is
<Vansteelandt & Daniel 2017 Under assumptions (AI-A3), the interventional effect
Σ
Σ
E(Y)A= M, m, o;) =o
M.
mp MpE
(F.10)
P(M, - mmaf=mm,A=aa))) P(M; =m:4A = =))
where M, the support of M, Og = [m] m. 11-1, M441 mpp. P(M = m A = a) is the
probability of M m when setting A
Theorem F,4 Under assumptions (AI-A3) and Model 7 we have
DM = E..
Remark F55 The proof can be found in Section G.5 Based on Definition 3.2 and Equation F.10 both the proposed DM and the effect defined in Vansteelandt & Daniel (2077) contain the information of
the causal effect of A on the mediator M. iee. P(M, m a + 1) P(M; m4; A = a).
G TECHNICAL PROOFS
G.1 PROOF OF THEOREM 3.1
Proof G.I In this proof, will give the explicit expressions of causal effects defined under the LSEM. First, Equation 3 is equivalent to
= CA: M aA + BLM +EM;
(G.1)
Y = 7A+B M +
Based M aA B eM, by moving B M to the left-hand side, we have
(4, 4p-BLM aA + Eas.
Suppose the mediators are sorted in the topological order (a series of elementary transformation of the matrix), then the matrix BY, is strictly upper triangular with the diagonal element 0. Thus, have 1, B is invertible, then I, B under its original order should be also invertible (any invertible matrix after elementary transformation is still invertible).
Therefore, we can rewrite M as a purely function of plus the error term as follows.
M (Tp BT- la + (Tp BT) teM-
(G.2)
28

---

Published as a conference paper at ICLR 2021
Then replace mediators in Equation
G.1 with Equation G.2 and obtain
A EA:
M (Ip B -1 -
(G.3)
Y = = + 3 M +
= + ( (Tp - - + (Tp - BT) + cy)-
Next, show how to get the explicit expressions of E a}} under the LSEM. Following the results in Rosenbaum & Rubin 1983), under the assumption (A2), we have P[M]do(( = a)} P(M a), and thus,
F[M)do(A = a)} E[M|A= a)-
Similarly, we can get E = a}} = E = a) under the assumption (AI), and FAY(do( A a, M =m=m]] = FEY a, M m) under the assumption (A3).
Based on above results and Equation
G.3 we have
E{Y = d)} = [[Y|A=A))
=E(7++ ey
=ya+ + 3 E(MIA= = a}
(G.4)
=7a+ 3 E( -B aAA (Ip B cul = a}
+ 3T (Ip -1)) B
where the first '=' is held under the assumption (Al), the second and forth '=' given by Equation G.3 that Y = 7A + _ M + cy and M (Ip BBs)) aA + (Ip -B -
Following the same calculation procedure of EdY|do( A a)}, give the natural direct effect under assumptions (A 1A3) and Model as
DE E(Y do = a
m(())) E(Y a)} = (a(a+ +)+ T (a) tra + ++33* mm)
7
where the first '=' is given by the definition of the DE. Similarly, the natural indirect effect is
IE = [do A = a, M = mlatt))) [ = a)}
(ya B m(at1)] {7a + 3T m())
BT (Tp B - 1 B (Ip B -L
= -3 (Ip -
Thus, the total effect of A on Y is
TE F(Y do
1)} = IE B
Finally, we give the expression for the natural direct effect of M₁ on Y under the LSEM. Based on the assumption (A2) and Equation G.2 have
a -))
=E(M = = - a}
(G.5)
=((I, (a) ((Ip
-1
=((Ip -
where ((I, BTy -ta)i is the i-th element of the vector (T, -1
29

---

Published as a conference paper at ICLR 2021
Then, based on Y 3 M + cy and the assumption (A3), we have,
PyY]do(A = a, M,- m; (a)
+ 1, = = E[Y]do(A=
=E a, M₁ m
(0)
+ +1 0.01 E{Y|A= a}
(a)
1 (a)
1
=ya + 3
m (a)
+1
7/0
m (a)
Bi.
(a)
(a)
mip
map
where 1, a p X 1 vector with the i-th element as while others qual to O. and Be is the i-th element of the vector B.
Thus, we have
DM, = (M
= a]}}]
[[Y|do a, M, = + 1,0, = =0; - - =))]
= [1 - x
-3, - B
G.2 PROOF OF THEOREM
Proof G.2 1. If there is no directed path nt E (G)) such that the length of n² is larger than 2, iee. the length of n* E (9)) is either or 2. Here, the path with length 1 corresponds to
A + Y, and paths with length are A M, - Y with possibly i == 1, .p. Thus, there is interaction among mediators.
By the definition of the LSEM, we have BM Opxp. where 0pxp is a p zero matrix. Following the path method (the causal effect of X, on X, along a directed path from X, + X, in 9 can be calculated by multiplying all edge weights along the path) illustrated in Wright (922]) and Nandy et al. (2017). we could obtain Σ = -, 3,a; = IE. (See a toy example provided in section to illustrate how to use the path method to manually compute the causal effects.)
2. If there exists at least one directed path x³ E (9) such that the length of n²" is larger than 2, and the element in B is nonnegative, we have BM # Opxy Without loss of generality, suppose there exists M; E M with a set of directed path that starts with A, contains M, then goes through other mediators, and ends with Y. denoted each path in such set as Ttj = [AA+ - - M,.-.
Y} for 1.... , Tii where n, is the size of such path set for M, and the weights of edges in Tij is positive. Note the set {T..a) excludes the paths end with M, +
Let ex, denote the causal effect of A on Y through directed path Ta.j- Based on the path method in Wright 1921) and Nandy el al. (2077) with the definition of MM, we have its theoretical form as
TM:
Caaa'
(G.6)
By Equation G.5 and the definition of Tee we have its first multiplier
E EdMildo A = a + 1)} E[M])do( A a)} ((Tp - -B -aaii
which is also the first multiplier in both DM, and 1M,
And the second multiplier of n can be expressed as
E{Y do(M, = m4 + 1)} E{Y (o = m.)}
-EAY|do[MM mm.)) (a) + 1)} ERY [do(M, m (2))}}
=E{Y (do A a, M, = =m,") (a) 1)} - E(Y = =a, M, m( )
=E'YY |do[A = a, M; = = + +1} - E Ido(A a}}
30

---

Published as a conference paper at ICLR 2021
where mi (o) is the value of M, when setting do(A a). Here, the first '=' is valid since my can be
arbitrary number, and the second and third "=' are based on the equivalent interventions.
Based the technique of plus and minus the same term, decompose the second multiplier of di into two parts as follows
[{Y|do(A a, M4 m (0) + - - E(Y
E(Y (EYYdo((A-)) - E[Y]do(A a, M, m! 1, -.,02,0,01) 0.001)
the vecond multiplier of DM;
+ E[Y|do(A- = a, M, m, (0)
+ 1)} E( M, m, (a) +1, = of 0/
the second multiplier of Mr
(G77)
where Ω2 = M M; is the sets of mediators except M,. and o0 is the value of Ω when setting do(A = a). Here, the first term in the above equation corresponds to the second multiplier of DM, while the second term the second multiplier of 1M4
Thus, the summation of T): is
Σ Σ ((M,)d((A=aa1)) -(MA))))))M
[EtY do(M; my 1)} EYY m.))
E(DM +MM)+IMM)- DM. Σ IM IE+
where the first '=' is from Definition F.1 the second '=' is given by Equation G.7 and Definition 3.2 and 3.3 and the last "=" comes from Theorem 3.2 and the theoretical form of IM in Equation G.6
Here, have ex, 0 since the weights of edges in viij is positive based the path method in Wright 1921) and Nandy et al. 2017)). Then, Σ, E Cns, is also strictly larger than 0. Therefore, we have
Σ Σm>L IE.
G.3 PROOF OF THEBRRE E.2
Proof G.3 1. If there doesn't "resn't exist edge M₁ - Y in G, then by definition we have B. 0. Thus, LE, B.G. =.
2. If there exists edge M. +Y in G. Suppose there is a directed path set with size m, associated to the edge M; Y, where each directed path #.,j starts with node A and ends with M, Y, denoted as म.. = {A -
M, > Y} for 1,... m.,
(A.M) the causal effect of Let eas., denote the causal effect of A on Y through directed path Re.g, G
A on Mr through directed path Ja.j, and (AL,,Y) the causal effect of M, Y through directed edge M; + Y. Following the path method in Wright [921] and Nandy et al. (2017), we have (A.Mi) (M.Y).
ea,, Cx,,
Then the i-th last edge-specific effect is equal to the summation of the effect through each path R..j. ie.,
LF,-
CR..
(A.M)
(AMM)((MM.))
(M,.)
(A.M.)
1=1
1=1
Here, by the similar argument based on the path method, have e (M.Y) Br and Σ)= (A.M. the total causal effect of A Mr
Cris
31

---

Published as a conference paper at ICLR 2021
Recall that S {(Ip ) a is the causal effect of A M.. Therefore, the i-th LE is the product of the causal effect of A on M, and the causal effect of M, on Y, i.e,
LE;
BiGi-
G.4 PROOF OF THEOREM
Proof G.4 Given a general DAG 9 with nodes MA, let the union of all directed paths that contain the i-th last edge T; { A 7 + M, Y},i 1,.-.p. Here, have
Ti = { hssssm established in Section G.3 It is clear that the union set of Tq in 9 is equal to the of all directed paths start with A and end with node Y (except A Y) in 9 as
- [mAY()) (A-
Also, based the uniqueness of each last edge, Tú is pairwise disjoint, i.e.
Ts T -,
Vi + j.
Since the EE is defined as the total causal effect of on Y that goes through mediators, we have the TE equal to the causal effect that goes through the set AY (g) {A Y}, i.e. the TE equal to the causal effect that goes through set Up Tr. Based on the mutual disjoint property of To we have the IE is exactly the summation of the causal effect through Ty. Lastly, from the definition of LEd have
IE
LE:.
G.5 PROOF OF THEOREM
Proof G.5 The proof of the consistency of our defined DM the interventional effect &, can be completed based on Equation under assumptions (Al-A3) and Model
Recall the definition in Equation 1.10
we have
& = Σ
Σ MR,
E(Y a, = m,, = oi/A=
m EM,
{PM [ PMM m. m, at)
Σ
Σ
= 0. E(Y|AA a. 0, P00
a = m. a
MjEM,
my EM,
E(Y|A= a. m,, Ω = P(1, 0. - a) PPM mm/A= a)
Given A = a, the value of M, is (a) and Ωt takes (a) while when setting A = a + 1, the value (i)
of M, m Therefore, we have P(M, miIA a) = if m, mi otherwise is O, and 1 if o; g:a) otherwise is 0.
P(Ω A A
Under assumptions (AI-A3), we have
(a+1) ΩR: = 0000)) B(YIA a, M, m, (0)
Ω (a))
& = A M; m
Then, based on the LSEM that Y = 7A + 3T M + ev, we can further obtain that
m (a)
(0)
& = T
m
(1+1)
ya
T
m (
B Bm
(a+1)
m
(4)
mp
mp
32

---

Published as a conference paper at ICLR 2021
From Equation G.2 we have
= 3: [10, BM) a):(aa
B + 1) {(Tp - aa)a)aa = B.(lT, - BBy -1 'a)..
Thus, under assumptions (Al-A3) and Model 7 we have
DM = Et-
33


Reviews:
[{"review_id": "mvVpCpwvqIJ", "reviewer": "AnonReviewer3", "summary": "", "questions": "", "limitations": "", "rating": 5, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "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": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "##### 1. Summary\nThe paper proposes a framework for the analysis of causal inference. Its main contribution is to decompose the indirect effect by teasing out the causal contribution of a set of mediators. In a series of experiments with simulated data the authors show that the proposed method, ANOCE, outperforms other comparison partners. The manuscript also contains an analysis of real-world data that describes the causal effects of the lockdown of cities in the Hubei province (China) to reduce the spread of COVID-19. An extensive supplementary file is also part of the submission and it includes additional experiments and technical proofs.\n \n##### 2. Rationale for the score\nThe paper is well structured and very well written. The main document provides all explicit definitions of the causal effects under the linear structural equation model (LSEM). The appendix contains additional technical proofs. I find the notation to be consistent and easy to follow. Despite these pros, I find the results to be inconclusive and I suggest additional experiments that will make a more persuasive case for the method.\n\n##### 3. Positive aspects\n- the problem is well formulated and the manuscript is clear and well written\n- the authors show how factors that affect the performance of other methods, such as the number of mediators or sparisity, have no negative impact on the performance of ANOCE.\n- Tables 2 and 3 in the appendix provide a convincing argument on how ANOCE can disentangle the natural direct effect of the mediators as the sample size increases.\n- I commend the authors for using ANOCE to analyze real-world data. The analysis of COVID-19 spread in China is interesting, although some aspects of the results do not seem convincing (see section 4 below for more details).\n\n##### 4. Negative aspects\n- For the comparison partners, the text indicates that default hyper-parameters were used. It is not clear for which methods the parameters were fine tuned, if this was done at all. As shown in Table 4, NOTEARS and DAG-GNN perform worse than PC in some settings of scenario 3*, with respect to structural Hamming distance. This begs the following questions: i) for the VAE in DAG-GNN, what are the default parameters? ii) for NOTEARS, Zheng et al. performed an analysis of the sensitivity of the threshold. They report a choice of w = 0.3 for their experiments, making clear that the choice of threshold becomes problematic as the sample size shrinks. Although the authors of the current manuscript also chose 0.3 as threshold, the question is how appropriate is this threshold for the simulated dataset they use. Note: In contrast to the 500 samples generated in the current manuscript, 1,000 simulated samples were created in Zheng et al. 2018. A fair comparison would have been to compare against the same simulated datasets described in Zheng et al. 2018 (see details below).\n- In my opinion, some additional experiments are needed to create a convincing case for ANOCE. From the description in the Appendix (Section C.1) it seems that only one random graph was generated. Of course, it is understandable if the plots shown in Figures 2 and 3 only depict one true graph. But in order to compare against other methods, the experiments must consist of hundreds of different random graphs. Two considerations: i) the expected number of edges for each node is set to 2. Why not try other settings where the expectation is 1, or 4? ii) the random graphs are created using the Erdös-Rényi model. To make a stronger case for the method, the experiments should contain random graphs generated as scale-free networks, as detailed in [R1]. Combining i) and ii), these random networks are usually named ER1, ER2, ER4, SF1, SF2, and SF4.\n- With respect to the analysis of COVID-19 spread, there are some aspects that are not clear: i) how to interpret the values in Table 5. For example, Beijing and Changsha have a similar value in the “Ratio of Wuhan” column. The IM of both cities is also similar but the DM of Beijing is close to null. It is not clear if the DM and IM values can be interpreted properly. Therefore, this makes me question the added value of this table. ii) Figure 10.c is too small to see and no meaningful insights can be extracted from it. One suggestion: if each node is linked to a real city, plot the figure on a map of China. The nodes A and Y can be outside the map (similar to Figure 10.d). iii) what do the colors mean in Figure 10.d?\n\n##### 5. Questions to be addressed during rebuttal period\nUnless I’m mistaken, the current results seem to be for only one random ER2 network (one with 50 nodes and another with 500 nodes, to be precise). Please, clarify the generation process for random network(s).\n\nAdditionally, it will be useful to learn how ANOCE performs in other types of networks (ER1, ER4, SF1, SF2 and SF4). \n\nWhat hyper-parameters were used for the comparison partners?\n\n##### 6. Additional references\n[R1] Barabási, Albert-László, and Réka Albert. \"Emergence of scaling in random networks.\" Science 286, no. 5439 (1999): 509-512.\n", "extra_scores": {}, "extra_text": {}}, {"review_id": "azayf86BPNB", "reviewer": "AnonReviewer4", "summary": "", "questions": "", "limitations": "", "rating": 6, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "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": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "This paper combines two contributions to causal inference: a definition of and distinction between direct and indirect mediator effects to take account of cases of more than one mediator on a mediating path and how to identify them in linear structural equation models, and a modification of the DAG-GNN method for learning causal structure to incorporate certain background knowledge. In addition, the paper applies the concepts and the method to some interesting COVID19 data. \n\nI am not sure how significant the first two contributions are; they strike me as primarily incremental. However, I am inclined to recommend acceptance in my overall judgment.\n\nSome specific comments/questions:\n\n1. Remark 3.2 is a bit misleading. Although the background setup is supposed to be a linear model, it's better to make the qualification clear *within the remark* that the decomposition refers to decomposition in linear models rather than decomposition in general.\n\n2. The simulation uses primarily linear Gaussian models, whose structure is not always identifiable even given the background knowledge assumed in this paper. That is, sometimes the direction of an edge between two mediators may not be identifiable from observational data. It's unclear to me how this issue is handled in the experiments. An arbitrary member in the equivalence class is used to estimate the various effects?\n\n3. Similarly, with the COVID19 data, the determination of some of the directions in the discovered network may be quite arbitrary. Did the authors notice any causal direction that does not make much sense or accord with the known migration trend?\n\n4. The method is to said to also apply when the exposure is categorical. I would appreciate some clarification on how this variation affects the definitions and the identification formulas. And how about categorical mediators?  \n\n5. On p. 3 there is the following remark on methods based on conditional independence tests: \"However, these methods have no guarantee in the finite sample regime because of the hardness of testing conditional independence when the variables are continuous (Shah & Peters, 2018).\" This remark is a little misleading. The lack of finite sample guarantee in the absence of substantial assumptions is not just a problem with such methods. It is everyone's problem, as argued in, for example, Robins et al. (2003) \"Uniform consistency in causal inference\".  ", "extra_scores": {}, "extra_text": {}}, {"review_id": "aRIOxri-Gxz", "reviewer": "AnonReviewer1", "summary": "", "questions": "", "limitations": "", "rating": 6, "confidence": 2, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "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": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "### Summary\nThe authors propose a framework to analyse causal effects, considering interactions between mediators. The causal structure learning includes a temporal causal relationships of variables. Also, the authors validated the framework on real-world data set, in particular the Chinese coronavirus outbreak near Hubei.\n\n### Reasons for score\nThe paper proposes an interesting and useful framework. The results show it can outperform the baselines (figure 3). The problem is clear and well motivated. Also, the technical details (section 4) are well explained, allowing for the reproducibility of the paper.\n\nThe proposed approach seems to be an extension of the DAG-GNN. Also, the latter is shown to have similar results in the analysis (fig 3). One main difference proposed by the authors is that ANOCE-CVAE can handle non-linear cases, contrarily to DAG-GNN. In [1] the authors mention that they method is actually able to handle parametric models beyond linear, so I may have missed something. Probably a way to improve the paper would be adding more details on the contribution with respect to the state of the art. \n\n### Pros\n- Analysis on synthetic and real data sets shows proposed model to outperform baselines.\n- Detailed explanation on the technical part.\n\n### Cons\n- Not super clear what are the improvements with respect to the art.\n\n### References\n[1] Yu, Yue, et al. \"Dag-gnn: Dag structure learning with graph neural networks.\" arXiv preprint arXiv:1904.10098 (2019).\n", "extra_scores": {}, "extra_text": {}}, {"review_id": "97caEOKl7N4", "reviewer": "AnonReviewer2", "summary": "", "questions": "", "limitations": "", "rating": 8, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "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": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "This paper presents a novel approach to structure learning for cause-effect models, focussing on the influence of mediators.\n\nThis approach is supported with an algorithm, which is implemented and tested on various datasets. The authors summarize important preliminaries and recap related approaches. The derivation is driven mathematically, supported by definitions and theorems. \n\nIt is remarkable that the authors build on Rubin's potential outcome framework and Pearl's graphical models and the do-formalism, where these \"schools of thought\" have been developing in somewhat seperated strains.\n\nI recommend to accept the paper.\n\nOverall, the paper is very well-structured, the accompanying appendix is detailed and provides extensive supplemental material.\nThe approach is driven theoretically and tested experimentally. \\\nThe choice of COVID-19 data as experimental data adds a good final touch,  taking a recent and prevalent real-world example into account. \n\nOne downside is that no code implementation is provided (e.g. as supplementary material). ", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Meta-review: This paper proposes an approach to learn the causal structure underlying a dataset with acyclicity and other structure constraints, and then used the inferred structure to compute partial causal effects. The authors show that, on simulated data, the proposed method outperforms others in the literature. The manuscript also contains an analysis of real-world data that describes the causal effects of the lockdown of cities in the Hubei province (China) to reduce the spread of COVID-19. 

Overall, the reviewers think that this is a well structured and written paper. From a novelty viewpoint, the main contribution consists in formalising the causal contribution of mediators, as the method for computing the causal structure is based on a small modification to previous literature.

The main concern raised by the reviewers were on the experimental evaluation. Some of these concerns were addressed by the authors during rebuttal, whilst some on the number of nodes remained. We encourage the authors to consider these concerns in the final version of the manuscript.
Decision: Accept (Poster)