Analysis of Variance of Multiple Causal Networks
Authors: Zhongli Jiang
Dabao Zhang
Venue: NeurIPS 2023 poster (2023)
Keywords: causal inference
large graphs
multi-task learning
structural model
directed cyclic graph

Abstract: Constructing a directed cyclic graph (DCG) is challenged by both algorithmic difficulty and computational burden. Comparing multiple DCGs is even more difficult, compounded by the need to identify dynamic causalities across graphs. We propose to unify multiple DCGs with a single structural model and develop a limited-information-based method to simultaneously construct multiple networks and infer their disparities, which can be visualized by appropriate correspondence analysis. The algorithm provides DCGs with robust non-asymptotic theoretical properties. It is designed with two sequential stages, each of which involves parallel computation tasks that are scalable to the network complexity. Taking advantage of high-performance clusters, our method makes it possible to evaluate the statistical significance of DCGs using the bootstrap method. We demonstrated the effectiveness of our method by applying it to synthetic and real datasets.

Analysis of Variance of Multiple Causal Networks
Zhongli Jiang
Department of Statistics, Purdue University, West Lafayette, IN jiang548@purdue.dd
Dabao Zhang
Department of Epidemiology Biostatistics, University of California, Irvine, CA dabao @uci.edu
Abstract
Constructing directed cyclic graph (DCG) is challenged by both algorithmic difficulty and computational burden. Comparing multiple DCGs is even more difficult, compounded by the need to identify dynamic causalities across graphs. We propose to unify multiple DCGs with a single structural model and develop a method simultaneously construct multiple networks and infer their disparities, which can be visualized by appropriate correspondence analysis. The algorithm provides DCGs with robust non-asymptotic theoretical properties. It is designed with two sequential stages, each of which involves parallel computation tasks that are scalable to the network complexity. Taking advantage of high-pete biigh-prformance clusters, our method makes it possible to evaluate the statistical significance of DCGs using the bootstrap method. We demonstrated the effectiveness of our method by applying it to synthetic and real datasets.
limite-informmtion--ee
Introduction
Study of causal networks described by directed cyclic graphs (DCGs) becomes increasingly popular in variety of fields, such as biomedical and social sciences [10 3) Detecting dynamic structures across networks sheds light on mechanistic dynamics, with promising applications. For example, gene regulatory networks related to various types of cancer differ from that of healthy cells in complicated way and discovery of deviated regulations helps identify cancer-related biological pathways and hence design target-specific drugs [23] 18). Interpersonal networks in social media have become a favored way for people get informed breaking news and natural disasters nowadays. Understanding their evolution over time may provide promising information to businesses, educators, and governments 77. However, the complex structure and large number of nodes impose challenges on constructing even single network, not to mention modelling and comparing multiple networks.
Structural models are widely used deseribe causal networks and conduct causal inference [8] 15)
With available instrumental variables (IVs) from genomic variations, many methods have been developed to construct large causal networks of gene regulation by integrating transcriptomic and genotypic data. For instance, Liu et al. [14] proposed to assemble driver-respone identified from local structural models, formidable task due to enormous numbers of possible uuiver-responder pairs. Cai et al. [4] developed sparsity- aware maximum likelihood (SML) method by putting & penalty on causal effects. may reach local maximum, leading to the identifiability issue. Chen et al. [ proposed two-stage penalized least square (2SPLS) approach which allows parallel computing and shows superior performance.
relationships
To infer sparse differences between two large networks, Ren and Zhang [6] recently developed parallel algorithm based 2SPLS and achieved better performance than separately constructing the
37th Conference on Neural Information Processing Systems (NcurIPS 2023).

---

two networks using 2SPLS. Zhou and Cai [25] proposed a fused sparse SEM (FSSEM) method to learn the differences through maximum likelihood inference of joint networks which is computationally infeasible with large networks. Li et al. 133 employed a Bayesian scheme in their method BFDSEM, unrealistically assuming equal sample sizes from both cohorts.
To the best of our knowledge, there is no algorithm available for investigating causal effects varying across multiple cohorts while identifying stable ones. In this paper, we develop novel algorithm, designed for parallel computation, to construct a unified structural model for multiple causal networks. While each causal network may be depicted by a DCG, we deliver an analysis of variance (ANOVA) of these networks to identify causalities that are different across networks as well as important drivers and responders. Our algorithm is scalable in several aspects. Firstly, it is scalable to data size including sample size and number of variables. It would still be efficient and powerful under high dimension scheme. Secondly, is scalable to the computational environment with its parallel computation. The fast computation allows to attack a daunting task in studying multiple networks, i.e., calculating p-values via the bootstrap method and thus controlling the false discovery rate. Thirdly, the algorithm is scalable to model complexity as it is able to infer causal networks beyond directed acyclic graphs (DAGs), such as DCGs which are demanded to depict gene regulatory networks [14].
The rest of the paper is organized as follows. We first state the model and discuss its identifiability guaranteed by available instrumental variables in Section 2 In Section 3 we present our proposed algorithm ANalysis Of VAriance of directed Networks, termed as NetANOVA, and introduce mea- sures to quantify an endogenous variable's contribution as either drivers or responders. The theoretical justification is shown in Section 4 with detailed proofs in Supplemental Material. We demonstrate the feasibility and promise of our algorithm with a large-scale simulation study shown in Section [ and real data analysis to compare gene regulatory networks in healthy lung tissues and lung tissues We conclude the paper with a discussion in Section 7
with two types of cancer in Section 6
2 A Unified Structural Model of Multiple Causal Networks
2.1 Model specification
We consider K pertinent causal networks, each for a cohort. For each, say k-th, cohort we have n (k)
observations in (x()), (()), where (() is an (*()) p matrix including values of p endogenous variables, and X is an (2()) x q matrix including values of q exogenous variables. Each variable is assumed to have mean zero. For the k-th causal network, we consider each, say i-th, endogenous variable is causally affected by other variables as follows,
(k)
(k) (k)
(R)
+ E
(I)
where Y²!k) is the i-th column of Y and (K) is the submatrix of (k) excluding the i-th column.
The (p- 1) dimmennsill column vector (/
(k) includes the direct causal effects of other endogenous variables (drivers) on the i-th one (responder). The set I, includes the indices of exogenous variables
that serve as IVs of the i-th endogenous variable so X (k) is a ( x [Z,| matrix including values of corresponding IVs. Each IV has a direct causal effect on its corresponding endogenous variable and can be identified through domain knowledge. &1. includes causal effects from corresponding IVs.
(k) includes disturbance errors which are independently distributed with mean zero and standard <
deviation o7
Pooling together equation for each endogenous variable, we model the k-thh causal network as,
(k) (() (() (k)) (R)
which describes DCG as shown in Figure
The p x matrix (k) describes the causal relationships between each pair of endogenous variables and has zero diagonal elements to prohibit self-regulation. The q x matrix (R) encodes the causal effects of exogenous variables on endogenous variables. The (4) x p matrix (*) includes all disturbance errors. We assume that both ( and (k) have been appropriately centralized within the cohort, no intercepts are needed in the above model.

---

DCG
DCG II
DCG
Figure 1: DCGs of three causal networks. The nodes outside the shaded regions are exogenous variables with IVs in light blue and confounding variables (unobservable U and observed O) in white but disturbance errors are omitted. Causal relations in red differ across networks, and dashed ones cannot be revealed by available data due unobservable U unavailable IV for Ys..
2.2 A unified structural model
Many algorithms have been developed to construct structural model for single causal network [14 4 5]. We will show that, in the interest of comparing multiple causal networks, we can also unify them into single structural model and hence develop an appropriate algorithm to construct all networks at the same time to deliver an ANOVA of multiple networks.
Suppose have baseline, say K-th, network and are interested in others' deviation from the baseline. We first reparameterize the causal effects of all endogenous variables in (I with
B = 7/th) (K) for k = 1,2.
= -
1;
(3)
Therefore 3(k) includes baseline effects and, for l: * K, ((()) includes deviated effects of k-th network from the baseline. When it is of interest compare the networks with each other, we can similarly reparameterize to consider the deviated effects of each network from the average effects.
For any K sets of x matrices, say A, with 1,2, ... K, we define matrix-valued function, T(A .... AR) = (AI:K, (diag( (A,,
- (,(K))))"",
where AI:k is constructed by stocking A1, ... Ax row-wisely, and Om((-1) is m(K - 1)- dimensional column vector with all elements zero. Further denote
Y, = (Y(II))) y(2)r
(K)TT
B (g( (K)
(1T) (2)
[K
Φ1. = (ØT (1)T
OT.
(2)T
(K)T (KYTT
ΦT.
)
Ta
6; (c (L)T
(2)T
Notice that (1) = kkk] + (h) for lo = 1,2, .....
from
we have a unified structural model
for K causal networks described by (2).
Y; = T (1) (2) ..... Y (K) | diag(2 XX
X x +
(4)
2.3 Instrumental Variables, model identifiability, and confounding variables
As shown in Figure [ the ability reveal causal effects relies on the available IVs. An IV is an exogenous variable that only affects directly the driver but not the responder, except indirectly through its effect the driver. This means that IVs can be used to address the model identifiability issue, which is complicated by the model's inherent endogeneity, and allow us to isolate the causal effects of endogenous variables. We specify the following assumption on available IVs following the rank condition (20]. a necessary and sufficient condition for model identification.
; b. 30 c (1,2, ... .p) but + 0: Vi C. 3I, C (1,2, ...) Assumption 1. a. (k) Il (k);
with I, *f 0 and d (1) * 0: d. Vi.j C with $J.I 0: Vj & C, 1i 0.
(k)
3

---

Assumption La specifies IVs observed in ( that are uncorrelated to the disturbance errors. Assumption 1.b specifies the endogenous variables in C that have available IVs and Assumptions I.c and 1.d state the restrictions on the available IVs, i.e., their availability and uniqueness, respectively. While Mi is the j-th component of 7(,), Assumption 1.e states that we are not able to identify any causal relation with j-th endogenous variable as driver if it has IV available, e.g., Ys in Figure
striking advantage of the IV method is its robustness to confounding but exogenous variables, which may be observable or unobservable, as shown in Figure [ As will be demonstrated by our later algorithm, available IVs help disassociate an endogenous variable from disturbance errors of other endogenous variables in model fitting. Therefore, we may explicitly incorporate observable con- founding variables to improve causal inference. However, when there are unobservable confounding variables, can still effectively construct multiple networks with available IVs.
3 Model Building and Interpretation
The single structural model in (4) makes it possible to employ limited-information approach to construct multiple causal networks even in the case of a large number of endogenous variables. Therefore we develop two-stage algorithm NetANOVA which allows parallel computing for fast computation and hence bootstrapping the data for significance assessment. Another daunting but necessary task is to understand the multiple networks and catch the important variables in the networks and important dynamic causalities across them. We therefore propose coefficents of determination and cause as well as correspondence analysis.
3.1 The algorithm NetANOVA
The disassociation stage:
We first predict each endogenous variable in C, solely based on all
available exogenous variables, to disassociate it from the disturbance errors of other endogenous variables, following aforementioned assumption 1. For this purpose, we rearrange the terms (2) to obtain the reduced model,
( = x (1) π (N) + (*)) where (k) (b() (I (()) and (**) c ( (I- (
(5)
As shown in the real data analysis in constructing gene regulatory networks, the number of exogenous variables, i.e., the dimension g for (*), can be much larger than the sample size (2()) rising issues of prediction consistency and computational time. To address them, we first apply the iterative sure independence screening (ISIS; (9) to screen for exogenous variables and then perform regression with &g penalty. As later theoretical analysis shows, allows for q S (x)( for some 0 0.
Specifically for each endogenous variable i E C and E (1,..., K). we apply ISIS with computa- tional complexity (n to screen for set of exogenous variables, indexed by set (* Using these d = (M (k) exogenous variables, we apply ridge regression to obtain the ridge estimator,
π Mg
= (X
MM!h M( +
(8)
m
(6)
and predict y(K)) with
(1 (k) Mn( (000
म (6)
where (a(a) is a tuning parameter that can be selected via generalized cross validation
The inference stage: We use the predicted endogenous variables in the previous stage to identify and estimate the causal effects, hence constructing the causal networks. We first calculate the projection matrix,
P (k)
X
(k) I,
X
(K)T T
T
T,
Multiplying the projection matrix P, = diag (P p(2) can eliminate the exogenous variables from the model and get,
p(K)
to both sides of model (4 we
P,Y. P.Y B + P.e.,
(8)

---

where -- T Y -i.... Y( Since P,T-, and 6, are correlated, we instead regress P,Y, against P,Y-, with r= T Y ( which is disassociated from E,. With possibly high-dimensional , we apply adaptive lasso 226 to obtain its estimator,
3 arg min 9
|PYY
P no -
n
where 24 is a tuning parameter, and w; 1B - for some > 0 where (Badl contains the absolute
values of B which is preliminary estimator. The original networks can be recovered by calculating
(k) (k) 3 + B for k 1, .,2 -
%
In summary of these two stages, we summarize the algorithm Algorithm
Algorithm Aalysis Of VAriance of directed Networks (NetANOVA)
Input: (Y (()), & E {1, 2, K}, with each variable centralized within cohort but scaled 1-1)
according to the baseline cohort; Predefined index set I; for each i E C; <- (n STAGE
for E Co
1. Reduce the dimension of (() by ISIS get X (k)
setl Set X mnlh
x if if (k))
2. Estimate Y(!)) by regressing ( against X (4) with ridge regression.
end for STAGE 2
for i 1.2.....PP
1. Calculate projection matrices P,.
2. Predict Y_ from ( C. 3. Estimate B, by regressing P,Y, against P,Y with adaptive lasso. end for
Output: B which contains the baseline and differential causal effects.
3.2 Coefficient of determination and coefficient of cause
For each endogenous variable i, we can calculate its coefficient determination to measure the proportion of its variation due to the effects of its drivers, i.e., Y (k) (R)
1 for cohort k;
2(k)
| (1()
On the other hand, we can also calculate the coefficient of cause for each endogenous variable i which summarizes proportions that contributes the variation of its responders, for cohort k,
(2(R)
||X
(8) (k) (k)
(k)
While we name the above coefficient of determination as causal RD with its value between zero and one, we will simply call the coefficient of cause 22 which is positive but may be greater than one.
3.3 Correspondence analysis of causal effects
We can conduct correspondence analysis of causal effects in (k) to reveal clusters of drivers, responders, or driver-responder pairs which outstand from the rest in each cohort k. We can also conduct correspondence analysis of deviated causal effects to reveal these clusters which deviate the most from baseline or the rest. However, each causal effect or its deviation may vary differently and we need to standardize them, based on the bootstrap results, before correspondence analysis. For example, to compare cohorts k and / for their causal effects, should instead obtain the following
5

---

compared to that of learning in a single-cohort analysis, where the estimators (a for each cohort are estimated separately using the algorithm without re-parameterizatioo. For each k. summary statistic,
- 7ij -2//0)))///00000000 /(B-)) 1/2
(1) 26j 11g Yis
(k.1)
(5 (k)
(k.b) Ois
where Voj and (k.b) Vaf
are estimates of corresponding effects from the observed data and b-th set of when both of its numerator and denominator are sufficiently small. We can define similar statistics to investigate effects' variation from their means.
(00)
bootstrap data, respectively; we also let =
singular value decomposition (SVD) of 2(k.) (zij (((..1)))
"pxp can obtain its left and right singular
vectors, say (U/4:0} and (v,(k.)), respectively. While the left singular vectors help identify respon- ders that differ between the two cohorts, the right singular vectors help identify drivers that differ
between the two cohorts. We can overlay the plot of U,(h) vs. U7(h) with the plot of Vi(a) vs. v-(k) to identify responder-drivee pairs for their causalities varying the most across networks.
4 Theoretical Analysis
In this section, we will establish non-asymptotic guarantees and show that our constructed DCGs have good theoretical properties. We characterize the properties with a prespecified sequence g(k) exp(-o exp(-o(n(()) with g(k) + 0 (8) + oc, i.e., each g(k) approaches zero slower than (x)(-n(A)). We denote 8min minısAs k ((h) and only present main results here, with details including general notations, assumptions, and proofs in Supplementary Material.
Theorem 4.1. Let Ek=
supp(B,), and In is function of n and 8min specified in 9 in Supplemental Material. Then under Assumption I and Assumptions B.1 B.4 in Supplementary Material, we have that, with probability at least 6min 5 with 5
>> k=1 8(h),
(Bounded Errors) 1|3, - 3.11 S |S:] 9n (d log(pK/0) V T 112maax v log(d/ómnin)) /n:
2. (Causality Consistency) (f.) (3g)
We will next show that, with proper choice of { we can bound the estimation errors system- wiie with ultra high probability. We can pick { such that dwin X e-n where E (0, min 0, 0)). Then we have fn S (2 (220))2, and the restriction on the number of true signals can be reduced to |S S 9/2 which is requirement that can be fulfilled especially in sparse model. We have d d=n' 1-8 s0 9n S re (1-20/22 As result, the bound can be dominated by the order of (1,1n(1-40)/2
Note that the error for the whole system of p nodes can be controlled by similar bound, replacing each occurrence by |S,| with max [S,]. with probability at least - p(8mmn + 8). Following the 15i50
former calculation, we learn that, when max [S|] grow slower than n. (4/-1)/2 the error will approach 1S*Sp
zero when n - co. To achieve the bound with high probability, only need põ diminish. That is, could have the error well controlled even if the dimension grows up to e. where ().
Given the auspicious estimation performance, we will further discuss the benefit of pooling cohorts
separately conduct the analysis on each cohort and obtain the following property for single task estimation.
Corollary 4.1.1. Under the same conditions as Theorem 1 8min 8/K,
4.1
we have that, with probability at least
1 | 7; 10 1 15 (d log(pK/8) V
| 2mma (on(dd/mmin)///m
The proof directly follows by treating each cohort baseline cohort in Theorem 4.1 Note that the denominator making a crucial difference in inferring the bound. The algorithm indicates that when we aggregate the samples that does not differ too much, the estimator tends to converge faster,

---

especially when n >> (A)), where we have plethora of samples in addition to the samples for the base cohort.
With the same choice of smin and 5 as the above theorem, the probability for consistent variable selection can approach one, with less stringent requirement in the size of S, The above theorem implies that our proposed method can identify both baseline and differential regulatory effects, among all of the K networks with a sufficiently large probability, not only in terms of the set of true signals but also the sign of signals. In the case of gene regulatory networks, for instance, our method could correctly distinguish the up and down regulations between genes.
Next, we will present a theorem for the coefficient of determination and coefficient of cause defined in Section 3.2 Denote the statisties R and C2 calculated with real parameters as, respectively,
(() R
= 1
| (M) -
- (1) 1
Co
2 -/ IB)
(1-111
Then derive the following properties.
Theorem 4.2. (Coefficients of DeterminationnCause Consistency) Under the same conditions as in
Theorem 4.1 with hn vk(k) /- log + 1, we have that,
1. With probability at least p225min + 5/K),
12tk) R (()))
S fd log(p/0) n (A33/4
|| IT 2 V log(d/ (oddm)mn)))) (1 + ll2ht):
2. With probability at least p(8min + 5/
|Siggn (d log(p/8) V || Irlimax V Vlog(d/0min) (1 ||B]|] v V1).
28))
Coi
(())
As considered in Theorem 4.1 we set 6mén x 6/(KK) = enn" Then hn X Vn(k), and D [12(() | x max [S, (2-400//2 ( (max (S:1) -20-1 which would vanish if ISiSP
ISiSP
[Si] S (4401))/2 Similarly, E (02(6)
- 22(1) max |5 (1-m(1140)// 40)/4 and will converge ISVSP
mar 18iSP
to 0 when max [S|| S n (40-1)/4 1)/4, From the above discussion. both R² R²,C0 for all the nodes in the SSiSP
whole system is & consistent. It is worth noting that, although our algorithm fits model for each responder, estimate of each node's explanatory power as a driver is also well controlled over the system, but with a slightly larger bound shown in the above theorem.
5 Simulation Study
We examine the performance of NetANOVA by simulating data of sample size 200, 500 and 1000 from each of three pertinent DCGs, consisting of 1000 endogenous variables. However, only 50 endogenous variables are involved with causality, with each regulated by 3 others on average. Each of the three networks has unique causal effects. First two networks share five causal effects but with opposite signs to that of the third one, resulting in the first two networks having 15 different effects from the third one and the first network have 20 different effects from the second one. The size for causal effects is taken from uniform distribution over [-0.8, -0.3] U[0.3, 0.8). The IVs are generated from a multinomial distribution with 3 outcomes 0, 1, 2. with probabilities 0.25, 0.5, 0.25 respectively. The disturbance errоrs are sampled independently from N (0, 0.13).
For each sample size, we simulated 100 data sets and applied NetANOVA to each by bootstrapping 100 times to calculate p values and construct Receiver Operating Characteristic (ROC) curves as shown in Figure 2 For causal effects of all three DCGs, Figure a shows an almost perfect area under curve (AUC) for each sample size. When comparing causal effects of the first two DCGs vs. the baseline one, i.e., the third DCG, the AUC bottoms 0.867 with sample size at 200, but reaches one with sample size at 1000, showing excellent performance of NetANOVA.
7

---

Three DCGs
(b) DCG vs. DCG
DCG vs. DCG
(d) DCG DCG
0.75
Senaltiily
Sensltuity
Senelliity
Sanslfidity
0.000
200
0.50
220
8.50
220 1000
0.50
200
1000
0.25
0.25
0.25
0.25
0.00
0.00
0.00
0.50
0.00
0.20
0.50
0.75
1.00
0.00
0.25
0.50
0.75
1.00
0.00
0.25
0.50
75
1.50
0.00
0.75
0.50
0.75
1.0
1-Specflcty
1-Speaftety
1-Spoffiitt
1-Speallcty
Figure ROC curves of NetANOVA on simulated data. (a) ROC curves of causal effects constructed in all three networks; (b) ROC curves of deviated effects of the first network from the third one; (c) ROC curves of deviated effects of the second network from the third one; (d) ROC curves of deviated effects of the first network from the second one.
Taking one simulated dataset with sample size at 500, we evaluated our proposed statistics R2 and C2 by plotting estimated statistics against the statistics calculated with true causal effects Figure For each of the three networks, we observe linear trends with slope almost one for both Rº and C statistics, indicating their excellent performance.
(b)
0.6
+1.01 X3,
0.8
+1.01 X1.
+1.01 X1,
0.4-
2
netl netz
0.4-
nett netz
0.2-
R
C
Figure 3: Statistics R2 and C2 on simulated data with sample size 500. We plot each statistic (n²² or C2) vs. its value calculated with true causal effects (R) or CB) for each network. We simply label the three networks as and 3. Shown in the top-left corners are the fitted linear regression models between the pairs as well as the corresponding coefficients of determination.
also took one simulated data with sample size 500 and conducted our proposed correspondence analysis to identify drivers and responders that show important dynamic causality across the three DCGs. When comparing DCG vs. DCG III, Figure 4 shows a driver-responder pair (14, 45) with different causal effects, which is verified by observing 14 significantly up regulates 45 in DCG I vs. DCG III in Figure 4 (blue connections). However, we also see the driver 18 and responder 45 stay opposite the y-axis, because 18 significantly down regulates 45 implying in DCG DCG III in Figure 4a. a. On aaxiss, we see pairs (12,13) and (12,7) which correspond to another cluster of causal networks shown in the right of Figure a.
Real Data Analysis
applied NetANOVA to investigate the gene regulatory networks of lung tissues healthy individuals and patients with lung adenocarcinoma (LUAD) lung squamous cell carcinoma (LUSC). We obtained transcriptomic and genotypic data of healthy lungs (n 482) from the Genotype-Tisme Expression (GTEx) project [6] and of both LUAD (n 485) and LUSC (n 406) from The Cancer

---

(a) Part of the True Networks
25
13
49
3
19
45
18
(b) DCG vs. DCG III
(c) DCG vs. DCG
(d) DCG I vs. DCG II
I3
Dimerdoonn 2.0.2%
Dinerriio (18119
Dimension 2116 (16.28%)
13
17
23
3 [25] 21 24
14
36
11
19
mension (21.7%)
Dimension (36%)
Dimension 1 (20.7%)
Figure 4: Correspondence analysis of NetANOVA on simulated data with sample size 500. (a) Plot of some deviated causal effects among the three true networks (shown in Figure (6) with DCG vs. DCG III, DCG II vs. DCG III, and DCG vs. DCG II in blue, red, and green, respectively. There are
total of 1000 endogenous variables with the first 50 involved in causal relations. (b)-(d) Plots for correspondence analysis with drivers in red and responders in blue.
Genome Atlas (TCGA) project [22). After pre-processing, there are 15,135 genes and 427,820 single nucleotide polymorphisms (SNPs) being shared by three cohorts. Cis-eQTL mapping identified 7059 genes with at least significant SNPs inside their genetic regions (p (p-value<: 0.05), i.e., valid IVs. We bootstrapped the data 100 times to assess the significance of all effects with results shown in Table 1
Table 1: Summary of causal effects identified by NetANOVA. Shown in the columns are the results from the original data, different bootstrap cutoffs (80% 100%), and adjusted by Benjamini-HHochbee adjustment (BH), respectively.
Type of Effects
Original
80%
90%
95%
100%
BH
Healthy Tissue
79833
16594
11481
8760
4602
3165
LUAD Healthy Tissue
13711
1185
848
670
458
296
LUSC vs Healthy Tissue
13104
1385
980
768
477
274
LUSC vs LUAD
18139
1615
976
665
289
38
Controlling adjusted p-value at 0.1, we have the largest subnetwork bearing differential structures shown in Figure 5]a, which is verified via STRING [21] as in Figure 5... STRING reports protein-protein interaction (PPI) enrichment with p-value 10-16 16 implying significant causal effects between the genes shown in Figure Sla We identified many previous validated relationships, such as the connected pair RPS14 and RACK1 (GNB2L1) which were experimentally verified by [ Our results suggest new findings, particularly causal relationship rather than mere association, e.g., deviated regulations between ARL10 and CLTB from healthy lung tissues in both LUAD and LUSC. Our correspondence analysis of these deviated effects, shown in Figure c-e. confirms these important driver-responder pairs, and the calculated statistics R2 and C2 Table 2 also show such deviation between the three cohorts.

---

(a) Constructed Largest Subnetwork
(b) Enrichment in STRING
LUAD vs. Healthy
(d) LUSC vs. Healthy
(c) LUAD vs. LUSC
Dimereton (29.2%)
Dimsension (28.7%)
25
Dinemsionn (10.4N)
HIGD2A
1R110
[AR100] (LLT
[XHOCT [GAAT4
HISOZA, ELTD
6.0
GR814
[4001]
-15
0.0
Dimension (29.4%)
Dimension 1 (40.7%)
Dimension 1 (511%%)
Figure 5: Partial results of gene regulatory networks for human lung. (a) The largest subnetworks of gene regulation with deviated causal effects between healthy, LUAD, and LUSC lung tissues. Shown in black is the gene regulatory network for healthy lung tissue. Shown in blue, red, green are significantly deviated causal effects of LUAD vs. healthy, LUSC vs. healthy, and LUAD vs. LUSC, respectively. (b) Enrichment of the largest subnetwork in STRING. (c)-(e) Correspondence analysis of the deviated effects between the largest subnetworks with drivers red and responders in blue.
7 Discussion
Though having promising applications, construction and dynamic study of multiple causal networks computationally challenging in practice due to the involved large systems. We develop the algorithm NetANOVA which avoids optimizing full-information objective functions of the whole system and instead takes rimitedinffrrmation objective functions, each for revealing all drivers of single responder. Such a responder- focused method allows to deploy parallel computation in two sequential stages and makes possible be computationally scalable the number of involved responders. With available clusters of computers, we can further take the bootstrap method for the usually infeasible task, i.e., evaluating the significance of the constructed causalities in each network and dynamic causalities across networks. Our theoretical analysis and simulation study demonstrate the utility and efficiency of NetANOVA.
Table 2: Statistics R² and C22 for Important Genes
Cohort
ARL10
CLTB
BIF3E
GRB14
HIGD2A
KHDC1
LUAD
0.90
0.74
0.73
0.58
0.71
0.60
R²
LUSC
0.79
0.79
0.77
0.75
0.54
0.00
Healthy
0.00
0.00
0.56
0.00
0.30
0.00
LUAD
1.41
0.74
1.18
0.00
0.86
0.60
C
LUSC
0.79
0.79
0.75
0.75
0.16
0.00
Healthy
0.00
0.00
0.00
0.00
0.00
With the model complexity of multiple causal networks, especially DCGs that NetANOVA supports, it is a daunting task to visualize them and comprehend each network and variation across multiple networks. We have proposed statistics R² and C² to quantify the contributions of responders and drivers within a causal network, respectively. While comparing these two statistics across networks may help identify responders and drivers involved in network dynamics, we develop a correspondence analysis visualize the key players. While our simulation study and real data analysis show promising results, correspondence analysis could be further developed to realize its full potential.
10

---

Acknowledgments and Disclosure of Funding
This research is partially supported by the Office of Research of Purdue University and the National Cancer Institute under Grant R03 CA235363. The results shown here are in whole or part based upon data generated by the TCGA Research Network: https://www.cancer www.ee....eotttgg...... The Genotype-Tisse Tissue Expression (GTEx) Project was supported by the Common Fund of the Office of the Director of the National Institutes of Health, and by NCI, NHGRI, NHLBI, NIDA, NIMH, and NINDS. The data used for the analyses described in this manuscript were obtained from the GTEx Portal on 10/03/2021 and dbGaP accession number phs000424.v8.ppppppmm v8..4.44.....22. on 7/23/2019.
References
[1] Francisco Acosta-Reyes, Ritam Neupane, Joachim Frank, and Israel S Feráández. The Israeli acute paralysis virus IRES captures host ribosomes by mimicking a ribosomal state with hybrid tRNAs. The EMBO Journal, 38(21):e102226 2019.
[2] Peter Bickel, Yaacov Ritov, and Alexandre B Tsybakov. Simultaneous analysis of Lasso and Dantzig selector. The Annals of Statistics, 37(4):1705-177 2009.
[3] Stephen Borgatti, Ajay Mehra, Daniel J Brass, and Giuseppe Labianca. Network Analysis in the Social Sciences. Science, 323(5916):892-88)) 2009.
[4] Xiaodong Cai, Juan Andrés Bazerque, and Georgios B Giannakis. Inference of Gene Regulatory Networks with Sparse Structural Equation Models Exploiting Genetic Perturbations. PLoS Computational Biology, 95))e 1003068, 2013.
[5] Chen Chen, Min Ren, Min Zhang. and Dabao Zhang. A Two-Stage Penalized Least Squares Method for Constructing Large Systems of Structural Equations. The Journal of Machine Learning Research, 19(1):40-77) 2018.
[6] GTEx Consortium al. The GTEx Consortium atlas of genetic regulatory effects across human tissues. Science, 369(6509) 2020.
[7] Mário Cordeiro, Rui P Sarmento, Pavel Brazdil, and João Gama. Evolving Networks and Social Network Analysis Methods and Techniques. Social Media and Journalism: Trends, Connections, Implications, 101(2), 2018.
[8] Henriette Engelhardt, Hans-Peter Kohler, and Alexia Prskawetz. Causal analysis in population studies. Springer, 2009.
[9] Jianqing Fan and Jinchi Lv. Sure Independence Screening for Ultrahigh Dimensional Feature Space. Journal of the Royal Statistical Society: Series B (Statistical Methodology). 70(5): 241-911, 2008.
[10] Ilias Georgakopoulo--Soare Chengyu Deng, Vikram Agarwal, Candace SY Chan, Jingjing Zhao, Fumitaka Inoue, and Nadav Ahituv. Transcription factor binding site orientation and order are major drivers of gene regulatory activity. Nature communications, 14(1):2333, 2023.
[11] Gene H Golub, Michael Heath, and Grace Wahba. Generalized Cross- Validatan Method for Choosing a Good Ridge Parameter. Technometrics, 21(2):215-223, 1979.
[12] Roger A Horn and Charles R Johnson. Matrix Analysis. Cambridge University Press, 2012.
[13] Yan Li, Dayou Liu, Tengfei Li, and Yungang Zhu. Bayesian differential analysis of gene regulatory networks exploiting genetic perturbations. BMC Bioinformatiis, 21(1):1-13, 2020.
[14] Bing Liu, Alberto de La Fuente, and Ina Hoeschele. Gene Network Inference via Structural Equation Modeling in Genetical Genomics Experiments. Genetics, 178(3): 2008.
[15] Judea Pearl. Causality. Cambridge university press, 2009.
[16] Min Ren and Dabao Zhang. Differential Analysis of Directed Networks. arXiv preprint
arXiv:1807.10173, 2018.

---

[17] Phillippe Rigollet and Jan- Christian Hütter. High Dimensional Statistics. Lecture notes for course 188997, 8 81:8814, 2015.
[18] Da Ruan, Alastair Young, and Giovanni Montana. Differential analysis of biological networks. BMC Bioinformatics, 16(1):1-13, 2015.
[19] Mark Rudelson, Roman Vershynin, et al. Hanson-Wright inequality and sub- gaussian concen- tration. Electronic Communications in Probability, 18, 2013.
[20] Peter Schmidt. Econometrics. New York, Marcel Dekker, 1976.
[21] Damian Szklarezyk, Annika L Gable, David Lyon, Alexander Junge, Stefan Wyder, Jaime Huerta- Huerta-Cepas. Milan Simonovic, Nadezhda T Doncheva, John Morris, Peer Bork, al. String v v11: protein-pooted association networks with increased coverage, supporting functional discovery genome-wide experimental datasets. Nucleic Acids Research, 47601):D07--D613, 2019.
[22] John N Weinstein, Eric A Collisson, Gordon B Mills, Kenna R Mills Shaw, Brad Ozenberger, Kyle Ellrott, Ilya Shmulevich, Chris Sander, Joshua M Stuart, Cancer Genome Atlas Research Network, et al. The Cancer Genome Atlas Pan-Cancer analysis project. Nature Genetics, 45 10):
[23] James West, Ginestra Bianconi, Simone Severini, and Andrew E Teschendorff. Differential network entropy reveals cancer system hallmarks. Scientific Reports, 2, 2012.
[24] Peng Zhao and Bin Yu. On Model Selection Consistency of Lasso. The Journal of Machine Learning Research, 1.2541-2563, 2006.
[25] Xin Zhou and Xiaodong Cai. Inference of differential gene regulatory networks based on gene expression and genetic perturbation data. Bioinformatics, 36(1):197-204, 2020.
[26] Hui Zou, The Adaptive Lasso and Its Oracle Properties. Journal of the American statistical association, (4114776)1118-442)) 2006.
12


Reviews:
[{"review_id": "aTmX2JhIsN", "reviewer": "Reviewer_mopu", "summary": "This paper proposes an algorithm, called NetANOVA, for constructing a unified structural model for multiple causal networks with cycles. The algorithm is designed for parallel computation and is scalable to data size and network complexity. It is able to infer causal networks beyond directed acyclic graphs (DAGs), such as directed cyclic graphs (DCGs), which are demanded to depict gene regulatory networks. The paper provides theoretical justification for the algorithm and demonstrates its feasibility and promise with a large-scale simulation study and a real data analysis to compare gene regulatory networks in healthy lung tissues and lung tissues with two types of cancer.", "questions": "* The word \"variational\" is used frequently. However, this word usually refers to well-known methods/concepts in machine learning, such as variational inference. Is this word choice for \"variational causality\" validated?\n* Regarding the weakness, how would the existence of cycles affect identification of causal effects using IV? Since most results on IV are either discussed under the potential outcome framework, or discussed using do-calculus when the underlying causal structure can be described by DAGs.\n", "limitations": "The authors do not explicitly discuss the limitations of this work.", "rating": 4, "confidence": 3, "soundness": 2, "presentation": 2, "contribution": 2, "strengths": "* It proposes a method to unify multiple cyclic graphs with a single structural model and also accommodate their disparities, which is a promising method for analyzing gene regulatory networks where feedback loops are regularly encountered.\n* The algorithm is designed for parallel computation and is scalable to data size and network complexity.", "weaknesses": "* Instrumental variables-based identification of causal effect are proposed for DAGs. There is a lack of discussion on how this extends to causal structures that contains cycles. It is unclear to me if the identification results in 2.3 is correct when cycles exists.\n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "“Instrumental variables-based identification of causal effect are proposed for DAGs. There is a lack of discussion on how this extends to causal structures that contains cycles. It is unclear to me if the identification results in 2.3 is correct when cycles exists.”  \n\nWe agree that IV-based methods are mainly proposed for DAGs. However, the econometrics field has observed much development in systems of equations, including model identification, under the name \"simultaneous equation models [20]. Our specification of model identification is adopted from the rank condition developed for “simultaneous equation models” but in terms of instrumental variables. A rather recent work on identification can be found in the work by Matzkin (2008).  \n\nMatzkin, Rosa L. \"Identification in nonparametric simultaneous equations models.\" Econometrica 76.5 (2008): 945-978.  \n\n“The word \"variational\" is used frequently. However, this word usually refers to well-known methods/concepts in machine learning, such as variational inference. Is this word choice for \"variational causality\" validated?”  \n\nThanks for highlighting the potential ambiguity. In our manuscript, we employed the term “variational” to refer to variation/deviation across multiple networks. To avoid any confusion and ensure clarity, we will replace “variational” with “perturbational” in the revised version of our paper.   "}, {"review_id": "jWxXmXMaJ4", "reviewer": "Reviewer_ufKZ", "summary": "This paper presents a unified structural model that describes multiple DCGs in one model and develops a limited-information-based method to simultaneously infer networks and their disparities. Furthermore, it provides robust non-asymptotic theoretical properties. And it is applied to synthetic and real datasets to show its performance.\n", "questions": "See the weakness section. We do not have any further technical questions.", "limitations": "We do not see any potential negative societal impact of this paper.", "rating": 5, "confidence": 4, "soundness": 3, "presentation": 2, "contribution": 3, "strengths": "The paper proposes a new model that manages to describes multiple DCGs within a single model, and further performs a comprehensive study on it, including algorithmic development, theoretical analysis and synthetic and real-data experiments.\n", "weaknesses": "The paper proposes such a model that describes multiple DCGs by stacking K-network-specific structural models in matrix forms with one network chosen as the baseline. However, we fail to see the benefits of proposing such a \"unified\" model, nor any mathematical novelty in model construction. In the problem setup, there is no relation between or information shared by K networks. And in algorithm development, the algorithm is almost like dealing with network-specific parameters independently. In experiments, we did not see any benefits this approach manifests.\n\nDue to this straightforward way of construction, it can be foreseen that the algorithmic development and theoretical analyses naturally follow without particular technical issues. For algorithm development, it seems more or less doing network-specific inference and then stacking results. The theoretical results are mainly these from high-dimensional statistics (specifically l1-related theory on variable screening and selection). Considering the almost independent treatment to each network in algorithm, the theoretical results naturally follow by repeating these theories from analysis of structural models; thus, the introduction of these high-dim stats tools may not be highly appreciated.\n\nThe simulation study and real-data applications only show the \"success\" application of the proposed model. Without the basic comparison to the state-of-the-art or the independent treatments to each network using classical structural model methods, we have no idea on how good the performance of the proposed method is and what additional benefits we gain by modelling multiple networks via a single model.\n\nThe literature review in Introduction is not good, which takes less than 1/3 space of Introduction. It is not a critical review, and it basically lists 6 papers (line 24-38) and tells what they did. These work are neither being explained on their relations to this paper nor commented on their strengths/weakness/debates. Just a list of work, which however is not sufficient to cover important work that closely relates to methods and theories in this paper. What's worse, the sentence at lien 25-26 reads like something irrelevant copied from other biological papers.\n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "“… fail to see the benefits of proposing such a \"unified\" model, nor any mathematical novelty in model construction. … no relation between or information shared by K networks.”  \n\nWe disagree with the reviewer’s claims on the benefit and novelty. To the best of our knowledge, we are the first to study and compare multiple causal networks, which allows us to address two types of problems: i) Deviation of other networks to a single baseline (e.g., network of the normal/healthy population); ii) Variation (i.e., difference from each other) across multiple networks.   \n\nOur method is novel in its strategically unifying multiple causal networks, shown from Line 86 to Line 95, which allows developing the algorithm NetANOVA to directly identify and estimate deviated causal effects when compared to a single baseline. We are also the first to explicitly define coefficients of determination and coefficients of cause for causal networks. Our proposed correspondence analysis of causal effects across multiple networks is also new with its defined summary statistics which take advantage of parallel computation to obtain bootstrap results.  \n\nThe theoretical benefit of using a unified model is further demonstrated in theorem 4.1. As shown in Theorem 4.1, the error bound of each hat{β}\\_i is related to n=sum_{k=1}^K n^{(k)}  due to the deal with the unified model. Instead, an independent construction of the k-th network will result in an error bound related to n^{(k)}.  \n\nPlease be aware that, i) the reparametrization in (3) requires a pooled model for all involved networks; ii) with K-th network as the baseline, the parameters in β\\_i^{(k)} show the difference between k-th network and the baseline; iii) although the algorithm allows heterogeneous instrumental variables across different networks and therefore may conduct the first stage independently, the second stage has to identify and estimate the deviated causal effects in hat{β}_i which are defined across networks. That is, at this key step, the unified model is directly constructed at Stage 2.3 of Algorithm 1. More details on the second stage can be found under “The inference stage” starting at Line 140.  \n\n“… the algorithmic development and theoretical analyses naturally follow without particular technical issues. … thus, the introduction of these high-dim stats tools may not be highly appreciated.”  \n\nWe disagree with the reviewer’s claim on the theoretical analysis (please refer to our above rebuttal on algorithmic development). For our first-ever algorithm that directly detects and compares multiple causal networks, the error bounds of the corresponding parameters, especially  {hat{β}\\_i, i=1,2,...,p}, is not a straightforward extension of the results from individual networks. As shown in Theorem 4.1, the error bound of each hat{β}\\_i is related to n=sum_{k=1}^K n^{(k)} due to the deal with the unified model. Instead, an independent construction of the k-th network will result in an error bound related to n^{(k)}.   \n\nTechnically, theoretical analysis of causal networks is challenged by inherent endogeneity and henceforth developed two-stage algorithms. Unifying multiple networks to enable direct comparison of multiple networks indeed makes it even more challenging in controlling error bounds while taking advantage of datasets from multiple populations.  \n\nWe focus our research on studying the deviation of other networks to a single baseline (e.g., network of the normal/healthy population), the first type which is popular in disease studies, especially cancer studies. As we stated in Line 89, our algorithm can be extended to study variation (i.e., differences from each other) across multiple networks. However, the parameters of a possibly unified model are involved with parameters from all networks, making it much more difficult (so further efforts are demanded) to develop an appropriate theoretical analysis.   \n\n“The literature review in Introduction is not good, …. It is not a critical review, and it basically lists 6 papers (line 24-38) and tells what they did.”   \n\nAs we are presenting a first-ever algorithm for constructing and analyzing multiple causal networks based on structural models, we focus on our review of related work. The three methods (Liu et al., Cai et al., and Che et al.) are the major development along this direction but for a single network.   \n\n“These work are neither being explained on their relations to this paper nor commented …”   \n\nWe disagree with the reviewer because we described the major idea of each method with comments on its advantages/disadvantages.\n“What's worse, the sentence at lien 25-26 reads like something irrelevant copied from other biological papers.”   \n\nWe again disagree with the reviewer. We have this sentence here to emphasize that constructing structural models with instrumental variables is mainly developed in the field of bioinformatics because natural instrumental variables are available in studying gene regulatory networks. This sentence leads to the review of the developed methods for structural models in bioinformatics. \n"}, {"review_id": "1ScrJRZmAd", "reviewer": "Reviewer_bLo4", "summary": "The paper introduces NetANOVA, an algorithm designed for parallel computation to construct a unified structural model for multiple causal networks, or DCFs. NetANOVA utilizes analysis of variance (ANOVA) to identify causalities that differ across networks, as well as important drivers and responders. It is scalable to large data sizes, computational environments, and model complexities, allowing for efficient analysis of multiple networks.", "questions": "1. In NetANOVA, most computations are matrix products and inversions. Can you briefly estimate the computation cost required when n and k scales?\n2. In Figure 4, are DCG(I/II/III) curated by humans? If so, what is the cost of curating each DCG?", "limitations": "See \"weaknesses.\"", "rating": 5, "confidence": 2, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "1. The paper introduces a two-stage parallelizable algorithm that scales when multiple cohorts exist.\n2. The paper gives solid theoretical results on the consistency for coefficients of determination and cause.", "weaknesses": "1. Experimental results on real data is limited to three DCGs, and it does not contain comparisons with baseline approaches.\n2. The theoretical result does not include the computation complexity of NetANOVA.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "“In NetANOVA, most computations are matrix products and inversions. Can you briefly estimate the computation cost required when n and k scales?”  \n\nAssuming bar{n}=sum_{k=1}^K n^{(k)} the average sample size,  we can break down the computational complexity as follows: The complexity associated with ISIS is O(K bar{n} p), complexity with ridge regression is O(K bar{n}^3 ), projection’s complexity is O(K bar{n}^3 ), complexity introduced by adaptive lasso is O(K^2 bar{n} p). Collectively, the computational cost for each node is O(K^2 bar{n} p) + O(K bar{n}^3).  \n\n“In Figure 4, are DCG(I/II/III) curated by humans? If so, what is the cost of curating each DCG?”  \n\nFor each sample size (200, 500, or 1000), we have randomly simulated 100 sets of networks (with each set including 3 DCGs) as described in Lines 229-236,  and accordingly simulated one dataset from each set. Figure 4.a shows a part of one simulated set with a sample size 500 (the whole network is shown in Figure 6). We applied our algorithm to the data simulated from this set and conducted correspondence analysis of the results.  Figure 4.b-d show the plot of the correspondence analysis. Figure 4.a was obtained via the freely available software Cytoscape. and Figure 4.b-d was from R with packages ggplot and ggrepel.\n"}, {"review_id": "37fn9CqZxr", "reviewer": "Reviewer_Wm1e", "summary": "In this paper, the authors proposed NetANOVA, an algorithm that simultaneously constructs multiple causal networks and infer their disparities. Theoretical justification of the proposed method is also derived. The paper further proposes measures for variable’s contribution to receivers and responders. Overall, the problem seems solid and experiment results reveal its effectiveness. ", "questions": "How is the proposed method compared with related works?\n\nThe model is derived on the basis of several assumptions. While in experiments, do those assumptions really hold on the datasets?\n\nDo the authors conduct runtime analysis of their proposed method? ", "limitations": "It will be good if the authors can provide a formal definition on the problem they are trying to solve.", "rating": 6, "confidence": 2, "soundness": 3, "presentation": 2, "contribution": 2, "strengths": "The task this paper is tackling is of interest and importance. The proposed method is unified on handling multiple causal networks simultaneously and is scalable in various aspects.\n\nThe derivation of the method is clear and reasonable. The authors discuss annotations and assumptions of their proposed method.  They show theoretical justification of the algorithm, which consolidates their work. \n", "weaknesses": "The paper is not really good to read. Section 2 starts deriving the method directly without further explanation and formulation of the problem. The authors should consider adding a subsection before the method to formally describe the problem.\n\nThe experiments seem good. However, no related methods are compared. The authors  mentioned several related works in the introduction and they highlighted advantages of the proposed method (scalability in various aspects). The experiments do not discuss how their methods compare with others. \n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "“The paper is not really good to read. Section 2 starts deriving the method directly without further explanation and formulation of the problem. The authors should consider adding a subsection before the method to formally describe the problem.”  \n\nThe main purpose is to construct and compare multiple causal networks to gain knowledge of cross-cohort changes in causal effects as well as commonly shared ones (with a focus on deviated causal effects of other networks from a baseline). We are sorry that we have to make the description concise to include all aspects of the algorithm development as well as the theoretical/empirical evaluation. We have a rather longer version of the manuscript and would be happy to make it available in arXiv.  \n\n“The model is derived on the basis of several assumptions. While in experiments, do those assumptions really hold on the datasets?”  \n\nAll networks are simulated according to Assumption 1. Assumptions taken for theoretical works are mostly weak and widely adopted in the literature. Specifically, in our experiment, both Assumptions B.1 and B.2 put assumptions on the spectrums of the variance and covariance matrix as well as the IV size by Assumption B,1,  and are satisfied by our simulation settings. Assumption B.3 holds because it is about augmenting the design matrix but the order of the eigenvalues remains unchanged in our simulation.  Assumption B.4 requires limited association between drivers and non-drivers, which also holds in our case.    \n\n“Do the authors conduct runtime analysis of their proposed method?”  \n\nFor our simulation study with one core on Rome CPU @ 2.0 GHz, analysis of one simulated dataset (without parallel computation) took 11,001s,  23,508s,  and 47,883s for sample sizes 200, 500, and 1000 respectively. The use of one node with 128 cores in HPC would cut the time to 86s, 183s, and 374s respectively, taking advantage of the parallel nature of the algorithm. \n"}, {"review_id": "Saa6rPg7j1", "reviewer": "Reviewer_Daft", "summary": "This paper proposes a single structural model to simultaneously construct multiple networks so as to identify causalities varying across multiple cohorts while identifying stable ones.  Each causal network is represented via a directed cyclic graph, and the authors propose an analysis of variance (ANOVA) algorithm (NetANOVA) to identify causalities that are different across networks and an SVD based correspondence analysis to identify the key drivers and responders. Theoretical properties of the algorithm for DCG construction are established and simulations and a data example illustrate the working of their methodology. \n", "questions": "Please answer questions mentioned under weaknesses above. \n\nSection 5: why is error variance of 0.1^2 a good choice? Can you provide the corresponding signal-to-noise ratio? \nSection 3.2: Would it help in interpretation to normalize the coefficient of cause as well so that it lies in (0,1] as R^2? \n\nPlease check the paper for typos (e.g. p8, line 257: n_2, n_3 instead?; p5, line 156: responders; p2, line 54: endogenous etc).", "limitations": "I did not find a discussion of limitations of the method. Please include a brief discussion of these in the paper.", "rating": 7, "confidence": 2, "soundness": 3, "presentation": 4, "contribution": 3, "strengths": "The paper deals with an important problem and appears to be a novel contribution as the problem of investigating causality with multiple networks as done here, has not been studied in the literature. It is presented in a very clear manner and is easy to read in most parts. Some technical aspects could be clarified further (as mentioned under Questions below). The simplicity of their framework; theoretical guarantees and computational scalability are the main strengths in my opinion.", "weaknesses": "The svd based technique employed for correspondence analysis to detect key responder drivers pairs, is not clear. It would be helpful if the authors could elaborate further on this while also providing the intuition behind this. Bootstrapping is used to achieve this- what sort of a bootstrap method is used and why? \nIntuitively, I would expect the accuracy of inferred DCGs  to depend on the range of cohort sizes (min of n^k's to max of n^k's) with lower accuracy as the range increases. Is this the case? If yes, is it evident from the theoretical results? ", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "“The svd based technique employed for correspondence analysis to detect key responder drivers pairs, is not clear. It would be helpful if the authors could elaborate further on this while also providing the intuition behind this. Bootstrapping is used to achieve this- what sort of a bootstrap method is used and why?”\n\nIt is intimidating to capture key information on vast information from multiple large networks. So we have proposed the correspondence analysis to reveal important clusters of drivers, responders, or even driver-responder pairs for their similarity/dissimilarity across causal networks. For the possible causal effect of node j on node i, its difference between cohort k and cohort l is estimated as hat{γ}_{ij}^{(k)}-hat{γ} _{ij}^{(l)}, however, this value may be data-sensitive and not stable. Therefore, we take advantage of bootstrap results to estimate its standard deviation and further standardize the difference to construct the matrix Z^{(k,l)} for SVD. With the top singular vectors catching the major variations cross responders/drivers, the correspondence analysis will reveal important clusters.  \n\nFor bootstrap, we took nonparametric bootstrap by randomly sampling the data with replacement, because it is difficult to develop asymptotic distributions on the estimated parameters for the high-dimensional networks.   \n\n“Intuitively, I would expect the accuracy of inferred DCGs to depend on the range of cohort sizes (min of n^k's to max of n^k's) with lower accuracy as the range increases. Is this the case? If yes, is it evident from the theoretical results?”  \n\nIt is a nice point to make. Yes, with fixed K cohorts, the smallest cohort size, which is denoted as n_min, and the sum of all the cohort sizes, which is denoted as n, determine the accuracy of inferred networks and deviated causal effects.  As shown in Theorem 4.1, the error bound is inversely proportional to n and is also proportionally to g_n which, as shown in Line 427, is inversely proportional to sqrt{n_min}. In summary, when n_min or n decreases, we have a larger error bound.   \n\n“Section 5: why is error variance of 0.1^2 a good choice? Can you provide the corresponding signal-to-noise ratio? Section 3.2: Would it help in interpretation to normalize the coefficient of cause as well so that it lies in (0,1] as R^2?”  \n\nWe chose the error variance of 0.1^2 following works on single causal networks [4, 13]. With respect to the signal-to-noise ratio (defined as var(γ_j Y_j )/var(Y_i |Y_{-i}) when node i’s value Y_i has causal contribution γ_j Y_j  from node j), it ranges from .02861 to 9.8534 with average at .7559 and median at .2335 (based on one simulated dataset).   \n\nIt is an interesting point on normalizing the coefficient of cause into [0,1]. One way is to calculate the averaged coefficient of cause via dividing it by the number of responders taking this driver. Alternatively, we may divide the total driver’s causal contributions in terms of variance by the sum of variances of all responders. Both measures should be read together with the total number of responders as a driver with too many responders may be of high interest even if the causal contribution to each responder is low.  \n\n“Please check the paper for typos (e.g. p8, line 257: n_2, n_3 instead?; p5, line 156: responders; p2, line 54: endogenous etc).”  \n\nThanks for pointing that out. We will correct the typos in our revised version.   \n\n“I did not find a discussion of limitations of the method. Please include a brief discussion of these in the paper.”  \n\nAlthough we have developed a limited-information likelihood method to avoid optimizing too many model parameters as the full-information likelihood method does, the proposed method may still be challenged by large K and massive total sample size n when there are too many cohorts to compare. When K is too large, each task in the algorithm (identifying and estimating causal effects for a single responder) has to estimate K(p-1) parameters with an n×(K(p-1)) design matrix, possibly demanding a large amount of memory.   \n\nWe have developed our algorithm for the case to compare all other networks to a single baseline network and provide theoretical analysis. In practice, we may be interested in the deviated effects of each network from the average effects. Our algorithm can be adopted for such a case. However, it is challenging to develop an appropriate theoretical analysis for this case, and deserves further study.\n"}]

Author rebuttal:
We sincerely appreciate all the reviewers for providing constructive comments, which provide us a chance to clarify some confusions and improve the quality of our work. We have carefully been through each comment and done our best to address each. While we have addressed each reviewer’s points in separate rebuttals, here we would like to address a common concern on “Lack of comparison to the state-of-the-art or the independent treatments to each network using classical structural model methods”.  

To the best of our knowledge, our algorithm is the first one developed for multiple causal network analysis and there is no other state-of-the-art method available to compare to.   

Pooling independently constructed networks seems appealing in terms of computation. However, as shown in Theorem 4.1, the error bound of each hat{β}\_i from our algorithm is inversely proportional to n=sum_{k=1}^K n^{(k)}. However, an independent construction of the k-th network will result in an error bound related to n^{(k)}. Furthermore, due to the high-dimensional property, the estimates of many parameters may follow mixture distributions. That is, unlike classical low-dimensional problems presenting asymptotic normal distributions, the estimated values of bootstrapped high-dimensional datasets may be mixed with many zeros and nonzero values. Therefore, even in the case that the computational challenge overtakes statistical efficiency so we have to pool independently constructed networks, it demands further development of appropriate strategies.  

Because of no state-of-the-art method to compare and difficulty in pooling results from independently constructed networks, we focus our simulation studies on evaluating the feasibility of our algorithm, and its efficacy over different sample sizes.


Meta-review: This paper proposes a methodology for estimating causal effects simultaneously in multiple directed cyclic graphs / structural models that are related, e.g., a baseline model and a few changed models. By pooling the data together, the paper shows that better estimation (with gains that depend on just how much the models are related.) Reviewers generally agree that this contribution is worth sharing with the community. However, several changes would significantly improve the paper. First, presentation needs to be tweaked. For instance, explain the model without equations first, and then instantiate it and, for the theory, include all key assumptions as part of the main text, so that the scope of when guarantees hold is clear. Second, the experiments currently do not have a baseline method to compare against. The authors are strongly encouraged to include one, e.g., perhaps even a method that estimates causal effects individually in each case. Ideally, this should be done along with highlighting the statistical and computational gains that the new methodology provides. Lastly, the authors are encouraged to incorporate all the minor suggestions provided by the reviewers, e.g., changing the choice of the word "variational", and better explaining the bootstrapping methodology and its dependence on the cohort sizes.
Decision: Accept (poster)