Interpretable Models for Granger Causality Using Self-explaining Neural Networks
Authors: Ričards Marcinkevičs
Julia E Vogt
Venue: ICLR 2021 Poster (2021)
Keywords: time series
Granger causality
interpretability
inference
neural networks

Abstract: Exploratory analysis of time series data can yield a better understanding of complex dynamical systems. Granger causality is a practical framework for analysing interactions in sequential data, applied in a wide range of domains. In this paper, we propose a novel framework for inferring multivariate Granger causality under nonlinear dynamics based on an extension of self-explaining neural networks. This framework is more interpretable than other neural-network-based techniques for inferring Granger causality, since in addition to relational inference, it also allows detecting signs of Granger-causal effects and inspecting their variability over time. In comprehensive experiments on simulated data, we show that our framework performs on par with several powerful baseline methods at inferring Granger causality and that it achieves better performance at inferring interaction signs. The results suggest that our framework is a viable and more interpretable alternative to sparse-input neural networks for inferring Granger causality.

Published as a conference paper at ICLR 2021
INTERPRETABLE MODELS FOR GRANGER CAUSALITY USING SELF-EXPLAININ NEURAL NETWORKS
Ričards Marcinkevičs Department of Computer Science
Julia E. Vogt
Department of Computer Science ETH Zürich Universitätstrasse 6
ETH Zürich Universitätstrasse 6
8092 Zürich, Switzerland
8092 Zürich, Switzerland julia. www.ggtti...f...o....mmomm eetzz.thcc.. ch
ricards mww.iineevicssinf.....
ABSTRACT
Exploratory analysis of time series data can yield a better understanding of com- plex dynamical systems, Granger causality is a practical framework for analysing interactions in sequential data, applied in a wide range of domains. In this paper, we propose a novel framework for inferring multivariate Granger causality un- der nonlinear dynamics based on an extension of self-explaining neural networks. This framework more interpretable than other neural- nettoll--twwoo-bbeeeee techniques for inferring Granger causality, since in addition to relational inference, it also allows detecting signs of Granger-causal effects and inspecting their variability over time. In comprehensive experiments on simulated data, we show that our framework performs on par with several powerful baseline methods at inferring Granger causality and that it achieves better performance at inferring interaction signs. The results suggest that our framework is a viable and more interpretable alternative sparse-input neural networks for inferring Granger causality.
INTRODUCTION
Granger causality (GC) Granger 1969)
a popular practical approach for the analysis of multivari-
ate time series and has become instrumental in exploratory analysis (McCracken) (0166 in various
disciplines, such as neuroscience (Roebroeck al., 2005), economics (Appiahl 2018), and clima- tology (Charakopoulos et al.) (018)). Recently, the focus of the methodological research has been on inferring GC under nonlinear dynamics (Tank et al.) 2018 Nauta et al.) 2019 Wu et al. 2020 Khanna & Tan 2020 Löwe et al. (220)), causal structures varying across replicates (Löwe et al. 2020). and unobserved confounding (Nuuta et al. 2019 Löwe et al. 2020).
To the best of our knowledge, the latest powerful techniques for inferring GC do not target the effect sign detection (see Section 2.1 for a formal definition) or exploration of effect variability with time and, thus, have limited interpretability. This drawback defeats the purpose of GC analysis as an exploratory statistical tool. some nonlinear interactions, one variable may have an exclusively positive negative effect another if it consistently drives the other variable up down, spectively. Negative and positive causal relationships are common in many real-world systems, for example, gene regulatory networks feature inhibitory effects (Inoue al. 2011) in metabolomics, certain compounds may inhibit promote synthesis of other metabolites (Rinschen et al.) (019).. Differentiating between the two types of interactions would allow inferring and understanding such inhibition and promotion relationships in real-world dynamical systems and would facilitate more comprehensive and insightful exploratory analysis. Therefore, we see need for a framework ca- pable of inferring nonlinear GC which is more amenable to interpretation than previously proposed methods (Tank et al. 2018: Nauta et al., 2019 Khanna Tan 2020). To this end, we introduce novel method for detecting nonlinear multivariate Granger causality that is interpretable, in the sense that it allows detecting effect signs and exploring influences among variables throughout time. The main contributions of the paper
follows:
We extend self-explaining neural network models Allarrz-Mells & Jaakkola 2018) to time series analysis. The resulting autoregressive model, named generalised vector autore-

---

Published as a conference paper at ICLR 2021
gression (GVAR), is interpretable and allows exploring GC relations between variables, signs of Granger-causal effects, and their variability through time.
2. We propose a framework for inferring nonlinear multivariate GC that relies on a GVAR model with igated by finding relationships that are stable across original and time-reversed (Winkler et al. 2016) time series data.
and time-smoothing penalties. Spurious associations are mit-
3. We comprehensively compare the proposed framework and the powerful baseline methods Tank et al. (2018), Nauta et al. (2019), and Khanna & Tan (2020) on a range of synthetic time series datasets with known Granger- relationships. We evaluate the ability of the methods infer the ground truth GC structure and effect signs.
2 BACKGROUND AND RELATED WORK
2.1 GRANGER CAUSALITY
Granger-causal relationships are given by a set of directed dependencies within multivariate time
series. The classical definition of Granger causality is given, for example, by Lütkepohl (2007). Be-
low we define nonlinear multivariate GC, based on the adaptation by Tank et al. (2018). Consider a
time series with p variables: {xt]eg+ {x] ] x2 ... ...x)) E= Assume that causal relationships
between variables are given by the following structural equation model:
1y: x = gi (x (xi(t-1)) 1) x)) + Eú, for 1 sssisp
(1)
where x) 1) is a shorthand notation for x1, x ot are additive innovation terms; and gi(-) are potentially nonlinear functions, specifying how the future values of variable x depend on the past values of x. We then say that variable x² does not Granger-cause variable x denoted as x² A x", if and only if g:(-) is constant in x x1s(t-1))
Depending on the form of the functional relationship gi()), we can also differentiate between positive and negative Granger-causal effects. In this paper, we define the effect sign as follows: if g:(-) is increasing in all x1- then we say that variable x has positive effect on x", if 9:(-) is decreasing in x (()(t-1) then x² has negative effect on x². Note that an effect may be neither positive nor negative. For example, x² can 'contribute' both positively and negatively to the future of x at different delays, or, for instance, the effect of x on x' could depend on another variable.
Granger-causal relationships can be summarised by a directed graph G = (V, referred to as sum- mary graph (Peters et al. 2017). where V {1, ... is set of vertices corresponding to variables, x is set of edges corresponding to Granger-causal relationships. Let A E 10, 1PPP denote the adjacency matrix of
and & = {i,,)
G. The inference problem is then to estimate A from observations {m=)²= where T the length of the time series observed. In practice, we usually fit a time series model that explicitly or implicitly infers dependencies between variables. Consequently, a statistical test for GC is performed. A conventional approach (Lütkepohl (007) used to test for
linear Granger causality is the linear vector autoregression (VAR) (see Appendix
2.2 RELATED WORK
2.2.
TECHNIQUES FOR INFERRING NONLINEAR GRANGER CAUSALITY
Relational inference in time series has been studied extensively in statistics and machine learning. Early techniques for inferring undirected relationships include time-varying dynamic Bayesian net- works (Song et al., 2009) and time-smoothed, regularised logistic regression with time-varying co- efficients (Kolar et al.] 2010). Recent approaches to inferring Granger causal relationships leverage the expressive power of neural networks (Montalto al. 2015 Wang al. 2018 Tank al. 2018 Nauta et al. 2019 Khanna & Tan 2020 Wu et al. 2020 Löwe et al. (020) and are often based on regularised autoregressive models, reminiscent of the Lasso Granger method (Amold et al. 2007).
2

---

Published as a conference paper at ICLR 2021
Tank et al. (2018) propose using sparse-input multilayer perceptron (cMLP) and long short-term memory (cLSTM) to model nonlinear autoregressive relationships within time series. Building on sample efficient economy statistical recurrent unit (eSRU) architecture with sparse input layer weights. Nauta et al. (0199) propose temporal causal discovery framework (TCDF) that leverages attention-based convolutional neural networks test B contains further details about these and other relevant methods.
this, Khanna & Tan 2020) introduce
for GC. GC. Appendix
2018 Nauta al. 2019 Khanna Tan 2020) and in Appendix B (Marinazzo et al. 2008 Ren et al. 2020; Montalto et al. 2015 Wang et al., 201 Wu
Approaches discussed above (Tank
al. 2020 Löwe et al., (020) focus almost exclusively on relational inference and do not allow easily interpreting signs of GC effects and their variability through time. In this paper, we propose more interpretable inference framework, building on self explaining- neuraal networks Melvss & Jaakkola 2018), that, as shown by experiments, performs on par with the techniques described herein.
2.2.2
STABILITY-BASER SELECTION PROCEDURES
The literature stability- -based model selection is abundant (Ben-Hur et al. 2002 Lange et al. 2003 Meinshausen Bühlmann 2010 [Sun et al.) (013). For example, Ben-Hur et al. (2002) propose measuring stability of clustering solutions under perturbations
assess structure in the data and select an appropriate number of clusters, Lange et al. (2003) propose a somewhat similar approach. Meinshausen Bühlmann (2010) introduce the stability selection procedure applicable to a wide range of -igenniinnii problems: their method guides the choice of the amount of regularisation based on the error rate control. [Sun et al. (2013) investigate similar procedure in the context of tuning penalised regression models.
2.2.3 SELP-EXPLAINING NEURAL NETWORKS
2018) introduce self-explaining neural networks (SENN) - class of intrinsically interpretable models motivated by explieimess,, faithfulness, and stability properties. A SENN with a link function g()) and interpretable basis concepts h(x) RD - R2 follows the form f(x) = 9 (0(z)) ((==)) ...., (()))))).,(,))))( f(z))
Alvarez-Melis & Jaakkola
where E IR are predictors; and 0() is neural network with & outputs. We refer 0(x) as generalised coefficients for data point æ and use them to 'explain' contributions of individual basis concepts to predictions. In the case of g()) being sum and concepts being raw inputs, Equation simplifies to
f(x) Σ (()))az
(3)
Appendix lists additional properties SENNs need to satisfy, defined by Alvarez-Melis & Jaakkola 2018).
A SENN is trained by minimising the following regulaitt loss function, which balances performance with interpretability:
(f))), y) + XLo (f(x)),
where C, (fx), y) is loss term for the ground classification regression task; 0 is regular- isation parameter; and (o (a=)) = (() (xx) J (=))l1 the gradient penalty, where Jh the Jacobian of h(-) w.r.t This penalty encourages f(-) to be locally linear.
3 METHOD
We propose an extension of SENNs (Alvarez-Melis & Jaakkola 2018) to autoregressive time series modelling, which is essentially a vector autoregression (see Equation i1 in Appendix A with gen- eralised coefficient matrices. We refer to this model as generalised vector autoregression (GVAR). The GVAR model of order K given by
K
X
Vo. (x(-k) XL-k + Et.
(5)
3

---

Published as a conference paper at ICLR 2021
where Vo, R" RPx is neural network parameterised by 0n. For brevity, omit the intercept term here and in following equations. No specific distributional assumptions are made on the additive innovation terms E Lo, (x++x) is a matrix whose components correspond to the generalised coefficients for lag k at time step f. In particular, the component (i. of Vo, (xx-+) corresponds to the influence of x/_x on x, In our implementation, we use K MLPs for Lo (-) with p imput units and p" outputs each, which are then reshaped into an Rpxp matrix. Observe that the model defined in Equation 5 takes on a form of SENN (sce Equation (3) with future time series values the response, past values basis concepts, and sum as link function.
Relationships between variables x x'....''' x" and their variability throughout time can be explored by inspecting generalised coefficient matrices. To mitigate spurious inference in multivariate time ries, we train GVAR by minimising the following penalised loss function with the mini-batch gradi- ent descent:
T Σ R(V:) 1=++1
A
Σ
| - 10:2-111:
P
vell2 (6)
T-K
T
K
T
t=++1
t=++1
where (x:) is single observed replicate of p-variate time series of length T; ã; Ek= &6, (#t-R) Ør-R is the one-step forecast for the 1-th time point by the GVAR model; V is a shorthand notation for the concatenation of generalised coefficient matrices at the 1-th time point:
6xx (xI-K) L 00x- (0I-K+1) Vo, (z)-1) E Kp; (-) is sparsity-induccin penalty term; and 7 2. 0 are regularisation parameters. The loss function (see Equation 6 consists of three terms: (i) the mean squared error (MSE) loss, (ti) a sparsity-inducing regulariser, and (iii) the smoothing penalty term. Note, that in presence of categorically-vaaue variables the MSE term can be replaced with e.g. the cross-entropy loss.
The poarsit--inducing term R(-) is an appropriate penalty on the norm of the generalised coefficient matrices. Examples of possible penalties for the linear VAR are provided in Table 4 in Appendix A These penalties can be easily adapted to the GVAR model. In the current implementation, we employ the elastic-net-style penalty term (ou Hastie, 2005 | + (1 a) | with a = 0.5.
Nicholson et al. 2017 R( =
The smoothing penalty term, given by T- Ei=k+1 | Lellz. is the average norm of the difference between generalised coefficient matrices for two consecutive time points. This penalty term encourages smoothness in the evolution of coefficients w.r.t. time and replaces the gradient penalty Le (f (æ)) from the original formulation of SENN (see Equation 4p. Observe that if the term is constrained to be O, then the GVAR model behaves as a penalised linear VAR on the training data: coefficient matrices are invariant across time steps.
Thus, the proposed penalised loss function (see Equation 6 allows controlling the (i) sparsity and (ii) nonlinearity of in- ferred autoregressive dependencies. As opposed to the related approaches of Tank et al. (2018) and Khanna & Tan 2220)) signs of Granger causal effects and their variability in time can be assessed as well by interpreting matrices V (æ:), for K 1 < T. Figure [ shows a plot of generalised coefficients versus time in a toy linear time series (see Ap- pendix D for details). Observe that for causal relationships, generalised coefficients are large in magnitude, whereas for non- links, coefficients are shrunk towards O. More- over, the signs of coefficients agree with the true interaction signs (a,). We further support these claims with empirical re- sults in Section 4 In addition, we provide an ablation study
chusal
0.75
0.50
150
200
300
400
900
Figure 1: GVAR generalised coeffi- cients inferred for a time series with linear dynamics.
for the loss function in Appendix D
3.1 INFERENCE FRAMEWORK
Once neural networks Vl Ve.. k 1, ..., K, have been trained, we quantify strengths of Granger- causal relationships between variables by aggregating matrices V a, (æx) across all time steps into
4

---

Published as a conference paper at ICLR 2021
summary statistics. We aggregate the obtained generalised coefficients into matrix S E RPxP as follows:
mediank+1sesr ((V (mi) in D. for s.ss S
St max ISASK
Intuitively, Si, are statistics that quantify the strength of the Granger-causal effect of on x using magnitudes of generalised coefficients. We expect St,j to be close to 0 for non-causal relationships and S > 0 if x + x². Note that in practice S is not binary-valued, as opposed to the ground truth adjacency matrix A, which we want to infer, because the outputs of L òb (-) are not shrunk to exact zeros. Therefore, we need a procedure deciding for which variable pairs S.gj are significantly different from 0.
To infer binary matrix of GC relationships, we propose heuristic stability- based procedure that relies time-reversed Granger causality (TRGC) (Haufe et al. [0122] Winkler et al.) 2016). The in- tuition behind time reversal is to compare causality scores obtained from original and time-reversed data: expect relationships to be flipped time-reversed data (Haufe et al., 2012 Winkler al., 2016). Winkler et al. (2016) prove the validity of time reversal for linear finite-order autoregres- sive processes. In work, time reversal is leveraged for inferring stable dependency structures in nonlinear time series.
Algorithm T summarises the proposed stability-based thresholding procedure. During inference, two separate GVAR models are trained: one on the original time series data, and another on time- reversed data (lines 3-4 in Algorithm D.. Consequently, we estimate strengths of GC relationships with these two models, as in Equation 7 and choose threshold for matrix S which yields the highest agreement between thresholded GC strengths estimated on original and time-reversed data (lines 5-9 in Algorithm D. A sequence of Q thresholds, given by & = (f1...., {q), is considered where the i-th threshold an S--quattle of values in S. The agreement between inferred thresholded structures is measured (line 7 in Algorithm[] using balanced accuracy score (Brodersen et al. 2010), denoted by BA (s))), equal-to the average of sensitivity and specificity, to reflect both sensitivity and specificity of the inference results. Other measures can be used for quantifying the agreement, for example, graph similarity scores (Zager & Verghese (008). In this paper, we utilise BA, because considered time series have sparse GC summary graphs and BA weighs positives and negatives equally. In practice, trivial solutions, such as inferring no causal relationships, only self-causal links or all possible causal links, are very stable. The agreement for such solutions is set O. Thus, the procedure assumes that the true causal structure is different from these trivial cases. Figure in Appendix E contains an example of stability-base baseed thresholding applied simulated data.
Algorithm 1: Stability-based thresholding for inferring Granger causality with GVAR.
Input: One replicate of multivariate time series {æ: )==) regularisation parameters À and 7 2 0; model order K 1: sequence E (E1, ..., EQ), S & < E2 < ... < EQ 1.
Output: Estimate À of the adjacency matrix of the GC summary graph.
Let { be the time-reversed version of { i.e. "i.,...,]] = {xT, .... #
Let (X, x) be the elementwise thresholding operator. For each component of X,
T *Xi,.x' X) =1, if [X,,, | 2 X. and 5 (X.,,, x) 0. otherwise.
Train an order K GVAR with parameters and y by minimising loss in Equation 6 on
{- and compute S as in Equation
Train another GVAR on { and compute S as in Equation
for i to Q do
Let fi; = 9e (S) and R; = (5), where 9 (X) denotes the {-quattlee of X.
Evaluate agreement
S = BA (S, ((T)))((TTT) (ST.R.)) + BA (T(S R..)) (T (S.Kk))
end
Let 1" arg maxi and g* = E...
10 Let Á (S, 9- (5)).
return Á.
5

---

Published as a conference paper at ICLR 2021
To summarise, this procedure attempts to find dependency structure that is stable across original and time-reversed data in order to identify significant Granger- relationships. In Section we demonstrate the efficacy of this inference framework. In particular, show that it performs par with previously proposed approaches mentioned in Section 2.2
3.1.1 COMPUTATIONAL COMPLEXITY
Our inference framework differs from the previously proposed eMLP, cLSTM (Tank et al., 2018), TCDF (Nauta et al. 2019), and eSRU (Khanna & Tan] 2020) w.r.t. computational complexity. Men- tioned methods require training p neural networks, one for each variable separately, whereas our inference framework trains 2K neural networks. A clear disadvantage of GVAR is its memory com- plexity: GVAR has many more parameters, since every MLP it trains has p² outputs. Appendix F provides a comparison between training times simulated datasets with p E {4. 15, 203. In prac- tice, for a moderate order and a larger p. we observe that training a GVAR model is faster than a cLSTT and eSRU.
EXPERIMENTS
The purpose of our experiments is twofold: (i) to compare methods in terms of their ability to infer the underlying GC structure; and (it) to compare methods in terms of their ability to detect signs of GC effects. We compare GVAR to 5 baseline techniques: VAR with F-tests for Granger causality" and the Benjamini- Hochbbrrg procedure (Benjamini & Hochberg 1995) for controlling the false
discovery rate (FDR) (at 9 0.05): cMLP and cLSTM (Tank et al. 2018 TCDF (Nauta al.
2019 3 and eSRU (Khanna & Tan 2020 [
We particularly focus on the baselines that, similarly GVAR, leverage yparsity-inducing penalties, namely cMLP, cLSTM, and eSRU. In addition,
provide a comparison with dynamic Bayesian networks (Murphy & Russell 2002) in Appendix [
The code is available in the GitHub repository: https: /ttpss://gittuut com/ 6092467 GVAR.
4.1 INFERRING GRANGER CAUSALITY
We first compare methods w.r.t. their ability to infer GC relationships correctly on two synthetic datasets. We evaluate inferred dependen- cies on each independent replicate/simulation separately against the adjacency matrix the ground truth GC graph, an example is shown in Figure 2 Each method is trained only on sequence. Unless otherwise mentioned, we use accuracy (ACC) and balanced accu- racy (BA) scores evaluate thresholded inference results. For cMLP,, cLSTM, and eSRU, the relevant weight norms are compared to O. For TCDF, thresholding is performed within the framework based the permutation test described by Nauta et al. (2019). For GVAR, thresh- olded matrices are obtained by applying Algorithm[] [ In addition, we look at the continuously- valued inference results: norms of relevant weights, scores, and strengths of GC relationships (see Equation (7). We compare these scores against the true structure using areas un- der receiver operating characteristic (AUROC) and precision-recat (AUPRC) curves. For all evaluation metrics, we only consider off- diagonal elements of adjacency matrices, ignoring self-causal rela- tionships, which are usually the easiest to infer. Note that evalua- tion approach is different from those of Tank et al. (2018) and Khanna & Tan (2020); this partially explains some deviations from their sults. Relevant hyperparameters of all models are tuned to maximise the BA score or AUPRC (if model fails to shrink any weights zeros) by performing a grid search (see Appendix H for details about hyperparameter tuning). In Appendix M we compare the prediction error of all models on held-out data.
Effect
Cause
Figure 2: GC summary
graph adjacency matrix of the Lorenz 96 system with = 20. Dark cells corre- spond the absence of GC relationship: light cells de- note a GC relationship.
As implemented in the statsmodeis library Seabolld & Perktold 2010).
httpa:///iitbbb..oo///aaaaaaaa// http:://gitthbb.com///aaamamma/ _or_[C_o____

---

Published as a conference paper at ICLR 2021
4.1.1 LORENZ 96 MODEL
A standard benchmark for the evaluation of GC inference techniques is the Lorenz 96 model (Lorenz 1995). This continuous time dynamical system in variables is given by the following nonlinear differential equations:
dx'
(x x)- x--1 x² + F. for Sisp.
dt
where x = x=, = xP-1, and x == x²; and F is a forcing constant that, in combination with p. controls the nonlinearity of the system (Tank al. 2018 Karimi & Paul 20T0). As can
be seen from Equation [ the true causal structure is quite sparse. Figure 2 shows the adjacency matrix of the summary graph for this dataset (for other datasets, adjacency matrices are visualised in Appendix G). We numerically simulate R 5 replicates with = 20 variables and T = 500 observations under 2018) and Khanna & Tan (2220. but includes more variables.
10 and F 40. The setting similar to the experiments of Tank et al.
Table summarises the performance of the inference techniques on the Lorenz 96 time series under F == TO and F = 40. For F = 10, all of the methods apart from TCDF are very successful at inferring GC relationships, even linear VAR. On average, GVAR outperforms all baselines, although performance differences are not considerable. For F 40, the inference problem appears to be more difficult AaapeendixJiivestiaate
performance of VAR and GVAR across a range of forcing constant values). In this case, TCDF and cLSTM perform surprisingly poorly, whereas cMLP, eSRU, and GVAR achieve somewhat comparable performance levels. GVAR attains the best combination of accuracy and BA scores, whereas cMLP has the highest AUROC and AUPRC. Thus, on Lorenz 96 data, the performance of GVAR is competitive with the other methods.
Table 1: Performance comparison on the Lorenz 96 model with F = 10 and 40. Inference is performed on each replicate separately, standard deviations (SD) are evaluated across 5 replicates.
Model
ACC(+SD)
BA(+SD)
AUROC(+SD)
AUPRC(+SD)
VAR
0.918(10.012)
0.838(20.016)
0,940(20.016)
0.825(10.029)
cMLP
0.972(+0.005)
0.95((=0.016)
0.963(#0.018)
0.908(+0.049)
cLSTM
0.970(+00010)
0.950(=0.028)
0.958(=0.029)
0.925(+0.050)
10
TCDF
0.871(#0.012)
0.709(=0.044)
0.857(=00027)
0.601(+0.053)
eSRU
0.966(10.011)
0,951(10.021)
0.963(10.020)
0.936(10.034)
GVAR (ours)
0.982(10.003)
0.902(10.006)
0.900(10.001)
0.976(10.016)
VAR
0.864((F000088)
0.585(=0.028)
0.745(+0.047)
0.474(+0.036)
cMLP
0.683(+0.027)
0.80((+0.017)
0.909
0.956(+0.033)
cLSTM
0.844(+0.012)
0.656(#0.037)
0.661(=00038)
0.005)+00063)
40
TCDF
0.775(+0.003)
0.597(=00029)
0.679(+00021)
0.015(+0.050)
eSRU
0.8.7(L00009)
0.886(10.016)
0.934(10.021)
0.834(10.033)
GVAR (ours)
0.945(10.010)
0.885(10.044
0.009)
0.916(L0.024)
4.1.2 SIMULATED FMRI TIME SERIES
Another dataset we consider consists of rich and realistic simulations of blood-oxygen-level- dependent (BOLD) time series (Smith al. 2011) that were generated using the dynamic causal modelling functional magnetic resonance imaging (fMRI) forward model. In these time series, vari- ables represent 'activity' in different spatial regions of interest within the brain. Herein, we consider R = 5 replicates from the simulation no. 3 of the original dataset. These time series contain p = 15 variables and only T 200 observations. The ground truth causal structure is very sparse (see Ap- pendix G Details about hyperaarameter tuning performed for this dataset can be found in Appendix
H.2 This experiment is similar to one presented by Khanna & Tan (2020).
Table 2 provides comparison of the inference techniques. Surprisingly, TCDF outperforms other methods by a considerable margin (cf. Table I.. It is followed by our method that, on average, outperforms cMLP, cLSTM, and eSRU in terms of both AUROC and AUPRC. GVAR attains BA score comparable to cLSTM. Importantly, eSRU fails to shrink any weights to exact zeros, thus, hindering the evaluation of accuracy and balanced accuracy scores (marked as 'NA' in Table This

---

Published as a conference paper at ICLR 2021
experiment demonstrates that the proximal gradient descent (Parikh & Boyd 2014). implemented by eSRU (Khanna & Tan 2020). may fail to shrink any weights to or shrinks all of them, even in relatively simple datasets. cMLP seems provide little improvement over simple VAR w.r.. AUROC or AUPRC. In general, this experiment promisingly shows that GVAR performs on par with the techniques proposed by Tank et al. (2018) ananna & Tan 2020) in more realistic and data-scarce scenario than the Lorenz 96 experiment.
Table 2: Performance comparison simulated IMRI time series, eSRU fails to shrink any weights to exact zeros, therefore, we have omitted accuracy and balanced accuracy score for it.
Model
ACC(+SD)
BA(=DD)
AUROC(#SD)
AUPRC(#SD)
VAR 0.910(=0.006)
0.013(H0.015)
0.615(H0.044)
0.175(#0.054)
cMLP
0.846(+000025)
0.614(E0.0688)
0.616(#0.068)
0.191(+0.058)
cLSTM
0.830( 020)
0.655(+0.0053)
0.234(#0.0588 0.665(=0.051)
TCDF
0.899(+0.023)
0.728(+0.003))
0.802(1.0.01)
0.368(+.0126)
eSRU
NA
NA
0.654(L0.005)
0.190(+00095)
GVAR (ours)
0.806(++0.070)
0.652(+00.055)
0.687(+00066)
0.289(+...116
4.2 INFERRING EFFECT SIGN
So far, we have only considered inferring GC relationships, but not the signs of Granger- causal effects. Such information can yield a better understanding of relations among variables. To this end, we consider the -otke-Voltera model with multiple species (Bacaër (201 provides a definition of the original two-species system), given by the following differential equations:
dx' ax' -3x² jEPa(X))
Σ y nx for siss.
(9)
dt
dy' dt = Syd
x² -y-, for 1 sjsp.
(10)
kEPaly))
where x correspond to population sizes of prey species; yd denote population sizes of predator species; a, B, n, 6 Paly") are sets of Granger-causes of x² and y², respectively. According to Equations 9 and 10 the population size of each prey species x' is driven down by |P((x') predator species (negative effects), whereas each predator species yd is driven up by |Pa(y') prey populations (positive effects).
p 0 are fixed parameters controlling strengths of interactions; and Px(x*)
Tredatoo-gnyy Dynamics
We simulate the multi-species ically. Appendix contains details about simulations and the summary graph of the time series. To infer effect directions, we inspect of median generalised coefficients for trained GVAR models. For cMLP, cLSTM, TCDF, and eSRU, we spect signs of averaged weights in relevant layers. For VAR, examine coefficient signs. For the sake of fair compar- ison, we restrict all models to a maximum lag of K (where applicable). In this experiment, we focus on BA scores for positive (BAp..) and negative (BAneg) relationships. Ap- pendix L provides another example of detecting effect signs with GVAR, a trivial benchmark with linear dynamics.
system numer-
Generatioed Confiic enbs
1.0
15
Predesor Prodaco
Dayy
1000
1500
2000
Figure 3: Simulated two-species Lotka-Volterra time series (top) and generalised coefficients (bot-
tom). Prey have positive effect predators, and vice versa.
Table 3 shows the results for this experiment. Linear VAR does not perform well at inferring the GC structure, however, its coefficient signs are strongly associated with true signs of relationships. cMLP provides a considerable improvement in GC inference, and surprisingly its input weights are informative about the signs of GC effects. cLSTM fails to shrink any of the relevant weights to zero; furthermore, the signs of its weights are not associated with the true signs. Although eSRU performs better than VAR inferring the summary graph, its weights not associated with effect signs
8

---

Published as a conference paper at ICLR 2021
all. TCDF performs poorly in this experiment, failing to infer any relationships apart from self- causation. Our model considerably outperforms all baselines in detecting effect signs, achieving nearly perfect scores: it infers more meaningful and interpretable parameter values than all other models.
These results not surprising, because the baseline methods, apart from linear VAR, rely in- terpreting weights of relevant layers that, in general, do not need to be associated with effect signs and are only informative about the presence or absence of GC interactions. Since the GVAR model follows form of SENNs (see Equation 2). its generalised coefficients shed more light into how the future of the target variable depends on the past of its predictors. This restricted structure is more intelligible and yet is sufficiently flexible to perform on par with sparse-inpus neural networks.
Table 3: Performance comparison on the multi-species Lotka-Volterra system. Next to accuracy and balanced accuracy scores, we evaluate BA scores for detecting positive and negative interactions.
Model
ACC(+SD)
BA(=SD)
papos
BAneg
VAR
0.383((-0.005)
0.635(H0.060)
0.845(+-0.02)
0.781(=0.042)
cMLP
0.825(+000035)
0.834(+0.043)
0.889(+00033)
0.846(+00008)
cLSTM
NA
NA
0.491(10.026)
0.604(+:0002)
TCDF
0.832((10.013)
0.500(10.012)
0.504(10.090) 0.538(10.045)
eSRU
0.703(110.048)
0.755(10.010)
0.501(10.025)
0.650(100078)
GVAR (ours)
0.977(+-0.005)
0.961(+00014)
0.932(+0.027)
0.999(=0.001)
In addition to inferring the summary graph, GVAR allows inspecting variability of generalised coef- ficients. Figure 3 provides example of generalised coefficients inferred for two-species Lotka- Volterra system. Although coefficients vary with time, GVAR consistently infers that the predator population is driven up by prey and the prey population is driven down by predators. For the multi- species system used to produce the quantitative results, inferred coefficients behave similarly (see Figure 12 in Appendix K)
CONCLUSION
In this paper, focused two problems: (i) inferring Granger-causal relationships in multivari- ate time series under nonlinear dynamics and (ii) inferring signs of Granger-causal relationships. We proposed novel framework for GC inference based on autoregressive modelling with self- explaining neural networks and demonstrated that, on simulated data, its performance is promisingly competitive with the related methods Ta al. (2018) and Khanna & Tan (2020). Proximal gra- dient descent employed by cMLP, cLSTM, and eSRU often does not shrink weights to exact zeros and, thus, prevents treating the inference technique as statistical hypothesis test. Our framework mitigates this problem by performing stability- based selection of significant relationships, find- ing a GC structure that is stable on original and time-reversed data. Additionally, proposed GVAR model is more amenable to interpretation, since relationships between variables can be explored by inspecting generalised coefficients, which, as we showed empirically, are more informative than in- put layer weights. To conclude, the proposed model and inference framework are a viable alternative to previous techniques and are better suited for exploratory analysis of multivariate time series data.
In future research, plan thorough investigation of the stability- based thresholding procedure (see Algorithm [ and of rimerreversal for inferring GC. Furthermore, we would like to facilitate more comprehensive comparison with the baselines real-world data sets. It would also be interesting to consider better-informed link functions and basis concepts (see Equation Z Last but not least, we plan to tackle the problem of inferring time-varying GC structures with the introduced framework.
ACKNOWLEDGMENTTS
We thank Djordje Miladinovic and Mark McMahon for valuable discussions and inputs. We also acknowledge Jonas Rothfuss and Kieran Chin-Cheong for their helpful feedback on the manuscript.

---

Published as a conference paper at ICLR 2021
REFERENCES
D. Alvarez-Melis and T. Jaakkola. Towards robust interpretability with self-explaining neural net- works. In Advances in Neural Information Processing Systems 31, pp. 7775-7784. Curran Asso- ciates, Inc., 2018.
M. O. Appiah. Investigating the multivariate Granger causality between energy consumption, eco- nomic growth and CO2 emissions in Ghana. Energy Policy, 112:198-2000 208, 2018.
A. Amold, Y. Liu, and N. Abe. Temporal causal modeling with graphical Granger methods. In Proceedings of the 13th ACM SIGKDD international conference on Knowledge discovery and data mining- - KDD'07. ACM Press, 2007.
N. Bacaër. Lotka, Volterra and the predator-prey system (1920-1926). In A Short History of Math- ematical Population Dynamics, pp. 71-76. Springer London, 2011.
A. Ben-Hur, A. Elisseeff, and I. Guyon. A stability based method for discovering structure in clustered data. Pacific Symposium on Biocomputing, pp. 6-17, 2002.
Y. Benjamini and Y. Hochberg. Controlling the false discovery rate: A practical and powerful approach multiple testing. Journal of the Royal Statistical Society. Series B (Methodological), 57(1)22893300 1995.
K. H. Brodersen, C. S. Ong. K. E. Stephan, and J. M. Buhmann. The balanced accuracy and its posterior distribution. In 2010 20th International Conference Patlern Recognition, pp. 3121- 3124, 2010.
A. K. Charakopoulos, G. A. Katsouli, and T. E. Karakasidis. Dynamics and causalities of atmo- spheric and oceanic data identified by complex networks and Granger causality analysis. Physica A: Statistical Mechamics and its Applications, 495:436-453 2018.
C. W. J. Granger. Investigating causal relations by econometric models and cross-spectral methods. Econometrica, 37(3):424-438, August 1969.
Haufe, V. V. Nikulin, and G. Nolte. Alleviating the influence weak data asymmetries on Granger-Causes Causal analyses. In Latent Variable Analysis and Signal Separation, pp. 25-33. Springer Berlin Heidelberg, 2012.
K. Inoue, A. Doncescu, and H. Nabeshima. Hypothesizing about causal networks with positive and negative effects by meta-level abduction. In Inductive Logic Programming, pp. 114-129. Springer Berlin Heidelberg. 2011.
A. Karimi and M. R. Paul. Extensive chaos in the Lorenz-96 model. Chaos: An interdisciplinary journal of nonlinear science, 20(4):043105, 2010.
S. Khanna and V. Y. F. Tan. Economy statistical recurrent units for inferring nonlinear Granger causality. In International Conference on Learning Representations, 2020.
T. Kipf, E. Fetaya, K.-C. Wang, M. Welling, and R. Zemel. Neural relational inference for interacting systems. In Proceedings of the 35th International Conference on Machine Learning, volume 80, pp. 2688-2697, PMLR, 2018.
M. Kolar, I.. Song, A. Ahmed, and E. P. Xing. Estimating time-varying networks. Annals of Applied Statistics, 4(1):94-123, 03 2010.
Lange, M. L. Braun, Roth, and J. M. Buhmann. Stability-based model selection. In Advances in Neural Information Processing Systems, pp. 633-642, 2003.
E. N. Lorenz. Predictability; problem partly solved. In Seminar Predictability, volume 1, pp- 1-18, Shinfield Park, Reading, 1995.
S. Löwe, D. Madras, R. Zemel, and M. Welling. Amortized causal discovery: Learning to infer causal graphs from time- series data, 2020. arXiv:2006.10833.
H. Lütkepohl. New Introduction to Multiple Time Series Analysis. Springer, 2007.
10

---

Published as a conference paper at ICLR 2021
D. Marinazzo, M. Pellicoro, and S. Stramaglia. Kernll method for nonlinear Granger causality. Physical Review Letters, 100(14), 2008.
J. M. McCracken. Exploratory causal analysis with time series data. Synthesis Lectures on Data Mining and Knowledge Discovery, 8(1):1-147, 2016.
N. Meinshausen and P. Bühlmann. Stability selection. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 72(4):417-473, 2010.
A. Montalto, Stramaglia, L. Faes, G. Tessitore, Prevete, and D. Marinazzo. Neural networks with non-uniform embedding and explicit validation phase to assess Granger causality. Neural Networks, 71:159-171, 2015.
K. Murphy and Russell. Dynamic Bayesian networks: representation, inference and learning. 2002.
M. Nauta, D. Bucur, and C. Seifert. Causal discovery with attention-based convolutional neural networks. Machine Learning and Knowledge Extraction, 1:312-340, 2019.
W. B. Nicholson, D. S. Matteson, and J. Bien. VARX-L: Structured regularization for large vector autoregressions with exogenous variables. International Journal of Forecasting, 33(3):627-651, 2017.
N. Parikh and S. Boyd. Proximal algorithms. Foundations and Trends optimization, 133):177-239, 2014.
J. Peters, D. Janzing, and B. Schölkopf. Causal inference on time series using restricted structural equation models. Advances in Neural Information Processing Systems 26, pp. 154-162. Curran Associates, Inc., 2013.
Peters, D. Janzing, and B. Schölkopf. Elements of Causal Inference Foundations and Learning Algorithms. The MIT Press, 2017.
D. Quesada. dbnR: Dynamic bayesian network learning and inference, 2020. URL https://
CRAN. R- pe-trteect.ooo wtgg R package (v. 0.5.3).
Core Team. R: A language and environment for statistical computing, 2020. URL https:
7/www. R- R-project org/
W. Ren, B. Li, and M. Han. A novel Granger causality method based on HSIC-Lasso for revealing nonlinear relationship between multivariate time series. Physica A: Statistical Mechanics and its Applications, 143:423 2020.
M. M. Rinschen, J. Ivanisevic, M. Giera, and G. Siuzdak. Identification of bioactive metabolites using activity metabolomics. Nature Reviews Molecular Cell Biology, 20(6):33--677, 2019.
A. Roebroeck, E. Formisano, and R. Goebel. Mapping directed influence over the brain using Granger causality and fMRL. Neurolmage, 25(1):230-242, 2005.
S. Seabold and J. Perktold. statsmodels: Econometric and statistical modeling with python. In 9th Python in Science Conference, 2010.
S. M. Smith, K. L. Miller, G. Salimi-Khorshidi, M. Webster, C. F. Beckmann, T. E. Nichols, J. D. Ramsey, and M. W. Woolrich. Network modelling methods for FMRL. Neurolmage, 54(2):875- 891, 2011.
L. Song, M. Kolar, and E. Xing. Time-varying dynamic Bayesian networks. In Advances in Neural Information Processing Systems 22, pp. 1732-1740. Curran Associates, Inc., 2009.
W. Sun, J. Wang, and Y. Fang. Consistent selection of tuning parameters via variable selection stability. Journal of Machine Learning Research, 14(1)3319-3440, 2013.
A. Tank, I. Covert, N. Foti, A. Shojaie, and Fox. Neural Granger causality for nonlinear time series, 2018. arxiv: 1 05842.
11

---

Published as a conference paper at ICLR 2021
1. Tsamardinos, L. E. Brown, and C. R. Aliferis. The max-min hill-climbing Bayesian network structure learning algorithm. Machine Learning, 65(1):31-78, 2006.
Y. Wang, K. Lin, Y. Qi, Q. Lian, S. Feng, Z. Wu, and G. Pan. Estimating brain connectivity with varying-length time lags using recurrent neural network. IEEE Transactions on Biomedical Engineering, 155991993,, 2018.
I. Winkler, D. Panknin, D. Bartz, K.-R. Muller, and S. Haufe. Validity of time reversal for testing Granger causality. IEEE Transactions on Signal Processing, 64(11):2746-2760, 2016.
Wu, T. Breuel, M. Skuhersky, and J. Kautz. Discovering nonlinear relations with minimum predictive information regularization, 2020. arXiv:2001.0188.
H. Xing- Cheg, Zheng, T. Lei, and S. Li-Ping. Research structure learning of dynamic Bayesian networks by particle swarm optimization. In 2007 IEEE Symposium on Artificial Life, pp. 85-91, 2007.
A. Zager and C. Verghese. Graph similarity scoring and matching. Applied Mathematics Letters, 21(1):86-94, 2008.
H. Zou and T. Hastie. Regularization and variable selection via the elastic net. Journal of the Royal Statistical Society: Series (Statistical Methodology), 67(2):301-322 2005.
12

---

Published as a conference paper at ICLR 2021
A LINEAR VECTOR AUTOREGRESSION
Linear vector autoregression (VAR) (Litkepohl 2007 is time series model conventionally used test for Granger causality (see Section 2.1p. linear:
VAR assumes that functions gi(-) in Equation are
X2 =
Σ
v x + E::
(11)
where v E R is the intercept vector; 1 E RPxx are coefficient matrices; and E, ~ Ni (0, E Gaussian innovation terms. Parameter K the order of the VAR model and determines the maximum lag at which Granger-causal interactions occur. In VAR, Granger causality is defined by zero constraints on the coefficients, in particular. x' does not Granger-cause x² if and only if. for all lags h E {1,2, ...., K}, (V*),, O. These constraints can be tested by performing, for example, F-test or Wald test.
Usually VAR model is fitted using multivariate least squares. In high-dimensional time series, regularisation can be introduced to avoid inferring spurious associations. Table 4 shows various sparsity-inducing penalties for linear VAR model of order K (see Equation [1p. described by Nicholson et al. (2017). Different penalties induce different sparsity patterns in coefficient matrices These penalties can be adapted to the GVAR model for the sparsity-inuucing term R(.) in Equation
11, 422 2... Vg.
Table 4: Various sparsity-inducing penalty terms, described by Nicholson al. (017). for lin- ear VAR of order K. Herein, Y 1 L2
xx] E REXRP (cf. Equation and Øg. Different penalties induce different sparsity patterns coef-
MA:K [ L k+
ficient matrices.
Model Structure
Penalty
Basic Lasso Elastic net Lag group Componentwise
a | a)
ΣK | 12 E (0,1)
ΣE=1 Σ
(
((Pa.))all
Elementwise
Σ Σ1=1 Eド=1 -k==
Lag-weighted Lasso
Σ ka Ux
E (0.1)
B INFERRING GRANGER CAUSALITY UNDER NONLINEAR DYNAMICS
Below provide more detailed overview of the related work inferring nonlinear multivariate Granger causality, focusing on the recent machine learning techniques that tackle this problem.
Kernel-based Methods. Kernel-based GC inference techniques provide natural extension of the VAR model, described in Appendix A to nonlinear dynamics. Marinazzo et al. (2008) leverage reproducing kernel Hilbert spaces to infer linear Granger causality an appropriate transformed feature space. Ren et al. (2020) introduce a kernel-based GC inference technique that relies on regularisation Hilbert- Hilbert-Schmide
Schmidt independence criterion (HSIC) Lasso GC.
Neural Networks with Non-uniform Embedding, Montalto et al. (2055) propose neural networks with non-uniform embedding (NUE). Significant Granger causes are identified using the NUE, feature selection procedure. An MIP is 'grown" iteratively by greedily adding lagged predictor components as inputs. Once stopping conditions are satisfied, a predictor time series claimed a input. This significant cause of the target if at least one of its lagged components was added technique is prohibitively costly, especially, in a high-dimensional setting, since requires training and comparing many candidate models. Wang al. (2018) extend the NUE by replacing MLPs with LSTMs.
Neural Granger Causality. Tank al. (2018) propose inferring nonlinear Granger causality using structured multilayer perceptron and long short-term memory with sparse input layer weights, cMLP and cl STM. To infer GC, p models need be trained with each variable as response. cMLP and
13

---

Published as a conference paper at ICLR 2021
cLSTM leverage the group Lasso penalty and proximal gradient descent (Parikh & Boyd 2014 to infer GC relationships from trained input layer weights.
Attention-based Convolutional Neural Networks. Nauta al. (0119) introduce the temporal causal discovery framework (TCDF) that utilises attention-basee convolutional neural networks (CNN). Similarly cMLP and cLSTM (Tank et al. 2018), the TCDF requires training neural network models to forecast each variable. Key distinctions of the TCDF are (i) the choice of the temporal convolutional network architecture over MLPs or LSTMs for time series forecasting and (ūi) the use of the attention mechanism to perform attribution. In addition to the GC inference, the TCDF can detect time delays at which Granger-causal caausll interactions occur. Furthermore, Nauta al. (2019) provide permutation-basee procedure for evaluating variable importance and identifying significant causal links.
Economy Statistical Recurrent Units. Khanna & Tan 2020 propose an approach for inferring
nonlinear Granger causality similar cMLP and cLSTM Tank et al. 20188. Likewise, they penalise norms of weights in some layers to induce sparsity. The key difference from the work of Tank et al.
(2018) is the use statistical recurrent units (SRUs) as a predictive model. Khanna & Tan (2020) propose new sample-efficient architecture economy-SRU (eSRU).
Minimum Predictive Information Regularisation. Wu et al. (2020) adopt an information- theoretic approach to Granger-causal discovery. They introduce learnable corruption, c.g. additive Gaussian noise with learnable variances, for predictor variables and minimise loss function with minimum predictive information regularisation that encourages the corruption of predictor time ries. Similarly to the approaches of Tank et al.) (2018); ; Nauta et al.] (2019): Khanna & Tan 20200) this framework requires training p models separately.
Amortised Causal Discovery & Neural Relational Inference. Kipf et al. (2018) introduce the neu- ral relational inference (NRT) model based on graph neural networks and variational autoencoders. The NRI model disentangles the dynamics and the undirected relational structure represented explic- itly as a discrete latent graph variable. This allows pooling time series data with shared dynamics, but varying relational structures, Löwe et al. (2020) provide natural extension of the NRI model the Granger-causal discovery. They introduce more general framework of the amortised causal discovery wherein time series replicates have varying causal structure, but share dynamics. In con- trast to the previous methods (Tank et al. 2018 Nauta et al. 2019: Khanna Tan 2020 Wu et al.) 2020), which in this setting, have to be retrained separately for each replicate, the NRT is trained on the pooled dataset, leveraging shared dynamics.
C PROPERTIES OF SELF-EXPLAINING NEURAL NETWORKS
As defined by Alvarez-Melis & Jaakkola (2018). g()). 0()), and h(-) in Equation 2 need satisfy:
1. g()) is monotonic and additively separable in its arguments;
2. 8a > 0 with 2; h(()), for all i;
3. (()) is locally sifference--oundee by h(-). i.c. for every wa, there exist 5 >> 0 and L e R
s.t. if ||æ: zo < 8, then | (xx))) < T. [.)))()) hooo))l;
4. {h(=);)<= are interpretable representations of æ;;
5. A is small.
D ABLATION STUDY OF THE LOSS FUNCTION
We inspect hyperparameter tuning results for the GVAR model on Lorenz 96 (see Section 4.1.1
and synthetic fMRI time series (Smith et al. 2011) (see Section 4.1 as an ablation study for the loss function proposed (see Equation (6). Figures and S show heat maps of BA scores (left) and AUPRCs (right) for different values of parameters X and for Lorenz 96 and fMRI datasets, respectively. For the Lorenz 96 system, sparsity-inducing regularisation appears to be particularly important, nevertheless, there is also increase in BA and AUPRC from a moderate smoothing penalty. For fMRI, we observe considerable performance gains from introducing both the sparsity- inducing and smoothing penalty terms. Given the sparsity of the ground truth GC structure and
14

---

Published as a conference paper at ICLR 2021
the number of observations (T 200), these gains not unexpected. During preliminary experiments, we ran grid search across wider ranges of À and y values, however, did not observe further improvements from stronger regularisation. In summary, these results empirically motivate the need for two different forms of regularisation leveraged by the GVAR loss function: the sparsity- inducing and smoothing penalty terms.
6.55
05
08
0.75
6.75
2
1.5
1.5
2.25
2.25
0.5
3.0
2.
0.9
0.006
0,012
0.019
0.045
0.8
0.006
0.032
0.039
0.025
Figure 4: GVAR hyperparameter grid search results for Lorenz 96 time series (under F 40) across 5 values of À E 3.00 and E [0.0, 0.02]. Each cell shows average balanced accuracy (left) and AUPRC (right) across replicates (darker colours correspond to lower performance) for one hyperparameter configuration.
0.30
0.0
0.0
0.62
0.26
0.75
0.75
0.24
1.5
0.10
15
0.20
2.25
2.25
0.18
3.0
0.54
10
5.5
0.02
0.65
0.08
0.1
0.6
0.02
0.05
0.02
0.1
0.14
Figure 5: GVAR hyperparameter grid search results for simulated fMRI time series across values of X E 3.00 and E [0.0. 0.1]. The heat map on the left shows average BA scores, and the heat map on the right average AUPRCs.
15

---

Published as a conference paper at ICLR 2021
E STABILITY-BASED THRESHOLDING: EXAMPLE
Figure shows an example of agreement between dependency structures inferred on original and time-reversed synthetic sequences across range of thresholds (see Algorithm []. In addition, plot the BA score for resulting thresholded matrices evaluated against the true adjacency matrix. As can be seen, the peak of stability agrees with the highest BA achieved. both cases, the procedure described by Algorithm chooses the optimal threshold, which results in the highest agreement with the true dependency structure (unknown the time of inference).
1.0
1.0
0.8
tel. Accuncy
Bal. Accurcy
0.6
0.2
(c) BA ground bruth
Agreement (c)
BA ground bruth
0.0
0.0
15
(a) Lorenz
10.
(b) Multi speciisspeiesssss Lotka Volterra.
Figure 6: Agreement (A) between GC structures inferred the original and time-reversed data across a range thresholds for one simulation the Lorenz 96 (a) and multi-species Lotka-Voltera b) systems. BA score (x) is evaluated against the ground truth adjacency matrix.
F COMPARISON OF TRAINING & INFERENCE TIME
To compare the considered methods in terms of their computational complexity, measure training and inference time across three simulated datasets with E {4, 15, 20) variables and varying time series lengths. This experiment was performed
Intel Core i7-7500U CPU (2.70 GHz 4) with a GeForce GTX 950M GPU. All models were trained for 1000 epochs with a mini-batch size of 64. In each dataset, the same numbers of hidden layers and hidden units were used across all models. When applicable, models were restricted to the same order (K). 5 average training and inference time in seconds with standard deviations. Observe that for the fMRI and Lorenz 96 datasets, GVAR is substantially faster than cLSTMM and eSRU.
Table 5: Average training and inference time, in seconds, for the methods. Inference was performed on time series generated from the linear model (see Appendix L. simulated fMRI time series (see Section 4. 1.2 and the Lorenz 96 system (see Section 4.1.1..
Linear 4, T (p
MMRI
Lorenz 96, F
Model
500,
(p 15, 200,
(p 20, T 500,
K =
K =
K 5)
VAR
0.018(=00001)
0,27(H-.01)
8.5(+:0.6)
cMLP
19.7(13.55)
15.1(1).0)
94.7(11.44
cLSTM
999.2(1177.22)
1763.3(1235.3)
2023.8(111.0)
TCDF
11.4(11.)
24.6(11.7)
50.1(12.1)
<SRU
62.9(+2.6)
258.3(+311.)
671.4(+9.7)
GVAR (ours)
75.3(+8.3)
20.4(+0.7)
197.7(+24.8)
16

---

Published as a conference paper at ICLR 2021
G GC SUMMARY GRAPHS OF SIMULATED TIME SERIES
Effect
Cause
Cause
(a) Lorenz 96.
(b) fMRI.
Cause
(c) Multi-species Lotka-Voltera.
Figure 7: Adjacency matrices of Granger-causal summary graphs for Lorenz 96 (see Section 4.1.1 simulated fMRI (see Bection[4.1.2) and multi-species Lotka-Volterra (see Section 4.2) time series. Dark cells correspond to the absence of GC relationship, i.e. Arg = O; light cells denote GC relationship, i.e. A. 1.
H HYPERPARAMETER TUNING
In our experiments (see Section 4. for all of the inference techniques compared, we searched across a grid of hyperparameters that control the sparsity inferred GC structures. Other hyperparameters were fine-tuned manually. Final results reported in the paper correspond the best hyperparameter configurations. With this testing setup, our goal was to fairly compare best achievable inferential performance of the techniques.
Tables and provide ranges for hyperparameter values considered each experiment. For cMLP and cLSTM (Tank et al. 2018), parameter is the weight of the group Lasso penalty: for TCDF Nauta et al. 2019), significance parameter a is used to decide which potential GC relation- ships significant; eSRU (Khanna Tan 2020) has three different penalties weighted by A1:3- For the stability-based thresholding (see Algorithm in GVAR, we used Q = 20 equally spaced values [0, 1] as sequence 5 For Lorenz 96 and fMRI experiments, grid search results are plotted in Figures and Figure contains GVAR grid search results for the Lotka-Voltera experiment.
'We did not observe high sensitivity of performance w.r.t. as long sufficiently many evenly spaced sparsity levels considered.
17

---

Published as a conference paper at ICLR 2021
H.1 LORENZ 96
Table 6: Hyperparameter values for Lorenz 96 datasets with F = 10 and 40. Herein, K denotes model order (maximum lag). If hyperparameter is not applicable model, the corresponding entry is marked by 'NA'
hidden layers
# hidden units NA
# training epochs NA
Learning rate NA
Mini-batch size NA
Sparsity hyperparam-s
Model
K
VAR
NA
NA
F 10:
À E 2.0];
cMLP
2
50
1,000
1.0e-2
NA
À E [0.0, 1.0]
F 10:
À E [0.1, 0.6]:
cLSTM
NA
2
50
1,000
5.0e-3
NA
À E [0.2, 0.25]
F 10, 40:
TCDF
2
50
1,000
1.0u-2
64
E [0.0,
10, 40:
eSRU
NA
2
10
2,000
5.0e-3
64
A1:3 E 0010
10, 40:
GVAR
2
1,000
1.0u-4
E [0.0,3.0],
Y E [0.0, 0.025
0.0
0.0
09
0:75
0.75
1.5
0.50
1.5
0.7
0.75.
2:25
2.25
5.70
3.0
3.0
0.65
0.3
0.9
0.006
0.012
0.019
0.025
0.0
0.006
0.12
0.019
0.025
Y
Figure 8: GVAR hyperparameter grid search results for Lorenz 96 time series, under F 10, across 5 values of À [0.0,3.0 3.0) and y E [0.0, 0.02]. Each cell shows average balanced accuracy (left) and AUPRC (right) across replicates.
18

---

Published as a conference paper at ICLR 2021
H.2 FMRI
Table 7: Hyperparameter values for simulated fMRI time series.
hidden layers
# hidden units
training epochs
Learning rate
Mini-batch
Sparsity
Model
K
size
hyperparam-s
VAR
NA
NA
NA
NA
NA
NA
cMLP
50
2,000
1.0e-2
NA
A 0.0011 0.75]
cLSTM
NA
50
1,000
1.0e-2
NA
E [0.05, 0.3
TCDF
I
50
2,000
1.0e-3
64
E [0.0, 2.0]
E [0.01, 0.05).
eSRU
NA
2
10
2,000
1.0e-3
64
Àz E [0.01, 0.05),
Àa E [0.01,
À [.00,
GVAR
1
50
1,000
1.0e-4
64
7 0.0,
H.3
LOTKA- VOLTERRA
Table 8: Hyperparameter values for multi-species Lotka-Volterra time series.
hidden layers
# hidden units
training epochs
Learning rate
Mini-batch size
Sparsity hyperparam-s
Model
K
VAR
NA
NA
NA
NA
NA
NA
cMLP
2,000
5.0e-3
NA
X
0.2,, 0.4
cLSTM
NA
2
50
1,000
5.0e-3
NA
À E 0.0, 1.0
TCDF
2
50
2,000
1.0e22
256
E [0.0,2.0 2.0]
À1 E [0.01, 0.05);
<SRU
NA
2
10
2,000
1.0e-3
256
λ2 E 10.01,0.05,
Àa E 1.0
0.0, 1.0),
GVAR
2
50
500
1.0e-4
256
0.00
19

---

Published as a conference paper at ICLR 2021
0.55
AUPHIC
1.0
0.6
0,005
0.005
(a) BA
(b) AUPRC
0.90
0.90
0.75
0.25
1.0
0.65
0.005
0.008
0.002
0.005
0.008
0.01
(c) BA pos
(d) BA
Figure 9: GVAR hyperparameter grid search results for multi-species VolterVolterra time series
across 5 values of X E [0.0, 1.0] and E [0.0, 0.01). Heat maps above show balanced accuracies
AUPRCs (b). and balanced accuracies for positive and negative (d) effects.
I COMPARISON WITH DYNAMIC BAYESIAN NETWORKS
We provide a comparison between GVAR and linear Gaussian dynamic Bayesian networks (DBN). DBNs are classical approach to temporal structure learning (Murphy & Russell 2002). We use R
(R Core Team 2020) package dbnR (Quesada] 2020) to fit DBNs on all datasets considered Sec- structure learning algorithms: the max-min hill-climiie cill-limbiigg (MMHC) (Tsamardinos ct al. 2006) and the particle swarm optimisation (Xing-Chee et al., 2007). Table conttains average balanced accuracies achieved by DBNs and GVAR for inferring the GC structure. Not surprisingly, DBNs outperform GVAR on the time series with linear dynamics, but fail to infer the true structure on Lorenz 96, fMRI, and Lotka- datasets.
tion 4 We
Table Comparison of balanced accuracy scores for GVAR and DBNs. Standard deviations (SD) are taken across 5 independent replicates.
Lorenz 96, F = 10
Lorenz 96,
Model
Linear
MMRI
Lotka-Volterra
= 40
GVAR (ours)
0.988(=0.084)
0.046)00.04) 0.982(=0.006)
0.652(=000055)
0.911(=0.014)
DBN (MMHC)
0.900 0.950(100.088)
0.821(#0.009)
0.687(=0.033)
0.522(=0.003))
0.586(20.337)
(PSO)
0.627(+0.011)
0.514(=0.025)
0.473(+-0.066
0.534(=0.02))
20

---

Published as a conference paper at ICLR 2021
J
THE LORENZ 96 SYSTEM: FURTHER EXPERIMENTS
In addition to the experiments in Section 4.1.1 we examine the performance of VAR and GVAR mod- els across range of forcing constant values F 0,5, 10, 25, 50 for the Lorenz 96 system with p = 20 variables. Figure [0] shows average AUPRCs with bands corresponding to the 95% CI for the mean. appears that for both models, inference is more chal- lenging for lower (< 10) and higher values of F (> 20). This observation in agreement with the sults in Section 4.1.1 where all inference techniques performed worse under F' 40 than under F 10. Note that herein same GVAR hyperparameters were used across all values of F. It is possible that better inferential performance could be achieved with GVAR after comprehensive hyperparameter tuning.
VAR
AUPRI
0.2
10
50
Forcing Constant Value (F)
Figure 10: Inferential performance of GVAR across a range of forcing constant
values.
K THE LOTKA- VOLTERRA SYSTEM
The original voltera system Bacaer 2011) includes only one predator and one prey species, population sizes of which are denoted by and y, respectively. Population dynamics given by the following coupled differential equations:
dx
ax Bxy,
(12)
dt
dy 5yx -yy
(13)
where a, B, 6, > 0 are fixed parameters determining strengths of interactions.
In this paper. we consider multiple species version of the system, given by Equations and 10
1.1, 3 = = 0.2, = 2.75 x 10- Pa(y') 2. p 10, i.. 2p 20 variables in total, with T 2000 observations. Figure II depicts signs of GC effects between variables in a multi-species
in Section 2 We simulate the system under
Pa(x')
with 2p 20 species and parents per variable. We simulate this system numerically by using the Runge- Kutta method We make few adjustments to the state transition equations, in particular: we introduce crrmally-distiibteer innovation terms to make simulated data noisy; during state transitions, we clip all population sizes below O. Figure 12 shows traces of generalised coefficients inferred by GVAR: magnitudes and signs coefficients reflect the true dependency structure.
Negative
Effect
No GC relaticsnhip Positive
Figure 11: Signs of GC relatiosshipss between variables in the Lotka Volterra system given by Equa- tions Lond TO with p 10. First ten columns correspond to prey species, whereas the last ten corre- spond to predators. Each prey species 'hunted' by two predator species, and each predator species 'hunts" two prey species. Similarly to the other experiments, we ignore self-causal relationships.
"Simulations are based on the implementation available at https: thub. com/ amka amm
21

---

Published as a conference paper at ICLR 2021
0.100
Non-causal
Predator - Prey
Generalized Coefficient
Pray
0.050
0.025
0.000
-0.025
-0.050
-0.075
250
500
750
1000
1250
1500
1750
2000
Time
Figure 12: Variability of GVAR generalised coefficients throughout time for a simulation of the multi-species Lotka-Volterra system. Coefficients for Granger non-causal relationships fluctuate around O; for Granger-causal relationships. coefficients are consistently different from O: positive for prey predator interactions and negative for predator - prey.
L EFFECT SIGN DETECTION IN A LINEAR VAR
Herein we provide results for the evaluation of GVAR and our inference framework on a very simple synthetic time series dataset. We simulate time series with = 4 variables and linear interaction dynamics given by the following equations:
01 + E
W₁ = 02Wi-1 + 03X4-1 + Et
(14)
y, a4y,-1 + a5 W + E
Z₁ a6Zt-1 (-1 asYt-1 E
where coefficients a; (1-0.8, 0..2] U [0.2. 0.8)) are sampled independently in each simulation; and E, ~ N (0, 0.16) are additive innovation terms. This is an adapted version of one of artificial datasets described by Peters al. (2013), without instantaneous effects.
The GC summary graph of the system is visualised in Figure 13 It is considerably denser than for the Lorenz 96, fMRI, and Lotka-Volterra time series investigated in Section Similarly to the experiment described in Section [.2, infer GC relationships with the proposed framework and evaluate inference results against the true dependency structure and effect signs. Table 10 contains average performance across 10 simulations achieved by GVAR with hyperparameter values K addition, we provide results for some of the baselines (no systematic hyperparameter tuning was performed for this experiment).
0.2, and 7 0.5. In
Effect
GVAR attains perfect AUROC and AUPRC in all 10 simulations. In some cases, stability-based thresholding fails to recover a completely correct GC structure, nevertheless, average accuracy and balanced accuracy scores are satisfactory. Signs of inferred generalised coef- ficients mostly agree with the ground truth effect signs, as given by coefficients á1.8 in Equation [4]
Cause
Figure 13: The adjacency matrix of the GC summary graph for the model given by Equation 14
Not surprisingly, linear VAR performs the best on this dataset w.r.t. all evaluation metrics. Both cMLP and eSRU successfully infer GC relationships, achieving results comparable to GVAR. However, neither infers effect signs well as GVAR. Thus, similarly to the experiment in Section 4.2 we conclude that generalised coefficients more interpretable than neural network weights leveraged by cMLP, TCDF, and eSRU.
To summarise, this simple experiment serves as a sanity check and shows that our GC inference framework performs reasonably in diw-dimeniional time series with linear dynamics and rela- tively dense GC summary graph (cf. Figure Generally, the method successfully infers both the dependency structure and interaction signs.
22

---

Published as a conference paper at ICLR 2021
Table 10: Performance on synthetic time series with linear dynamics, given by Equation Aver- ages and standard deviations are evaluated across 10 independent simulations. eSRU failed to shrink weights exact Os, therefore, we omit accuracy and BA scores for
VAR
cMLP
TCDF
eSRU
GVAR (ours)
ACC
0.975(H0.038)
0.867(20.085)
0.791(20.056)
NA
0.950(10.067)
BA
0.981(=:0009) 0.900(=0.064)
0.688(#0.169)
NA
0.938(+0.084)
AUROC
1.000(+-0.00) 1.000(+00.00)
0.866(+0.114)
0.972(+0.043)
1.000(+00000)
AUPRC
1.000(+-0.00)
1.000(+00000)
0.81((+0.128)
0.967(+0.046)
1.000C+0.000)
BA
0.995(+-0.014)
0.761(+00206)
0.574(+0.235)
0.613(+0.168)
0.920(+0.158)
B.Acg
0.990(1:0.019)
0.746(10.232)
0.550(10.169)
0.622(10.215)
0.928(10.151)
M PREDICTION ERROR
Herein we evaluate the prediction error models on held-out data. Last 20% of time series points were held out to perform prediction on Lorenz 96, fMRI, and Lotka-Volterra datasets. Root-mean- square error (RMSE) was computed for predictions across independent replicates:
Σ +==
f2 2
RMSE
(15)
P
T
where 2/ is the one-step forecast made by model for the t-th point of the j-th variable, and the length of the held-out time series segment. Table II contains average RMSEs for all models across the considered datasets. In general, RMSES are not associated with the inferential performance the models (cf. tables
and 3.. For example, while TCDF achieves the best inferential performance on IMRI (see Table 2) its prediction error is higher than for cMLP. This 'misalignment" between the prediction error and the consistency of variable selection is not surprising and has been discussed before, e.g. by Meinshausen & Bühlmann (2010).
Table 11: RMSEs of models held-out data. Averages and standard deviations taken across 5 independent replicates.
Lorenz 96
Lorenz 96
Model
[MRT
VolterVolterra
10
F 40
VAR
0.378(10.008)
1.08((10.021)
0..04210.042)
0.202(1000255)
cMLP
0.336(10.030)
0.795(10.011)
0.724(10.037)
0.006(10.016)
cLSTM
0.592(+00013)
0.983(+0.014)
0.874(#0.045)
0.691(+0.035)
TCDF
0.536(+0.015)
1.514(+00.30)
0.879(+0.022)
0.111(#0.016)
<SRU
1.000(+0.010)
1.006(+0.015)
1.000(+00048
0.720(+0.035)
GVAR (ours)
0.572(10.015)
1.00(1000018)
0.966(10.044)
0.119(10.009)
23


Reviews:
[{"review_id": "wXFx_7a0zZt", "reviewer": "AnonReviewer3", "summary": "", "questions": "", "limitations": "", "rating": 6, "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": "# Introduction\n\nThe introduction is solid. My comment here is that there ought to be more discussion on sign detection. You leave even the definition until section 2.1, where a lay reader may struggle to understand what you are talking about if they are merely skimming your paper (as many do). Hence, for readability and for people to actually use this work, I suggest you focus more on sign detection early on and why it is important. You say that it is \"important to differentiate between the two types of interactions\" - well why is that?\n\n# Background and related work\n\n## 2.1\n\nStay clear of using $\\mathbb{V}$ and $\\mathbb{E}$ for vertices and edges. These type-settings are usually reserved for expectations and variances respectively.\n\n# Method\n\nPlease give the form of SENN, in the same form as you've given your GVAR model, i.e. by way equation 4. They are difficult to compare otherwise. Right now what you claim to be an extension of SENN is not easily perceptible by reading section 2.2.3 and looking at equation 4.\n\nThe noise/innovation term $\\boldsymbol{\\epsilon}_t \\in \\mathcal{\\mathbf{0},\\mathbf{1}}$? As standard? Please specify.\n\nAt the bottom of this section you say \"signs of Granger causal effects and their variability in time can be assessed as well by interpreting matrices\" -- it would be helpful if you could provide a toy-example, of matrices, at this point to show us how precisely we are to interpret Psi to get to GC.\n\n## 3.1 \n\nIt is not clear _why_ you \"aggregate the obtained generalised coefficients into matrix $\\mathbf{S}\"?\n\n# Algorithm 1\n\nWhat is $\\boldsymbol \\alpha$?\n\nOn line 7 you introduce the balance accuracy score. What is that? Why is it used in GVAR? Why did you choose that? What are the alternatives? My impression was that BA is the average of recall obtained on each class, in a classification problem. But here you are using it in a much more novel way. It would be very interesting to understand your reasoning and logic here.\n\nOtherwise, well-written algorithm. Easy to follow and understand. My main pointer, as noted, is that you introduce far too many new things in the algorithm rather in the body itself. The algorithm is meant to summarise your method exposition of your method in the body, not introduce new items such as $\\boldsymbol \\alpha$ and BA.\n\n# Experiments\n\nWhat you are doing, from a purely ML point of view, is structure learning. This is very interesting of course, but I note that you have mainly compared to methods specifically focused on causality. It would be helpful/interesting if you also included more vanilla ML methods used for structure learning. \n\nSpelling in sec 4.1 \"releveant\" should be 'relevant'.\n\nIf you can move some of your adjacency matrices into the main paper, from appendix F, that would help the paper a lot. Especially for the Lorenz-96 model. They visually demonstrate what you are trying to achieve, you barely need any legend to go with it. Right, for space issues presumably, they are wasted in the appendices. Please try to move at least one into the body as they a good opticals.\n\nWhy is the inference problem more difficult? F=8 is known to cause deterministic chaos but what about F=10 and 40? Would it not make more sense to evaluate a small ablation study on F from say 0 to 50 in increments of 10? It is not clear to me how computationally heavy GVAR is, hence this may be asking a lot. But since you have not added any such metrics, it is hard to know (computational complexity that is).\n\nYour experiments shadow those of Khanna and Tan pretty well, yet you did not deploy GVAR on the DREAM-3IN SILICO NETWORK INFERENCE CHALLENGE? Why is that? That is the experiment where their model performs the worst, and since they do quite well (at least in Lorenz 96 as shown your paper) I would have thought that a good comparison too. Naturally, you only have so much space, but a comment would be helpful.\n\nPlease correct me if I am wrong here, but not one of your experiments contains an identfiable noise model, except for the model in appendix I? Why is the innovation N(0,0.16) that seems oddly specific? Why not standard white noise N(0,1)? And under this regime your models performs as well as VAR and CMLP by the looks of it.", "extra_scores": {}, "extra_text": {}}, {"review_id": "cre2OujuVuo", "reviewer": "AnonReviewer4", "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": "Interpretable ML approaches are very important to advance a wide variety of fields in healthcare and science. This work makes a significant contribution in this direction by successfully applying SENNs to multi-variate non-linear time series data. The method which the authors named GVAR is generic enough and can be applied to many different problems involving multi-variate time series data thereby making the paper quite impactful and significant. The authors look into GC relationships between variables including the sign of the relationship which they can successfully infer using their approach. The authors provide an array of favorable results across multiple simulated non-linear dynamical systems which corroborate their findings convincingly and also compare their results to a comprehensive set of related and SOTA approaches. There are no obvious inconsistencies or errors that I can see in the paper to the best of my knowledge. The paper is well written and conveys the information clearly enough. An exciting next step is to apply this approach to real data, could be a great addition to this paper but not necessary for publication.  ", "extra_scores": {}, "extra_text": {}}, {"review_id": "bEiOLmRcWok", "reviewer": "AnonReviewer2", "summary": "", "questions": "", "limitations": "", "rating": 6, "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": "Summary: \nThe paper proposes an interpretable method to detect Granger causality (GC) under nonlinear dynamics which can detect signs of Granger-causal effects (positive and negative) and inspect their variability over time. The novelty of this paper is 1 and 2 in the Pros below, and the methods utilized the existing stable frameworks such as a heuristic stability-based procedure that relies on time-reversed Granger causality (TRGC) (Winkler et al., 2016). \n\nReasons for score: \nAlthough the motivation and the experiments were clear and the code was reproducible, the presentations such as in the results were sometimes unclear (below).\nI think this would be a valuable paper in this community, the presentation and investigation may be required to obtain a higher rating. \n\nPros:\n1. The method called generalized vector autoregression (GVAR) generalizes the VAR which is the basis of linear GC methods. The method is interpretable based on a self-explaining neural network (SENN) framework and allows exploring signs of Granger-causal effects and their variability through time.\n2. The method utilizes a framework for inferring nonlinear multivariate GC that relies on a GVAR model with sparsity-inducing and time-smoothing penalties.\n3. The method outperformed the baseline methods on a range of synthetic Lorenz 96 and Lotka–Volterra datasets in inferring the ground truth GC structure and effect signs.\n\nCons:\n1. The detailed explanation and the difference from the previous TRGC work were unclear (below). \n2. There was little discussion about the reason why TCDF outperformed the proposed methods on the simulated fMRI dataset. \n3. There were no results in cost function (e.g., prediction error) even in the appendix. This detecting Granger causality would be unsupervised learning, thus the learning results based on the cost function are unknown. The information might help us understand the reason for 2. \n\nOther comments:\n\nIntroduction: \n“varying causal structures (Lowe et al., 2020)”: Did this paper focus on the varying causal structures?\n\nMethod:\nThe alpha in Alg.1 (threshold) is confusing with an elastic net regularization parameter.  \n\nThere is no explanation of Q of Alg. 1 in the main text. Is this a variable threshold in the analogous way of the ROC curve? How was the Q determined?\nIs the Alg. 1 the proposed method? The title “Stability-based thresholding” may be one of the components of the proposed method.\n\nExperiments:\n\nThe number of sequences to evaluate the methods in each dataset was unknown. In my understanding, the proposed methods learn each model for each sequence like the methods of Tank et al. (2018) and Khanna & Tan (2020). In the shared code, for example, Lorenz 96 used 5 sequences (in the code, called datasets). This may be critical information to evaluate the method.", "extra_scores": {}, "extra_text": {}}, {"review_id": "GL25BZf1vqH", "reviewer": "AnonReviewer1", "summary": "", "questions": "", "limitations": "", "rating": 4, "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 primarily deals with learning Granger-causal relationships in multivariate time series in the nonlinear dynamics setting. The core method uses vector autoregressive modeling with sparsity inducing regularizers (elastic net and smoothness based fused lasso) along with the recently proposed with self-explaining neural networks (for interpretability). The authors also augment the framework by learning Granger-causal structures that are stable on original and time-reversed data. Exhaustive empirical analysis is done with recent GC baselines.  Some of my concerns with the paper are the following\n\na) Paper does seem a bit incremental in terms of using elastic-net and fused lasso based smoothness regularization within VAR framework. The idea of bringing SENN within this is useful but not enough to claim significant novelty.\n\nb) Paper does not have scope of developing any useful theory ? Authors cite stability selection related work which has some decent guarantees on false discovery rate control. Can such theory be developed for the method proposed here ?\n\nc) On the positive side the empirical analysis and coverage of related work is very exhaustive. Of particular significance is the result presented in Table 3. However, given the double computation and higher number of parameters involved in the method, I was expecting some discussion on time complexity and how to improve scalability of the method in practice.", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Meta-review: The proposed approach is interesting and is differentiated enough from the recent body of work on Neural Network Granger causal modeling as it offers a mechanism for detecting signs of causality. 

The authors have satisfactorily addressed the points raised in the reviews. In particular relationship with prior work and novelty of the contributions are now clearly articulated. The added discussion on the superiority of TCDF on simulated fMRI experiments is insightful. Though prediction error is only a proxy for the task at hand, the readers will appreciate the added evaluation. 

The proposed approach to stability evaluation leveraging the time-reversal trick is novel and particularly pertinent, and could motivate some interesting follow-up work on this topic. It is also important that the authors have characterized the computational advantage of the approach. 
Decision: Accept (Poster)