Directed Cyclic Graph for Causal Discovery from Multivariate Functional Data
Authors: Saptarshi Roy
Raymond K. W. Wong
Yang Ni
Venue: NeurIPS 2023 poster (2023)
Keywords: Causal Embedding
Causal Discovery
Multivariate Functional Data
Directed Cyclic Graph
Causal Structure Learning
Bayesian Inference

Abstract: Discovering causal relationship using multivariate functional data has received a significant amount of attention very recently. In this article, we introduce a functional linear structural equation model for causal structure learning when the underlying graph involving the multivariate functions may have cycles. To enhance interpretability, our model involves a low-dimensional causal embedded space such that all the relevant causal information in the multivariate functional data is preserved in this lower-dimensional subspace. We prove that the proposed model is causally identifiable under standard assumptions that are often made in the causal discovery literature. To carry out inference of our model, we develop a fully Bayesian framework with suitable prior specifications and uncertainty quantification through posterior summaries. We illustrate the superior performance of our method over existing methods in terms of causal graph estimation through extensive simulation studies. We also demonstrate the proposed method using a brain EEG dataset.

Directed Cyclic Graph for Causal Discovery from Multivariate Functional Data
Saptarshi Roy
Raymond K. W. Wong Department of Statistics Texas A&M University College Station, TX 77843 raywong@tamu.ce edu
Yang Ni Department of Statistics Texas A&M University College Station, TX 77843 yni@stat. tamu.
Department of Statistics Texas A&M University
College Station, TX 77843 roys8001 1stttt. tamu. edu
Abstract
Discovering causal relationship using multivariate functional data has received a significant amount of attention very recently. In this article, we introduce a functional linear structural equation model for causal structure learning when the underlying graph involving the multivariate functions may have cycles. To enhance interpretability, our model involves low-dimensional causal embedded space such that all the relevant causal information in the multivariate functional data is preserved in this lower-dimensional subspace. We prove that the proposed model is causally identifiable under standard assumptions that are often made in the causal discovery literature. To carry out inference of our model, we develop fully Bayesian framework with suitable prior specifications and uncertainty quantification through posterior summaries. We illustrate the superior performance of our method over existing methods in terms of causal graph estimation through extensive simulation studies. We also demonstrate the proposed method using a brain EEG dataset.
Introduction
Motivation.
Multivariate functional data arise in many fields such as biomedical research Wei and
Li 2008 Chiou and Müller 2016), environmental science (Korte-Stapff al., 2022 finance Kowal
c al. 2017]. plant science Wong et al., 2019, Park et al.) 2022]. and sport science Volkmann et al.
202 where multiple variables are measured over time other domains. The increasing availability
of functional data in these fields provides us with great opportunities to discover causal relationships among random functions for the better understanding of complex systems, which helpful for various
machine learning and statistics tasks such as representation learning (Schölkopf et al., (202T) fairness
[Tang et al. 2023). transfer learning (Rojas-Carulla et al., 2018), and reinforcement learning (Zeng et 2023] One motivating example is electreencephalogrphhy (EEG) where electrical activity from the brain recorded non-invasively from electrode channels by placing them on the scalp or directly on the surface of the brain. Given its continuous nature and the short time separation between the adjacent measuring points, is natural to treat the data at each brain location/region as a function over time. A relevant scientific goal is to estimate brain effective connectivity among different regions, which will potentially allow us to make better decisions, design more effective interventions, and avoid unintended consequences. However, existing structural equation model (SEM) based causal discovery methods assume acyclic relationships among the random functions by imposing a directed acyclic graph (DAG) structure, which may be too restrictive for many real applications. For example, there are strong indications that brain effective connectivity studies, due to reciprocal polysynaptic connections, the brain regions far from exhibiting acyclicity [Friston] [011] Markov et al.) 20122. and that in genetic pathways, due to the presence multiple upstream regulators and downstream targets for every signaling component, feedback loops/directed cycles regular motifs [Brandman
37th Conference on Neural Information Processing Systems (NcurIPS 2023).

---

and Meyer, 2008]. Thus, in light of the prevalence of cycles in complex systems, it is desirable to have a flexible model for causal discovery among random functions that can account for such eyclic causal structures.
Challenges. Causal discovery for multivariate functional data in the presence of cycles is an inherently difficult problem th is not yet well understood. We highlight three prominent challenges. (i) Functional data are infinite-dimensional in nature. It may so happen that the low-frequency spectrum of one curve might causally influence the high-frequency spectrum of another curve. This demands identification of pertinent features that can be used to create a finite-dimensional representation of the data, which casier to work with and analyze. However, the challenge that we may not know a priori what the relevant features are when dealing with infinite-dimensionn objects. Blind adoption of standard -ao--aaualladaativve) low-dimensional features can lead to errors or inaccuracies. (ii) Although the identifiability of causal models for multivariate functional data in the absence of cycles has been established in recent works [Zhou et al., 2022b Lee and Li [022]] showing identifiability of causal models from multivariate data, let alone multivariate functions, is still a challenging and complex task in cases where causal relationships are obscured by the presence of cycles. (iii) It is common that functional data are only observed over discrete time points with additional noises. Such incomplete and noisy observations of the functions add another layer of difficulty in probing the causal relationships of interest.
Related work. Causal discovery from multivariate functional data has been studied by a few recent works [Zhou et al. 2022b Lee and Li, 2022 Yang and Suzuki, 2022]. which have already shown some promising results in discovering causality in, c.g., EEG data and longitudinal medical record data. However, all of them are limited to DAGs, which do not allow inference of cyclic causality. While there has been a surge of research on causal discovery methods for scalar random variables in the presence of feedback loops/eycles over the last few decades [Richardson, 1996 Lacerda et al.
2008 Mooij et al., 2011 Hyttinen et al. 2012 Huang et al., [ and Heskes 2013] Mooij and Claassen 2020 Zhou et al. 2022a], none these approaches have been extended to discovering causal dependencies among random functions in multivariate settings. Therefore, how to handle cyclic causal relationships among multivariate functional data while addressing the aforementioned challenges remains largely unsolved problem.
Contributions. In this paper, we propose an operator-based non-recursive linear structural equation based novel causal discovery framework that identifies causal relationships among functional objects in the presence of cycles and additional measurement/sampling noises. Our major contribution is four-fold.
1. We consider causal embedding of the functional nodes into a lower-dimensional space for dimension reduction that adapts to causal relationships.
2. We prove that the causal graph of the proposed model is uniquely identifiable under standard causal assumptions.
3. We capture fithin-function dependencies using data-driven selection of orthonormal basis that is both interpretable and computationally efficient.
4. To perform inference and uncertainty quantification from finite-sample data, we adopt a fully Bayesian hierarchical formulation with carefully selected prior distributions. Posterior inference is performed using Markov chain Monte Carlo (MCMC). We demonstrate the effectiveness of the proposed method in identifying causal structure and key parameters through simulation studies and apply the framework to the analysis of brain EEG data, illustrating its real-world applicability. Code will be made available on the project's website on Github
2 Model definition and causal identifiability
2.1 Notations
Let [p] = (1,....,p) for any positive integer p. causal directed cyclic graph (DCG) a graph
G = (V,g), which consists of set of vertices nodes V = [p] representing a set of random variables

---

and set of directed edges & = fe + j-i V} representing the direct causal relationships among the random variables. In a DCG, we do not assume the graph to be acyclic. A causal DCG model is an ordered pair (G, P)) where IP is a joint probability distribution over V (more rigorously, the random variables that V represents) that satisfies conditional independence relationships encoded by the causal DCG g. A simple directed cycle is sequence of distinct vertices (V1...., ukk such that the induced subgraph by these vertices is d + ... 7 2k + dj. For a vertex j V, we use pa(j) to denote the set of parents (direct causes).
2.2 Model framework
Consider multivariate stochastic process Y (I,,...,2, YT) where each Y, is defined on compact domain Tj C R. Without loss of generality, we assume 'Ti =... = Tp [0, 1). Suppose Yj E Hi where Hj is Hilbert space functions defined on Tj- We let (::)) denote the inner product of HHj. We propose a causal model that captures the relationships among ,y....,,;
Our proposed model considers an operator-based non-recursive linear structural equation model on the random functions Y as
Y2(-) = Σ (B,)))((++ Vj E [p].
(1)
(EPa(S)
where Bi linear operator that maps Hy Hj, and /g HH is an exogenous stochastic process. Clearly, for any j, E E V such that the edge f + j E & Bye is not null operator. Now by stacking the p equations in ( we obtain
BY + f,
(2)
where B = (Bye),,e= is a matrix of operators and f = fp) is a p-variate stochastic process. DAGs, the causal effect matrix can be arranged into lower block triangular structure given topological/causal ordering. But since our model allows for cycles, we have no such restriction on the structure of the operator matrix except that Bz5, Vj E [p]. is null, i.e., no self-loops.
Model I is infinite-dinee dimensional and hence challenging to estimate and interpret. To alleviate such difficulties. we consider oiw-imensioonl causal embedding structure. Specifically, we assume that the causal relationships are preserved in unknown low-dimensional subspace DD of Hy. Denote the dimension of Dj by Kj. Let Py and Q; be the projection onto D, and orthogonal complement in H, respectively. We assume Bye PPBrePP which implies that causal effects can be fully described within the oow-dimensional subspaces (D,} As such, can be split into
PjY, = Σ Byr(PPY))
(3)
(εP)()
2, = 2jlr-
We assume that Pjf, and f, are independent of each other. Now, by defining 0j P,,, and E; = Pyl,. [p]. (3) can be compactly written as
Bax + +,
where (aI,..... .T) and (p)) TT) with aj, E, E D,, Vj E [p].
In practice, the random functions in Y can only be observed over finite number of (input) locations, possibly with measurement errors. More specifically, for each random function Y,. we observe [(tju:Xjin { where Xju E RR is the measurement of Y, at location tju Ty and m, is the number of measurements obtained from Yy. Defining B, 2,Y,. we consider the following measurement model:
Xju = Y,(tuu) + teu aa(tju)) +B((tju) Cjui Yu E [myli [p].
(5)
with independent noises Eju N(0, αy), Vu E [my]-
More compactly, (5) can be written
X = x=) + +(t(t) +e,
(6)
3

---

where X
....pppp)) a(t) (ai(t1) ..... ap(tp) ((t)
=
(Bi(ti) ..... ,, T) and (e ..... eT with X, (Xji, "ij....;jjm)) (a,(()))...., az ()(()mn))) (3,(,) (8,(4)))... 3,(n))...,,,,)m)m)). 3,(1mm))) and ej (eg1....., Ejmy)
We call our proposed model, FENCE, which stands for FFunctional Embedded Nodes for Cyclic causal Exploration', reflecting its purpose.
2.3 Causal identifiabilty
In this section, we shall show that the graph structure of the proposed FENCE model is identifiable for functional data measured discretely with random noises under several causal assumptions. We start by defining causal identifiability and state assumptions.
Definition 2.1. (Causal Identifiability) Suppose Y p-variate random function and is the observed noisy version of Y given by (6). Assume X follows FENCE model S = (G, IP) where G is the underlying graph and P the joint distribution of X over G. We say that causally identifiable from X if there does not exist any other S (g", IP*) with g* * g such that the joint distribution 1 on X induced by G is equivalent to induced by g.
In other words, for causal graph to be identifiable, there must not exist any other graph such that the joint distributions induced by the two different graphs equivalent. Next, we list and discuss few assumptions to establish the causal identifiability of the proposed model.
Assumption 1. (Causal Sufficiency) The model S = PG is causally sufficient, i.e. there are no unmeasured confounders.
Assuming no unmeasured confounders keeps the causal discovery task more manageable especially for cyclic graphs with purely observational data.
Assumption 2. (Disjoint Cycles) The cycles in are disjoint, i.e., no two cycles in the graph have two nodes that are common to both.
Assuming disjoint cycles induces a natural topological ordering and forms a directed acyclic hypergraph-like structure within the DCG. The same assumption was made in Lacerda et al. [2008)
Assumption 3. (Stability) For the model S, the moduli of the eigenvalues of the finite rank operator are less than or equal to 1, and none of the real eigenvalues are equal to 1.
According to Fisher 1970 the SEM in can be viewed as being in state equilibrium, where the finite rank operator B represents coefficients in set of dynamical equations that describe a deterministic dynamical system observed over small time intervals as the time lag approaches zero. The eigenvalue conditions are deemed necessary and sufficient for the limiting behavior to hold, as argued by Fisher 1970 Such an assumption widely adopted in e.g., econometrics, and Lacerda et al. 2008 made this assumption well.
Assumption 4. Won--aussianity The exogenous variables have independent mixture of Gaussian distributions. i.e., Ejk
ind Myr TylkeN(Pjkmm Tykm) with Myx 2 2.
The assumption of non- Non-Gaussianity on the exogenous variables has been proven useful in causal discovery as it induces model identifiability in the linear SEM framework (Aapo and Petteri 1999 Shimizu al. 2006 Lacerda al. 2008 Spirtes and Zhang 2016). Mixture of Gaussian can approximate any continuous distribution arbitrarily well given a sufficiently large number of mixture components (Titterington et al. 1985 McLachlan and Peel 2000, Rossi 2014). It is also easy to sample, which facilitates our posterior inference.
C(t),, where C(t) diag( (CuuCt),......... Cpp(ty) and 7 represent another exogenous component of the model. Here C3)(t)) is mixing matrix that mixes the independent entires in y to generate temporal dependence
Assumption 5. (Non-causal dependency) We assume B(t)
ind JM Jkkm T Mi' (km) with Mss 2 1.
within the j-th block. We assume Yikk
Since the model assumes that all causal information in Y is preserved in the lower-dimensional space DD and not in its orthogonal complement, it is apparent that while each (3,(t)) within a block can have temporal dependence, it is independent of (tt)) when # and j. E E [p].

---

For some basis {d)x) Kg that spans the low-dimensional causal embedded space Dj. og in can be further expanded by, (i((iu)) = E (t Therefore 6 can be written more compactly as
5
+t( + ((t) +e,
where (t) = diag( 11 (ti)...... Pp(tp)) with + ,,(tt) = (tjo(tsu)) m,KK,
Assumption 6. (Sufficient sampling locations) The basis matrix + (t) of size Σi=1 my Σ Kj has full column rank.
This assumption implies enough sampling locations, over which each random function Yj is observed, to capture all the causal information that Y contains.
Given these six assumptions, main theorem establishes the causal identifiability of the proposed model.
Theorem 2.1. Under Assumptions 6 S PP) is causally identifiable.
The proof essentially involves two steps as shown in Figure On the left-hand side (LHS) of the diagram, we depict the hypergraph-like structure that emerges when assuming the existence of disjoint cycles (Assumption 2 whereas, on the right-hand side (RHS), we offer magnified view of the hypernodes (nodes containing simple directed cycle). Our approach to proving causal identifiability progresses from the LHS to the RHS. That is, we first prove the identifiability structure depicted on the LHS of Figure [ and then we proceed to establish the identifiability of each simple directed cycle within every hypernode in the hypergraph. The detailed exposition of the proof can be found in Section A of the Supplementary Materials.
the hypergraph-
Q
Figure Two important components of causal identifiability proof: (I) identifiability of directed acyclic hypergraph induced by disjoint cycles, and (II) identifiability of each disjoint cycle.
Bayesian model formulation
In this section, will describe the inference procedure of the proposed model. A straightforward approach would be a two-step procedure where the first step performs functional principal component analysis on each function marginally reduce the dimension, and then the second step learns causal structure based on the principal components. However, this simple approach has several disadvantages. First, the estimated functional principal components that explain the most variation of each individual function marginally may not optimally capture the cause-effect dependence relationships among different functions. Second, this procedure is unreliable since estimation uncertainty fails to propagate correctly from the first step to the second step. As such, we propose a fully Bayesian approach, which reduces the dimension of functional data adaptively for causal structure learning.
3.1 Model parameters
Let E (E,0,,== denote the adjacency matrix where Eze = indicates the existence of a directed edge from node & to node j. and 0 otherwise. Let {dv 1 be set of S common unknown basis
5

---

functions that approximate each random function Y, i... = denote the set of basis coefficients. Note that {ór} not the basis for the lower-dimensional causal embedded subspace Dj. However, we assume that the first Kj of them actually spans Dj and our goal is to hunt for them through a properly designed inference procedure. Moreover, according to our assumption, we build our SEM on the first K of the basis coefficients á = (ăjI,"" 1 "ar""""nkk) ,ūkK,)
āzāăz where {āzz 1
ki=
Defining a; (a,, aas)
with a; =Yj =-jjoyy they can be written as
à == Bā - E, where (āT. ... a a aT
(8)
(ET E 7 ..... Yp with Ej
)T
B
(E)2, {EjKK) and Y = (Ti,Kt1I..... ,,s)) Here B =
where
0
(BB ())-1),)))) with Bjj- 0 since we assume the absence of self loops. To carry out inference, we assume Ezks Yikk ~
ind D jkkm N(ukkm Tykm).
3.2 Adaptive basis expansion
As the dh''s 's are specifically useful for restricting the original function space for each Y, to a lower- dimensional causally embedded smooth space of dimension Kj, we make the basis { adaptive for causal structure learning by further expanding them with known spline basis functions (Kowal al. 2017]. (hx(t) Σ Akrbr(-), where b (b1,.... :....,,b)) is the set of fixed cubic B-spline basis functions with equally spaced knots and A (ARI,..., ARR) T are the corresponding spline coefficients. Since we do not fix A&'s priori, the basis functions de 's can be learned from data a posteriori and hence are adaptive to both data and causal structure (i.e., the basis functions, the functional data, and the causal graph are dependent in their joint distribution).
3.3 Prior specifications
Prior spline coefficients. The prior Ag is chosen serve multiple purposes. (i) It sorts the basis functions by decreasing smoothness and therefore helps to identify the spanning set size K2 for the underlying smooth causally embedded space Dj- (ii) Although a strict requirement for modelling purpose, it forces Φą 's to be orthonormal, i.e. f dx (w)( dw = I(k = k)). As such, the orthogonality constraints help eliminate any information overlap between the basis functions, which keeps the total number necessary basis functions that actually contribute to the causal structure learning to a minimum. (iii) It regularizes the roughness of Φu 's to prevent overfitting.
For (iii), more specifically, we restrict the roughness of the basis functions (-) by assigning a prior that penalizes its second derivatives Gu 1992 Wahba 1978 Berry al. 2002
A (N A 22)),
where 5 is the pseudoinverse of Ω = f [ (t) 6 (t) dt. Let Ω U be the singular value decomposition of Q. Following Wand and Ormerod 2010 to facilitate computation,
reparameterize (bx(-) Σ A-) with b(.) (1,4, 6² ()ÜD - where D is the (R 2) x (R 2) submatrix of D corresponding to non-zero singular values (note that the rank of Ω is R - 2 by definition) and Ü is the corresponding R x (R 2) submatrix of U. This induces prior on Å given by
r=1
A NO0, Sx) with S₁ diag( 00, A , A
In other words, the intercept and the linear term unpenalized but the non-linear terms are penalized, the degree of which is controlled by Ag- In practice we set the first two diagonal elements of S₁ as 10 We constrain the regularization parameters A1 >
As > 0 by putting uniform prior:
A ~ Uniform Lx. UR). V [8].
U, 100 La A V k E [S 1], U² A 4 k f2, ... Lt 10
which implies that the smoothness of (x))) decreases as k gets larger.

---

Priors on the adjacency matrix. We propose to use an independent Beroolllli prior on each entry Eye of E, ie.. Eşelp ind ~ Bernoulli(p) and p ~ Uniform(0, 1). of E with p integrated out is given by
The marginal distribution
E / P(ElP))(O) dp = Beta
E j7t
Now, for example, Eg denotes the null adjacency matrix and E₁ denotes the adjacency matrix with only edge, then over graph with one edge by a factor of p² p. and, importantl, this penalty increases with p. Thus, the
see that p( /p)
p. Therefore, an empty graph is favored
prior prevents false discoveries and leads to sparse network by increasing the penalty against additional edges as the dimension grows.
Prior on the causal effect matrix. Now given E, we assume an independent spike and slab prior on the entries of (B,0),e=
E ~ (1 -1- (.-Eee)MVV(N 0 s, Ix,) ( + (B, Ix,), where MVN 0, 1I,, Ik.) is a matrix- normal distribution with row and column co- variance matrices k, and Ik, respectively. We assume a conjugate inverse-gant -arse- prior on the causal effect size, InverseGamma(u,, by). We choose dyy by 1. We fix 0.02 so that when E = 0, Bjr negligibly small.
Priors on the parameters of the Gaussian mixture distribution. the parameters of the Gaussian mixture distribution:
We choose conjugate priors for
( (Tkk,...., TzkM,t) Dirichlet(3...... B), Vj E [p]. k E [S]
Pijkm N Na bn), Tjkm ~ InverseGamma(uu by), V E (p), E [S]. E [M,a] We have fixed values for the hyperparameters, = 1, 0g 0, by = 100, a, = b, = 1.
Prior on the noise variances. We assume a conjugate prior for oj InverseGamma(o. bg), j [p]. We choose ag = b, 0.01.
We simulate posterior samples through Markov chain Monte Carlo (MCMC). Details are given in Section B of the Supplementary Materials. Sensitivity analyses will be conducted test the hyperparameters including (a., (a,, by), (ao, ba), R, S, M and B.
Simulation study
Data generation The data were simulated according to various combinations of sample size (n), number of nodes (p), and grid size (m, = d jj E (p)) where n E (75, 150, 300], p E (20, 40, 60). and d E {125, 250). The grid evenly spans the unit interval [0, 1]: the results with unevenly spaced grids are presented in Section of the Supplementary Materials. The true causal graph G was generated randomly with edge formation probability 2/p. Given G, each non-zero block By of causal effect matrix was generated from the standard matrix-variate normal distribution. We set the true number of basis functions to be K = 4. In order to generate K = 4 orthonormal basis functions, we first simulated unnormalized basis functions by expanding them further with 6 cubic B B-spline basis functions where the coefficients were drawn from the standard normal distribution and then empirically orthonormalized them. The basis coefficients were generated following (8) with the exogenous variables E, drawn independently from Laplace distribution with location parameter = and scale parameter = 0.2. We have also considered other -an-Gussiin distributions for the exogenous variables; the corresponding results provided in Section Of the Supplementary Materials. Finally, noisy observations were simulated following (6) with the signal-to- noise ratio, i... the mean value of ()) sample, where [n)].
(t)l/a, across all and t, set to 5. Here, superscript (t) denotes the ith
For the implementation of the proposed FENCE, we fixed the number of mixture components to be 10 and ran MCMC for 5,000 iterations (discarding the first 2,000 iterations burn-in and retaining every 5th iteration after burn-in). The causal graph G was then estimated by using the median probability model Barbieri and Berger 2004), i.c.. by thresholding the posterior probability of inclusion at 0.5.
7

---

Methods for comparison. We compared our method with fLiNG (Zhou et al. 2022b] recently proposed directed acyclic graph (DAG) for multivariate functional data. Codes for Lee and Li 2022].
Yang and Suzuki [022] are not publicly available. Hence for more comparison, we considered two ad hoc two-step approaches. the first step of both approaches, we obtained the basis coefficients by carrying out functional principal component analysis (fPCA) using the fdapace (Zhou et al. [022c]
package in R. Then in the second step, given the basis coefficients, we estimated causal graphs using existing causal discovery methods, (i) LiNGAM [Shimizu et al. 2006) (ii) PC (Spirtes and Glymour 1991] and (iii) CCD [Richardson] 1996]: we call these three approaches IPCA-LiNGAN, IPCA-PC and fPCA-CCD respectively. Note that we did not use SEM-based eyclie discovery algorithm, LiNG-D [acerda et al. 2008) in the second step due to the unavailability of the code.
LiNGAM estimates a causal DAG based on the linear non non- Gaussian assumption whereas PC generally returns only an equivalence class of DAGs based on conditional independence tests. CCD algorithm is constraint- causal discovery method, which yields an equívalence class of cyclic causal graphs. LiNGAM and PC are implemented in the pcalg package [Kalisch et al., 2018 in R. CCD algorithm is implemented in the py-tetrad [Ramsey et al. 2018 package in python.
Performance metrics. To assess the graph recovery performance, we calculated the true positive rate (TPR), false discovery rate (FDR), and Matthew's correlation coefficient (MCC). For TPR and MCC, higher is better, whereas lower FDR is hetter.
Table I summarizes the results of 50 repeat simulations, demonstrating that the proposed FENCE model outperforms all competitors (fLiNG, fPCA-LiNGAM, and (PCA-CCD) across all combinations of n, p, and We provide the results of 1PCA-PC in Section Df the Supplementary Materials, which are similar to those of IPCA-LINGAME counting a non-invariant edge between two nodes as a true positive as long as the two nodes are adjacent in the true graph. The superiority of FENCE is not unexpected for three reasons. First, fLiNG, IPCA-LiNGAM, and IPCA-PC are not specifically designed for learning cyclic graphs. Second, two-step approaches like IPCA-LiNGAM, fPCA-PCC, and fPCA-CCD do not necessarily capture the causally embedded space through the functional principal components. Third, although fPCA- can handle cyclic graphs, it, being two-step approach, fails to capture true functional dependencies. Overall, these findings provide strong evidence of the effectiveness of FENCE compared to existing methods.
Results.
We favored iPCA-PC and 1PCA-CCD by
Additional simulations. We considered additional simulation scenarios with unevenly spaced grids, general exogenous variable distributions, true acyelic graphs and data generated using non- non-linear SEM, and also conducted sensitivity analyses of FENCE with respect to several hyperparameters; the results are presented in Section C of the Supplementary Materials. The performance of the proposed method is consistently better than competing methods and is relatively robust with respect to hyperparameter choices.
Table 1: Comparison of performance of various methods under 50 replicates. Since LiNGAM is not applicable to cases where q > n with q = Kp being the total number of extracted basis coefficients across all functions, the results from those cases are not available and indicated by "-" The metrics reported are based on 50 repetitions are reported; standard deviations are given within the parentheses.
FENCE
fLING
IPCA-LINGAM
IPCA-OCD
TPR
FDR
MCC
TPR
FDR
TPR
MCC
FDR
0.8501391
8.190.87)
7.7000155
0.350.19
0.100044
02500057
75
40
0.7H0000)
6.230.80)
0.86/8.94
0.37:0085
0.82/0.6)
0.35/005)
0.730002)
0.4700045
021:0.05)
75
0.75
0.27/0.05
0.83/8.04
0.34/OU07)
0.330005)
0.32/004)
0.6800.00)
0.610.05)
0.1900.3)
1:0
0.8801.071
6.1400.6)
6.998.05)
0.7500.05
0.39/005)
0.86/0.05)
0.0800.05)
0.710.03)
0.400.03)
0.200.44
150
0.8100.77
0.2100.6)
0.8710.55
0.750005)
0.37/0.04
0.35/0223
0.91/0..2)
0.0800.06)
0.730.04)
0.710005)
02250.05
1.0
60
0.79)0.6))
6.34/0.85)
0.866.04
3660004)
0500005)
0366005
0.7200.55
0 (4(00.0))
022/0003
300
20
125
0.91M188)
6.0998.04)
0.998.04
0.75(0055)
0.41:0.44)
0.30/0.19)
0.5410.25)
0.1100.29)
0.100.03)
0.390.04)
0.260.033
300
40
125
0.5770.94
6.15/0.85)
8.8708U5)
0.75(0055)
0.38/0.55
0.27/0.20
0.91/0.22)
0.08(0.06)
0.7710.03
04500.23
0.240.03)
100
60
12.5
(.85)0885
6.17H0.83)
0 (7760004)
D 03360005
0.78/0-17)
(..0
0.0510.00)
7.720.00
0 44000000
(222005)
75
20
250
0.8100.44
6.2348.02)
8.85
0L3RKOUTT)
0.000005)
0.39.0.4)
0.32/0.14)
0.52/0.03)
0.0900.04)
0.6770003)
0.4000033
0.220.04)
25
40
0.73/066)
6.2510.86)
0880005
0339005
0.600.00
051000
0210005
75
0.6770.03)
6.3460.8))
0.798.04
(4400044
0880055)
0.112004)
0.610.04
5.50000
(99004))
190
0.5300.06
6.17/0.85)
0.8668.05
0.73/007)
042/0.05)
0.32/0
0.7900055
0.1310.05)
0.730:04
0.360.03)
02400.3)
0710005)
(.81(0.00)
0.1300.00)
0.990.85)
0.490005
0.9900.00)
0216002)
300
10
0.86/.92)
6,1610.84)
6.85-8.04
0.72000
04472004)
0-450
0,3900.55
0. 1.100000
2960.00)
(.440000)
0,20003)
100
40
250
0.7910.6)
6.1610.85)
(.84+8.06
(740055
07710055
0.430005
0.3900
0.82(0.5)
1640.07)
2640.055
0 (.4000.0)
0:250005)
100
40
250
0.76/0.05)
6.21(0.83)
(.88-8.0)
(L7708.5)
0440005)
0.420055
0.2800.17 (7)
0.5000.04)
0. (.00044
0.7240.06)
0.53,0.03)
(222/004)
8

---

Real data application
Brain EEG data. We demonstrate the proposed FENCE model brain EEG dataset from an alcoholism study [Zhang et al. 1995]. This dataset was earlier used to demonstrate functional undirected graphical models [Zhu et al., 2016 Qiao et 2019] and functional Bayesian network [Zhou et al., [022b]] Data were initially obtained from 64 electrodes placed on subjects' scalps, which captured EEG signals at 256 Hz (3.9 ms epoch) during a one-second period. The study consists of 122 subjects, out of which 77 are in the alcoholic group and 45 are in the control group. Each subject completed 120 trials. During each trial, the subject was exposed to either single stimulus (a single picture) or two stimuli (a pair of pictures) shown on computer monitor. We particularly focus on the EEG signals filtered at a frequency bands between 8 and 12.5Hz using the eegf ilt function of the leglab toolbox of Matlab as band signals are associated with inhibitory control Knyazev 2007]. Given that the EEG measurements were recorded from each subject over multiple trials, these measurements are not independent of each other due to the time dependency of the trials. Moreover, since the measurements were obtained under various stimuli, the signals may have been affected by different stimulus effects. To mitigate these issues, band-filtered EEG signals for each subject across all trials under a single stimulus, resulting in a single event-related potential curve per electrode per subject. By doing so, eliminated the potential dependence between the measurements and the influence of different stimulus types. We performed separate analyses of the two groups to identify both the similarities and dissimilarities in their brain effective connectivity.
calculated the average of the
We conducted a Shapiro- Shaliro-Wilk Wilk normality test on the observed functions for each of the p 64 scalp positions each of the my 256 vj E [p] time points evaluate their Gaussianity. The results showed that for numerous combinations of scalp position and time point, the null hypothesis (which assumes that the observations are marginally Gaussian) was rejected. Thus, we conclude that the non-Gaussianity of the proposed model is appropriate. Next, for posterior inference, we MCMC for 20,000 iterations, discarded the first half as burn-in, and retained every 10th iteration after burn-in. The estimated causal graph by thresholding the posterior inclusion probability to 0.9 is given below in Figure 2
Lde
Lode
Figure 2: Estimated causal brain connectivity from EEG records by FENCE with posterior probability of inclusion 0.9, separately for the alcoholic (left) and control (right) group. The bi-directed edges are just directed cycles, i.c., ++ j means i + j and i +
Results. There are some interesting findings. First, for both groups (alcoholic and control), brain regions that are located in adjacent positions tend to be more connected than the brain regions that are far apart. Second, dense connectivity is observed in the frontal region the brain in both groups, with multiple cycles being formed. Third, compared to the control group, the alcoholic group has more connectivity across the left parietal and occipital lobes. Fourth, the same cycle of Iz, Cz, and RPA is observed in both groups

---

Validity. We now discuss the validity of our real data results. In Hayden et al. 2007 it was observed that alcohol-dependent subjects exhibited frontal asymmetry. distinguishing them from the control group. Our own investigation aligns well with these results, as we have identified denser connectivity across various brain regions in the middle and left areas the frontal lobe among alcoholic subjects, when compared to controls. Furthermore, Winterer et al.) [2033 documented coherent differences between alcoholics and controls in the posterior hemispheres, specifically in the temporal, parietal, and occipital lobes. In accordance with their findings, study provides additional support for this claim, as we have observed heightened activity with several cycles formed in those same regions within the alcoholic group when compared the control group.
Discussion
We briefly highlight here several potential avenues for future development of our current work. First, an intriguing and important direction would be to explore the relaxation of the causal sufficiency assumption in the model identifiability. Second, our current model is based on linear non- non-Gaussian assumption over the exogenous variables, but a nonlinear model could be considered as an alternative. Lastly, an alternative approach to determining the effective number of basis functions that span the liwer-dimensional causal embedded space would be to utilize the ordered shrinkage priors [Bhattacharya and Dunson 2011 Legramanti et al., 2020 order to adaptively eliminate redundant components, resulting in more flexible methodology.
Acknowledgements
Ni's research was partially supported by NSF DMS-211222 129443 and NIH IR0IGM148974-001.
References
Hyvärinen Aapo and Pajunen Petteri. Nonlinear independent component analysis: Existence and uniqueness results. Neural Networks, 12:429 439, 1999. URL hww.ddoo.oogg//0/ 50893-6080 (98)00000-0
Maria Maddalena Barbieri and James O. Berger. Optimal predictive model selection. The Annals of Statistics, 32(3);870 897, 2004. doi: 10. 00121400003334440000033 httpps:////doi org/ 10. 112114//0903336400000888
Scott M Berry, Raymond J Carroll, and David Ruppert. Bayesian smoothing and regression splines for measurement error problems. Journal of the American Statistical Association, 97 (457):160-166, 2002. doi: 10.11980066214502554777388 016214502753479801
URL httts:://oi..org//e/...... 1098////./9
A. Bhattacharya and D. Dunson. Sparse Bayesian infinite factor models. Biometrika, 98(2): 291-306, 06 1. ISSN 0006-3444. doi: 10. 10.1093//iomet/ast011.. URL https://doi... org/00 1093/biomet/asse asr013
Onn Brandman and Thomas Meyer. Feedback loops shape cellular signals in space and time. Science, 322(5900):390-333) 395, 2008.
Jeng-Min Chiou and Hans Hans-Georg Müller. A pairwise interaction model for multivariate functional and longitudinal data. Biometrika, 103(2):377-396 2016.
Franklin M. Fisher. A correspondence principle for simultaneous equation models. Econometrica, 38 (1):73-92, 1970. ISSN 00129682, 14680262. URL URL http:httt/////t jstor httttorr.rrr/stall////99999
Karl J. Friston. Functional and effective connectivity: A review. Brain Connectivity, 1(1):13-36. 201 doi: 10. 1089//0aain.200 1.0088. URL. httph::tt/////// 1r0//000.00////rra....
1009...989bbbra.. 2011.
PMID: 22432952.
Chong Gu. Penalized likelihood regression: A bayesian analysis. Statistica Sinica, 2:255-264, 1992.
10

---

Elizabeth Hayden, Ryan Wiegand. Eric Meyer, Lance Bauer, Sean O'Connor, John Nurnberger, David Chorlian, Bernice Porjesz, and Henri Begleiter. Pattems of regional brain activity alcohol- dependent subjects. Alcoholism, clinical and experimental research, 30:1986-91, 01 2007. doi: 10.1111/j.15000222 1530-0277.2005 X..4044.000..0
Biwei Huang, Kun Zhang. Jiji Zhang, Joseph D. Ramsey, Ruben Sanchez-Romero, Clark Gly- mour, and Bernhard Schölkopf. Causal discovery from helerngennousioossttinnyyy data. CoRR, ttss/1933.0672,, 2019. http://a/r/ http://wwwii.rr..raaaoo/. wwarg//bss//ie 1903 1903. 01672
Antti Hyttinen, Frederick Eberhardt, and Patrik O. Hoyer. Causal discovery of linear cyclic mod- els from multiple experimental data sets with overlapping variables. In Nando de Freitas and Kevin P. Murphy, editors, Proceedings of the Twenty-Eighth Conference on Uncertainty in Arti- ficial Intelligence Catalina Island, CA, USA, August 14-18, 2012, pages 387-390 AUAI Press. 2012. URL https:://sspppittee ottp/////a////apaaaaaaaaoooooooo
jsp?mmmu &smnu 2&
article_id-200116 2012p0occooooooenggg111
Markus Kalisch, A. Hauser, Marloes H. Maathuis, and Martin Mächler. An overview of the pealg package for r, 2018.
Gennady G Knyazev. Motivation, emotion, and their inhibitory control mirrored in brain oscillations. Neuroscience & Biobehavioral Reviews, 31(3:3777-395, 2007.
Moritz Stapff, Drew Yarger, Stilian Stoev, and Tailen Hsing. A multivariate functional-data mixture model for spatio-temporal data: inference and cokriging, 2022.
Daniel R Kowal, David S Matteson, and David Ruppert. A bayesian multivariate functional dynamic linear model. Journal of the American Statistical Association, 112(518) 2017.
Gustavo Lacerda, Peter L.. Spirtes, Joseph Ramsey, and Patrik O. Hoyer. Discovering cyclic causal models by independent components analysis. In Conference on Uncertainty in Artificial Intelli- gence, 2008
Kuang- Yao Lee and Lexin Li. Functional structural equation model. Journal of the Royal Statistical Society Series B, 84(2):600 629, April 2022. doi: 10.1111/rssb.... 12471. URL wttphhttt::/// //ideas repec. poppe..rg//oe a/ http:///orrrbb///y2222222ppp
629 .html
Sirio Legramanti, Daniele Durante, and David B Dunson. Bayesian cumulative shrinkage for infinite factorizations. Biometrika, 107 05 2020. ISSN 0006-3444. doi: 10.1093/biomet/ asaa008. URL https::::///o 10930....1999/
N. T. Markov, M. M. Eresey-Ravasz, A. R. Ribeiro Gomes, C. Lamy, L. Magrou, J. Vezoli, P. Misery, A. Falchier, R. Quilodran, M. A. Gariel, J. Sallet, R. Gamanut, C. Huissoud, S. Clavagnier, P. Giroud, D. Sappey-Marinics P. Barone, C. Dehay, Z. Toroczkai, K. Knoblauch, D. C. Van Essen, and H. Kennedy. Weighted and Directed Interareal Connectivity Matrix for Macaque Cerebral Cortex. Cerebral Cortex, 24(1):17-36, 09 2012. ISSN 1047-3211. doi: 10.1093/cercor//bs270. URL battps://doo org 10. 1093/ cercor /bh2270
G. J. McLachlan and D. Peel. Finite mixture models. Wiley Series in Probability and Statisties, New York, 2000.
Joris M. Mooij and Tom Claassen. Constraint-based causal discovery using partial ancestral graphs in the presence of cycles. In Ryan P. Adams and Vibhav Gogate, editors, Proceedings of the
Thirty-Sixth Conference on Uncertainty in Artificial Intelligence, UAI 2020, virtual online, August 3-6, 2020, volume 124 of Proceedings of Machine Learning Research, pages 1159-1168. AUAI Press, 2020. URL htttt::ppppoceeinnggeee .ml.. wwtpppeeee///aaeeeeeeeeee m00aa html
Joris M. Mooij and Tom Heskes. Cyclic causal discovery from continuous equilibrium data. In
Ann E. Nicholson and Padhraic Smyth, editors, Proceedings of the Twenty-Ninte Conference on Uncertainty in Artificial Intelligence, UAI 2013, Bellevue, WA, USA, August 11-15, 2013. AUAI Press, 2013. URL Utt htw.l//ow./aaaaat.o///////////aa fummunu

---

Joris M. Mooij. Dominik Janzing, Tom Heskes, and Bernhard Schölkopf. On causal dis- covery with cyclic additive noise models. In John Shawe-Taylor, Richard S. Zemel, Pe- ter L. Bartlett, Fernando C. N. Pereira, and Kilian Q Weinberger, editors, Advances in Neural Information Processing Systems 24: 25th Annual Conference on Neural Information Processing Systems 201 I. Proceedings of meeting held 12-14 December 2011, Granada, Spain, pages 639-6499-777 647, 2011. URL httpp:://poocceedigggg neurin. htappp.cc..poooo httpp/appr/21//hhaa///
Abstract .html
Yeonjoo Park, Bo Li, and Yehua Li. Crop yield prediction using bayesian spatially varying coefficient models with functional predictors. Journal of the American Statistical Association, pages 1-14, 2022.
Xinghao Qiao, Shaojun Guo, and Gareth M James. Functional graphical models. Journal of the American Statistical Association, 114(525):211-22 12222 2019.
Joseph Ramsey, Kun Zhang, Madelyn Glymour, Ruben Sanchez Romero, Biwei Huang, Imme Ebert-Uphohee Savini Samarasinghe, Elizabeth A Barnes, and Clark Glymour. Tetrad a toolbox for causal discovery. In 8th international workshop on climate informatics, 2018.
Thomas Richardson. A discovery algorithm for directed cyclic graphs. In Proceedings of the Twelfth International Conference Uncertainty in Artificial Intelligence, UAI'96, page 454 461, San Francisco, CA, USA, 1996. Morgan Kaufmann Publishers Inc. ISBN 155860412X.
Mateo Rojas-Carulla, Bernhard Schölkopf, Richard Turner, and Jonas Peters. Invariant models for causal transfer learning. J. Mach. Learn. Res., 19(1):309--134) jan 2018. ISSN 1532-4435.
Peter E. Rossi. Bayesian Non- and Semi- paraaneriamccccccc Methods and Applicatioons. Princeton University Press, 2014, ISBN 9780691 45327. htRlhhtp:://wwx. jstor wwwportagggge....ooooo...../ee.t
Bernhard Schölkopf, Francesco Locatello, Stefan Bauer, Nan Rosemary Ke, Nal Kalchbrenner, Anirudh Goyal, and Yoshua Bengio. Toward causal representation learning. Proceedings of the JEEE, 10965):612-634, 2021.
doi: 10. 10110///R..2222..335999999
Shohei Shimizu, Patrik O. Hoyer, Aapo Hyvrimen, and Antti Kerminen. A linear non-gaussian acyclic model for causal discovery, Journal of Machine Learning Research, 7(72):2003-2330, 2006. URL
http://jmlr http:///...ppppp///////////////a http///shiiiiim../ html
Peter Spirtes and Clark Glymour. An algorithm for fast recovery of sparse causal graphs. Social
Science Computer Review, 9(1)662-72, 1991. doi: 101177/899449331199900066 URL
//doi.org//i 11111.11777000444399000009000000
Peter Spirtes and Kun Zhang. Causal discovery and inference: concepts and recent methodological advances. Applied Informatics. 3. 2016. doi: 10.11866440335-01--0018844 URL https://doi
org/10. 100ee.116/ss/ 040534400550000000000
Zeyu Tang, Jiji Zhang, and Kun Zhang. What-is and how-to for faimess in machine learning: survey, reflection, and perspective. ACM Computing Surveys, 55(13s):1-37, jul 2023. doi:
10. N0.1145/359719 1999. URL
D.M. Titterington, A.F.M. Smith, and U.E. Makov. Statistical Analysis of Finite Mixture Distributions. Wiley, New York, 1985.
Alexander Volkmann, Almond Stöcker, Fabian Scheipl, and Sonja Greven. Multivariate functional additive mixed models, 2021.
Grace Wahba. Improper priors, spline smoothing and the problem of guarding against model
errors regression. Journal of the Royal Statistical Society: Series (Methodological), 40(3): 364-372, 1978. doi: wtwp////ec...///////////o...... 1tt.://00..11111....2666........ 1978. 1000000
1978 10010500.. U //rss 2ttp://.jj21171-111
onlinelibrary. www.eeeawbbrr..llle.ooo.ooo
M. P. Wand and J. T. Ormerod. Corrigendum: On semiparametrric regressson with o'sullivan penalised splines. Austratian & New Zealand Journal of Stulistics, 521)))2icc..552222229999999 2010. doi: https:/doi..og! 10.1111/j. 1467- 201 00578.x. URL htttpe::/omlinellieraa wiley. www.lyy...o.....ooooo.e /
10. 1111/11jjjjjp j. 1467 842X 2010. 00578.x
12

---

Zhi Wei and Hongzhe Li. A hidden spatial-temporal Markov random field model for network- analysis of time course gene expression data. The Annals of Applied Statistics, 2(1):408 429, 2008. doi: 10. 10.1214077-AAAS44.5
https::///// org/ 10. 1214/07 A0AS145
G Winterer, M.-A Enoch, K White, Mete Saylan, Richard Coppola, and D Goldman. Beg phenotype in alcoholism: Increased coherence in the depressive subtype. Acta psychiatrica Scandinavica, 108:51-60, 08 2003. doi: 10.1034/j. 16044/4j.600004444 00447.2003.00066.0
Raymond KW Wong, Yehua Li, and Zhengyuan Zhu. Partially linear functional additive models for multivariate functional data. Journal of the American Statistical Association, 114(525):406-118, 2019.
Tianle Yang and Joe Suzuki. The functional LiNGAM. In Proceedings of The 11th International Conference on Probabilistic Graphical Models, volume 186 of Proceedings of Machine Learning Research, pages 25-36. PMLR, 05-07 Oct 2022. URL https:://poceeddaiits
mlr mwr.m.r.rrr.... press/
v186/yang224 html
Yan Zeng, Ruichu Cai, Fuchun Sun, Libo Huang, and Zhifeng Hao. A survey causal reinforcement learning, 2023.
Xiao Lei Zhang, Henri Begleiter, Bernice Porjesz, Wenyu Wang, and Ann Litke. Event related potentials during object recognition tasks. Brain Research Bulletin, 38:531-5555 1995.
Fangting Zhou, Kejun He, and Yang Ni. Causal discovery with heterogeneous observational data. In James Cussens and Kun Zhang, editors, Uncertainty in Artificial Intelligence, Proceedings of the Entirty-Eighhte Conference on Uncertainty in Artificial Intelligence, UAI 2022, 1-5 August 2022, Eindhoven, The Netherlands, volume 180 of Proceedings of Machine Learning Research, pages 2383 2393. PMLR, 2022a. URL URL httptttt//pproeeaaggg.. mlr, grress/v180/zhoull html
Fangting Zhou, Kejun He, Kunbo Wang, Yanxun Xu, and Yang Ni. Functional bayesian networks for discovering causality from multivariate functional data, 2022b.
Yidong Zhou, Satarupa Bhattacharjee, Cody Carroll, Yaqing Chen, Xiongtao Dai, Jiaming Fan, Alvaro Gajardo, Pantelis Z. Hadjipantelis, Kyunghee Han, Hao Ji, Changbo Zhu, Hans Menn-Geerg Müller, and Jane-Ling Wang. fdapace: Functional Data Analysis and Empirical Dynamics, 2022c. URL httpp:://CAAl R- project. org wwwpppaaaaaapppppee
R package version 0.5.9.
Hongxiao Zhu, Nate Strawn, and David B Dunson. Bayesian graphical models for multivariate functional data. Journal of Machine Learning Research, 17(204):1-27, 2016.
13


Reviews:
[{"review_id": "l9Vy4dlPQE", "reviewer": "Reviewer_ucHi", "summary": "The paper considers the problem of learning causal structure from multivariate functional data that may involve cyclic interactions.\nIt employs an adaptive mapping to a lower dimensional space that retains the relevant causal information, in combination with a Bayesian framework that obtains posterior estimates through MCMC sampling. Effectiveness of the model is compared against several alternatives on synthetic and real-world data.", "questions": "- what is the relation/difference between your model assumptions and e.g. the simple SCMs in Bongers et al. [‘Foundations of structural causal models with cycles and latent variables.’,2018]?\n\n- how do you determine Kj for the lower dimensional space embedding? \n\n- is it true that you essentially treat multiple observations at different time points as independent observations on different system instances in equilibrium?\n\n- in. Table 1, how do you compare non-invariant edges in the equivalence class of e.g. PC/CCD vs. the ground truth? (as ‘half wrong’?)\n- given the causal sufficiency assumption, how come you find bidirected edges in Fig.2?\n\nOther remarks (for lack of a better place to put them):\n- title seems grammatically a bit weird \n- 5 ‘enhance interpretability’ => I understand the practical effectiveness in going via a lower dimensional internal representation, but the output does not reflect this, right? So how does it affect ‘interpretability’?\n- 88 ‘cyclic component’ => this definition seems a bit off as it allows for vertices in a subgraph that are not part of a cycle (so only the full graph G can be a ‘maximal’ cyclic component). I would suggest using the standard term ‘strongly connected component’.\n- 258 ‘simulate posterior samples through MCMC’ => actually I would not mind seeing a bit more on this step in the main paper\n", "limitations": "NA", "rating": 7, "confidence": 4, "soundness": 4, "presentation": 4, "contribution": 3, "strengths": "The paper is exceptionally well written and contains a comprehensive and powerful approach to cyclic causal discovery.\n\nMany of the steps involved are not new by themselves, but brought together effectively to produce an elegant and novel method that shows promising performance. The descriptions and theoretical derivations are concise, clear and consistent, and each of the steps involved follows logically to reduce a challenging problem to an effective, fully Bayesian solution. I particularly liked the way the adaptive basis expansion in 3.2 was incorporated into the overall method.\n\nAssumptions are strong, but clearly and carefully spelled out, leading to the desired identifiability conclusion in Thm2.1.\n\nThe synthetic data experimental evaluation is convincing (although due to the required assumptions inevitably biased in favour of the proposed approach), and the alcoholic EEG application is interesting (though hard to evaluate for non-experts).", "weaknesses": "- although I really like the paper, I think the main contribution lies mostly in the way the various steps are brought together in a clear and coherent way than in fundamentally new insights or approaches. (Nevertheless, I do consider the end result a valuable and novel contribution.)\n- assumptions are quite strong: to (still) see causal sufficiency feature so prominently is disappointing, the ‘disjoint cycles’ should be unnecessary (as admitted by the authors), and ‘non-Gaussianity’ may seem generic but implies identifiability may rely on weak distributional signals (data hungry) and it excludes the challenging complications encountered in the linear Gaussian case that e.g. CCD was specifically designed to handle. \n\n- one confusing aspect (at first) was that in the beginning of the paper the introduction of the time measurements seemed to suggest we are doing time-series analysis, but on closer inspection it is actually closer to standard observational causal inference (with cycles), right?\n\n- some questionable claims, e.g. l.245 ‘strongly prevents false discoveries’ could equally be stated as ‘heavily biased towards sparsity’, and l.307 ‘strong evidence of the effectiveness of FENCE compared to existing methods’ should contain the caveat ‘under the stated assumptions 1-6’. In particular, starting from PCA for LINGAM / CCD seems questionable, where LiNGAM is designed for DAGs, and independence based methods like CCD do not (need to) rely on functional assumptions that are essential for FENCE.\n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "**Response to your questions**:\n\n**[Q1]** Bongers introduces simple structural causal models (SCMs), which allows certain cyclic settings. An SCM is simple if any subset of its structural equations can be solved uniquely for its associated variables in terms of the other variables that appear in these equations. \n\n* **Similarity in assumptions**:  Like our Assumption 2, simple SCMs do not allow for self cycles.\n\n* **Dissimilarity in assumptions**: \n    1. As we focus on linear SEMs, Proposition C.3 in Appendix C of their paper states that a necessary and sufficient condition for a linear SCM to have unique solvability w.r.t. a subset $\\mathcal{L} (\\subseteq V)$ is the invertibility of the matrix $\\mathbf{A}\\_{\\mathcal{L}\\mathcal{L}} = \\mathbf{I} - \\mathcal{B}\\_{\\mathcal{L}\\mathcal{L}}$. A following remark also states that a sufficient condition for invertibility of $\\mathbf{A}\\_{\\mathcal{L}\\mathcal{L}}$ is that the spectral radius of the matrix $\\mathbf{B}\\_{\\mathcal{L}\\mathcal{L}}$ is less than $1$. So for a linear SCM to be simple, this condition needs to be true for all such subsets $\\mathcal{L}$. Although this condition implies the stability assumption (Assumption 3) of our paper, the converse is not true. So our stability assumption is less stringent.\n    2. They allowed confounders  whereas we assume causal sufficiency (Assumption 1).\n    3. One key property of simple SCMs is that the solutions always satisfy the conditional independencies implied by $\\sigma$-separation. So most nonparametric independence test-based methods can be used for simple SCMs by replacing $d$-separation with $\\sigma$-separation. This allows the graph to be identified up to an equivalence class. But, in our paper, additional parametric assumptions (Assumption 4 \\& 5) allows us to identify a unique graph.\n\n\n**[Q2]** \nTo determine $K_j$, we adopted a simple heuristic approach, following Kowal et al. (2017) and Zhou et al. (2022), which selects conservative global value for $K\\_{j} = K$ for all $j \\in {1, \\dots, p}$. Specifically, we impute the functional observations and arrange them into a $(n \\times p) \\times d$ matrix, where $d = |\\cup\\_{i,j} \\mathcal{T}^{(i)}\\_{j}|$ is the size of the union of the measurement grid over all functions. Next, we perform SVD and select the minimum value of $K$ to explain at least 90\\% of the variance. Note that although $K$ is fixed throughout MCMC, the basis functions are inferred (i.e., they are not functional principal components) adaptively to causal structure.\n\nNow, we've implemented another approach to tune $K$ based on WAIC. We tested it in simulations. We found that the $K$ value resulting in the lowest WAIC varies mostly around $\\pm2$ from the value we obtained using the heuristic approach. The graph recovery performance measured by Matthew's correlation coefficient is pretty similar to that of the heuristic approach. We will discuss these in the revised paper.\n\n**[Q3]** We assume our data are drawn from an equilibrium distribution of a dynamic system involving multivariate functions.\nLet $\\mathbf{Y}(\\cdot)[t]=(Y\\_1(\\cdot)[t],\\dots,Y\\_p(\\cdot)[t])^\\top$ denote a vector of functions at time point $t$ where the domain of each function $Y_j(\\cdot)[t]$ is not necessarily time. Note that we have used $()$ to denote the function input/domain and $[]$ to denote the time index of a dynamic system.\nWe consider an AR1-type dynamic system of those functions, \n\n$$\\mathbf{Y}(\\cdot)[t] = \\mathbf{\\mathfrak{B}}\\mathbf{Y}(\\cdot)[t-1] + \\mathbf{f}(\\cdot)$$\n\nwhere $\\mathfrak{B} = (\\mathcal{B}\\_{j\\ell})\\_{j,\\ell = 1}^{p}$ is a matrix of operators and $\\mathbf{f}(\\cdot)$ is a $p$-variate noise process, which has the same domain as $\\mathbf{Y}(\\cdot)[t]$. This recursively leads to,\n$$\\mathbf{Y}(\\cdot)[t] = \\mathbf{\\mathfrak{B}}^{t}\\mathbf{Y}(\\cdot)[0] + \\sum\\_{s = 0}^{t-1} \\mathbf{\\mathfrak{B}}^{s}\\mathbf{f}(\\cdot)$$\nFor the system to eventually reach equilibrium, it is crucial for both $\\mathbf{\\mathfrak{B}}^{t}$ and $\\sum\\_{s=0}^{t-1}\\mathbf{\\mathfrak{B}}^{s}$ to converge as $t$ approaches infinity. The condition essential for this convergence is that the moduli of the eigenvalues of $\\mathbf{\\mathfrak{B}}$ must be less than 1, and none of the real eigenvalues should equal 1. Given this criterion, as $t \\rightarrow \\infty$, $\\mathbf{\\mathfrak{B}}^{t}$ approaches 0, and $\\sum\\_{s=0}^{t-1}\\mathbf{\\mathfrak{B}}^{s}$ tends towards $(\\mathbf{I} - \\mathbf{\\mathfrak{B}})^{-1}$. Consequently, we can deduce that as $t\\to \\infty$,\n$$\\mathbf{Y}(\\cdot)[t] \\rightarrow \\mathbf{Y}(\\cdot) = (\\mathbf{I} - \\mathbf{\\mathfrak{B}})^{-1}\\mathbf{f}(\\cdot)$$ which is equivalent to Equation (2) in our paper. Therefore, we model $\\mathbf{Y}(\\cdot)$ at the equilibrium of a functional dynamic system. \n\n**[Q4]** We favored PC and CCD by counting a non-invariant edge between two nodes as a true positive as long as the two nodes are adjacent in the true graph. \n \n**[Q5]** The bi-directed edges shown in Figure 2 are directed cycles, i.e., $i\\leftrightarrow j$ means $i\\to j$ and $j\\to i$. We will explain it in the revised paper.  \n\n**Response to the weaknesses mentioned**:\n\n* (**Confusion regarding time measurements**) Yes, you are right. By viewing an entire function as one random object, we essentially discover causality based on cross-sectional observational data even though the functional object may be a function of time. In fact, they need not be functions of times; they could have other domains such as space and frequency domains. \n\n* (**Some questionable claims**) Yes, we agree to your point regarding the false discovery; we will change it to \"prevents false discoveries and leads to a sparse network\" in the revised paper. \n    We also agree with you that LiNGAM, PC, and CCD are not the best methods for the cyclic functional causal discovery problem. However, we had to resort to these methods for comparison because of the general lack of existing methods for this problem. \n"}, {"review_id": "nn277o8Fl1", "reviewer": "Reviewer_236D", "summary": "In this paper, the author proposes an operator-based non-recursive linear structural equation based novel causal discovery framework that identifies causal relationships among functional objects in the presence of cycles and additional sampling noises. Furthermore, author demonstrates the effectiveness of proposed method from experiment and theory.", "questions": "1. The proposed year of your baselines fLiNG, LiNGAM, PC, CCD are 2022, 2006, 1991, 1996 respectively, more recent baselines should be considered. According to line 282 and 283, the codes of Lee (2022) and Yang (2022) are not available. However,  fLiNG’s code is also not public (I tried my best to find it but failed). Hence, how do you obtain the results reported in the paper ? \n2. In line 18, you said casual discovery is popular in machine learning. But both in introduction and related work, I do not find relative description to support this standpoint.\n3. How to obtain equation (3) ? It appears to be left-multiplying Equation (1) by  and  respectively. If so, the term in first row should be .\n4. In line 116, Why does j vary from 1 to  ? The range of j should be from 1 to p.\n5. Many symbols are used without definition, like Beta in line 241, Dir and IG in line 254. Some symbols are easy to search for their meanings, while others are not as straightforward.\n6. The equation (7) are not mentioned in Section 3, is it really necessary ?\n7. The name of Section 3 is Model inference. However, nothing but the prior distribution of parameters are introduced, this is confusing.\n8. Reference [Fangting Zhou, Kejun He, and Yang Ni. Causal discovery with heterogeneous observational data] mentions that they do not restrict their model to be acyclic, which is in conflict with your statement in line 56 (Related work section).", "limitations": "No apparent negative societal impacts.", "rating": 5, "confidence": 4, "soundness": 3, "presentation": 2, "contribution": 2, "strengths": "The motivation is clear and the writing is well. The experiment results show that the proposed method significantly outperforms the baselines. The causal identifiability proof is impressive.", "weaknesses": "Please refer to the Questions.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "**Response to your questions**:\n\n**[Q1]** We obtained the fLiNG's code from the authors of that paper. Our research team has also reached out to Yang (2022) but we were unable to obtain their code. The paper by Lee (2022) mainly focused on the statistical theories and hence the code is not available.\n\n**[Q2]** Causal discovery has seen applications in many machine learning areas such as causal representation learning, causal fairness, causal transfer learning, and causal reinforcement learning. \n\n* Bernhard Schölkopf, Francesco Locatello, Stefan Bauer, Nan Rosemary Ke, Nal Kalchbrenner, Anirudh Goyal, Yoshua Bengio \"Toward causal representation learning.\" Proceedings of the IEEE 109.5 (2021): 612-634.\n\n* Tang, Zeyu, Jiji Zhang, and Kun Zhang. \"What-Is and How-To for Fairness in Machine Learning: A Survey, Reflection, and Perspective.\" ACM Computing Surveys 55.13s (2023): 1-37.\n\n* Mateo Rojas-Carulla, Bernhard Schölkopf, Richard Turner, Jonas Peters \"Invariant models for causal transfer learning.\" The Journal of Machine Learning Research 19.1 (2018): 1309-1342.\n    \n* Yan Zeng, Ruichu Cai, Fuchun Sun, Libo Huang, Zhifeng Hao. \"A Survey on Causal Reinforcement Learning.\" \tarXiv:2302.05209 (2023).\n\nWe will add these references in the revision of our paper.\n\n\n**[Q3]** \n\\begin{aligned}\n\\mathcal{P}\\_{j}Y\\_{j} &= \\sum\\_{\\ell}\\mathcal{P}\\_{j}\\mathcal{B}\\_{j\\ell}Y\\_{\\ell}~ + \\mathcal{P}\\_{j}f\\_{j} \\newline\n&= \\sum\\_{\\ell}\\mathcal{P}^{2}\\_{j}\\mathcal{B}\\_{j\\ell}\\mathcal{P}\\_{\\ell}Y\\_{\\ell}~ + \\mathcal{P}\\_{j}f\\_{j} \\newline\n&= \\sum\\_{\\ell}\\mathcal{P}\\_{j}\\mathcal{B}\\_{j\\ell}\\mathcal{P}^{2}\\_{\\ell}Y\\_{\\ell}~ + \\mathcal{P}\\_{j}f\\_{j} \\newline\n&= \\sum\\_{\\ell}\\mathcal{B}\\_{j\\ell}\\mathcal{P}\\_{l}Y\\_{\\ell}~ + \\mathcal{P}\\_{j}f\\_{j}\n\\end{aligned}\n    The second equality is obtained since we assume ${\\mathcal{B}}\\_{j\\ell}=\\mathcal{P}\\_{j}\\mathcal{B}\\_{j\\ell}\\mathcal{P}\\_\\ell$ (stated above equation (3)), which means that causal effects can be fully described within the low-dimensional subspace $\\mathcal{D}\\_{j}$. The third equality follows from that fact that $\\mathcal{P}\\_{j}$ is a projection onto $\\mathcal{D}\\_{j}$ and therefore we have $\\mathcal{P}^{2}\\_{j} = \\mathcal{P}\\_{j} \\forall j\\in$\\{$1, \\cdots, p$\\}. The fourth equality is obtained again by using ${\\mathcal{B}}\\_{j\\ell}=\\mathcal{P}\\_{j}\\mathcal{B}\\_{j\\ell}\\mathcal{P}\\_\\ell$.\n\n**[Q4]** It is a typo. Instead of $j$ in the subscript, it should be $u$, which runs from 1 to $m\\_{j}$. We will fix the typo in the revision of the paper.\n\n**[Q5]** By Beta, Dir and IG, we meant beta distribution, Dirichlet distribution, and inverse-gamma distribution, respectively. We will carefully read our paper and fix all the symbols and short-hands in the revision of the paper.\n\n**[Q6]** We referred to equation (7) in the proof for causal identifiability in Section A of the supplementary material. We also used the term $\\mathbf{\\Phi}(\\mathbf{t})$ defined in equation (7) for the description of Assumption 6 in the paper. \n\n**[Q7]** We agree and we will change the title of Section 3 from \"Model Inference\" to \"Bayesian Model Formulation\". \n\n**[Q8]** We believe the references you mentioned are from lines 54 and 56 in the paper. Since the first author has two papers in 2022, namely 2022a and 2022b, we made a reference to their paper titled 'Functional Bayesian Networks for Discovering Causality from Multivariate Functional Data' (2022b) in line 54. In this particular paper, the author constrains their model to be acyclic. Note that in their paper (2022a), which we cited in a different place, they did consider cyclic graphs."}, {"review_id": "1M3wiKHy7U", "reviewer": "Reviewer_nZkc", "summary": "The authors propose a causal model, i.e. a linear structural equation model for multivariate functional data. The authors' proposed model does not require the dependency structure to conform to a DAG, allowing for modeling cyclic cause-effect relationships. Given certain assumptions on exogenous variables and cycle structure, the authors show the identifiability of their model; a central assumption of which is that the causal interrelations between the variables pertain only to a subset of the function space which the observations inhabit. The authors propose an MCMC sampling procedure to sample from the posterior of the parameters given prior hyperparameters. The authors test their model under various conditions in synthetic data, in addition to a real world application.", "questions": "- In addition to the aforementioned assumptions, using a linear model is another disputable choice, as the authors note. Did the authors experiment with nonlinear generated data to examine the effects of this potential source of model misspecification?\n- Brain imaging is an obvious (and worthy) example for the utilization of such a modeling approach. Are there any other impactful application fields that can benefit from such work?\n- In L99, is it supposed to be $\\mathcal{H}_\\ell$ to $\\mathcal{H}_j$?\n- In L116, is it supposed to be $\\\\{(t_{ju}, X_{ju})\\\\}^{m_j}_{u=1}$?", "limitations": "The authors generally transparent about their modeling assumptions and potential effects thereof. However some of the points I mentioned above would benefit from more elaboration in the final paper.", "rating": 7, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "- The authors' research is soundly motivated, the central observation regarding the importance of being able to handle cyclic relationships in functional data is well-justified, especially given the application areas they consider.\n- The authors' central modeling assumption, that the causal interrelations pertain to a subset of the function space is interesting and has the potential to be utilized and expanded upon future research.\n- Minus some limiting assumptions, the authors' identifiability results and inference procedure is valuable.\n- The authors provide a rigorous exploration of various aspects of their model in comparison with baselines in synthetic data.", "weaknesses": "- The authors make a series of assumptions (1-6) in order to achieve identifiability, with the first one being causal sufficiency. Both the assumption of causal sufficiency and making simplifying assumptions in general is a widely used practice especially when expanding causal analysis to previously unexplored territory, and is usually acceptable to expedite the initial analysis of a new idea. Although it is still understandable that the authors make this assumption, for the kinds of data they are hoping to conduct causal analysis on this assumption unfortunately seems to be especially likely to be violated. For example, it sounds very improbable for measurements made in specific parts of the brain to not have latent confounding. A similar case against disjoint cycles can also be made, while the rest of the assumptions seem less offensive. That being said, I think the authors' contributions justify these limitations.\n- Given the presence of these limiting assumptions, some external validation through real data experiments is called for. The authors' experiment with real data present no external validation (e.g. discovering some known causal relationships among various brain regions).", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "**Response to your questions**:\n\n**[Q1]** Following your suggestion, we have carried out misspecification analysis for our model by considering nonlinearly generated data. We have considered a scenario with number of samples ($n$) $= 100$, number of nodes ($p$) $= 6$ and evenly spaced time grid over $(0, 1)$ of size $d = 100$. The summary measures corresponding to $10$ replicates are given by:\n\n| Method | TPR | FDR | MCC |\n| :--------: |  :--------: |  :--------: |  :--------: | \n| FENCE | 0.27(0.08) | 0.71(0.07) | 0.34(0.07) |\n\nThis poor performance is clearly expected because our modeling assumptions involve linear SEMs. We will add this misspecification analysis result with respect to non linearly generated data in the revised version of the paper. \n\n**[Q2]** Some other application fields that can benefit from our work are (1) finance (e.g., stock price data), (2) mobile health (e.g., wearable device data), and (3) electronic medical records (e.g., longitudinal studies).\n\n**[Q3]** It is a typo. We will correct it in the revision of the paper.\n\n**[Q4]** It is a typo. Instead of $j$ in the subscript, it should be $u$, which runs from 1 to $m\\_{j}$. We will correct it in the revision of the paper. \n\n**Response to the weaknesses mentioned**:\n\n1.(**Assumptions of the model**) We acknowledge the causal sufficiency assumption is limitation of our current work and relaxing it would be an interesting future direction.\n\nFor the disjoint cycle assumption, we provide an example below, which shows the reason why the assumption is needed to uniquely identify the underlying graph.  \nConsider the following two cyclic graphical models,\n\n**Model 1**:\n\\begin{aligned}\n    Y\\_{1} &= 0.95Y\\_{2} + \\epsilon\\_{1} \\newline\n    Y\\_{2} &= Y\\_{3} - Y\\_{4} + \\epsilon\\_{2} \\newline\n    Y\\_{3} &= -1.05 Y\\_{1} + \\epsilon\\_{3}\\newline\n    Y\\_{4} &= -0.1Y\\_{1} + \\epsilon\\_{4}\n\\end{aligned}\n\n**Model 2**:\n\\begin{aligned}\n    Y\\_{1} &= -0.95Y\\_{3} + \\epsilon\\_{1} \\newline\n    Y\\_{2} &= 1.05Y\\_{1} + \\epsilon\\_{2} \\newline\n    Y\\_{3} &= Y\\_{2} + Y\\_{4} + \\epsilon\\_{3}\\newline\n    Y\\_{4} &= -0.1Y\\_{1} + \\epsilon\\_{4}\n\\end{aligned}\n\nBoth Models 1 and 2 have overlapping cycles. Because independent component analysis (ICA) is identifiable up to the column permutation and scaling of the demixing matrix $(\\mathbf{I} - \\mathbf{B})^{-1}$, Models 1 and 2 are equivalent as their respective $\\mathbf{B}$ matrices are row permutation and scaling of each other.\nHence, this example shows that in the presence of overlapping cycles, the underlying true graph can only be identified up to (ICA) equivalence classes. \n\n2.(**Validity of the real data analysis results**) We now discuss the validity of our real data results. In Hayden et al. (2006), it was observed that alcohol-dependent subjects exhibited frontal asymmetry, distinguishing them from the control group. Our own investigation aligns well with these results, as we have identified denser connectivity across various brain regions in the middle and left areas of the frontal lobe among alcoholic subjects, when compared to controls. Furthermore, Winterer et al. (2003) documented coherent differences between alcoholics and controls in the posterior hemispheres, specifically in the temporal, parietal, and occipital lobes. In accordance with their findings, our study provides additional support for this claim, as we have observed heightened activity with several cycles formed in those same regions within the alcoholic group when compared to the control group.\n    \n**References**:\n\n* Hayden, Elizabeth P., et al. \"Patterns of regional brain activity in alcohol‐dependent subjects.\" Alcoholism: Clinical and Experimental Research 30.12 (2006): 1986-1991.\n\n* Winterer, G., et al. \"EEG phenotype in alcoholism: increased coherence in the depressive subtype.\" Acta Psychiatrica Scandinavica 108.1 (2003): 51-60.\n"}, {"review_id": "RNVZm1oXoU", "reviewer": "Reviewer_3JSM", "summary": "This paper discusses the inference of causal relationships between multivariate functions through the proposal of a causal discovery model. The model involves embedding the functional nodes into a lower-dimensional space and separating the structural equation model into two components: the projection onto the space and its orthogonal complement in the Hilbert space defined on the domain. The identifiability of the proposed model is proven based on several assumptions, including disjoint cycles, as described in the Supplementary Materials. The paper presents a model inference method utilizing a fully Bayesian approach. Experimental results using simulated data demonstrate that the proposed method outperforms both conventional causal discovery methods for multivariable functional data and other conventional causal discovery methods. Additionally, insightful observations are made from the application of the proposed method to brain EEG data.", "questions": "Is it possible to consider providing additional discussion on the results of the brain EEG data, in order to help readers ascertain the validity of the findings and comprehend their significance?", "limitations": "The limitation of the proposed method lies in the assumptions required to establish the identifiability of the model.", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "The strongest aspect of this paper is its introduction of a causal discovery model for multivariate functional data, which accommodates the presence of cycles. This is a significant contribution since many multivariate functional datasets, such as EEG data, inherently involve cycles, including feedback loops. Additionally, the paper proves the identifiability of the model and presents a fully Bayesian approach to model inference. The effectiveness of the proposed method is demonstrated through the utilization of both simulated data and real-world EEG data.", "weaknesses": "The discussion of the results of the brain EEG data is inadequate. While there are differences in connectivity between the alcoholic and control groups, readers are unable to determine the significance of these findings because most readers may not possess a background in brain science.\n\nThere may be a minor mistake: from my understanding, the subscript $j=1$ below the superscript $m_j$ in line 116 might actually be $u=1$.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "**Response to your question**:\n\nIn Hayden et al. (2006), it was observed that alcohol-dependent subjects exhibited frontal asymmetry, distinguishing them from the control group. Our own investigation aligns well with these results, as we have identified denser connectivity across various brain regions in the middle and left areas of the frontal lobe among alcoholic subjects when compared to controls. Furthermore, Winterer et al. (2003) documented coherent differences between alcoholics and controls in the posterior hemispheres, specifically in the temporal, parietal, and occipital lobes. In accordance with their findings, our study provides additional support for this claim, as we have observed heightened activity with several cycles formed in those same regions within the alcoholic group when compared to the control group.\n    \n* Hayden, Elizabeth P., et al. \"Patterns of regional brain activity in alcohol‐dependent subjects.\" Alcoholism: Clinical and Experimental Research 30.12 (2006): 1986-1991.\n\n* Winterer, G., et al. \"EEG phenotype in alcoholism: increased coherence in the depressive subtype.\" Acta Psychiatrica Scandinavica 108.1 (2003): 51-60.\n\n**Response to the weaknesses mentioned**:\n\nWe have provided additional validity of our results in the EEG data, as stated in our response above. In addition, we will correct the typo (in line 116) in the revision of our paper.\n"}, {"review_id": "JmQ99ZBW5a", "reviewer": "Reviewer_Y9Hq", "summary": "The paper develops a causal discovery method for directed cyclic graph. The proposed method is a two-step approach which utilizes a Bayesian approach to reduce the dimension of functional data and performs the causal structure learning on the learned embeddings. Experiments on the synthetic and real-world datasets well demonstrate the effectiveness of the proposed method.", "questions": "Please refer to weaknesses.", "limitations": "No\nI have not found the limitations and  broader societal impacts of this manuscript.", "rating": 5, "confidence": 3, "soundness": 3, "presentation": 1, "contribution": 2, "strengths": "1. There are seldom paper studying the causal discovery method on the directed cyclic graph. However, in some particular applications, there may exist cycles in the causal graphs.\n2. The proposed method outperforms baselines in various settings.\n3. The results of the brain EEG data provide some insights on the organization of brain regions.\n", "weaknesses": "1. It lacks some classical and SOTA baselines, e.g., NOTEARS[1], DAG-GNN[2], CSIvA[3]\n[1] Zheng X, Aragam B, Ravikumar P K, et al. Dags with no tears: Continuous optimization for structure learning[J]. Advances in neural information processing systems, 2018, 31.\n[2] Yu Y, Chen J, Gao T, et al. DAG-GNN: DAG structure learning with graph neural networks[C]//International Conference on Machine Learning. PMLR, 2019: 7154-7163.\n[3] Learning to induce causal structure. ICLR 2023.\n2. It is hard to set so many priors for the proposed model in practice.\n3. The presentation of the manuscript could be improved. Fig.1 is hard to understand.\n4. For the results of real data, it is better to also show the results of baselines.\n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "**[Q1]**  Following your suggestions, we have now compared the proposed method FENCE with DAG-GNN and NOTEARS (we could not find publicly available code for CSIvA yet). Both DAG-GNN and NOTEARS were applied to the functional principle components. The results are reported in the table below. We have considered the scenario where the number of samples $n = 150$, number of nodes $p = 40$, and evenly spaced time grid over $(0, 1)$ of size $d = 125$. These results are based on $10$ replicates.\n\n| Method | TPR | FDR | MCC |\n| --------- | :-----: | :-----: | :------: | \n| FENCE | 0.80(0.07) | 0.23(0.06) | 0.83(0.05) |\n| DAG-GNN | 0.34(0.03) | 0.63(0.07) | 0.41(0.08) |\n| NOTEARS | 0.38(0.05) | 0.58(0.05) | 0.47(0.09) |\n\nWe can see that FENCE has outperformed both DAG-GNN and NOTEARS for the above considered scenario, again indicating the ineffectiveness of approaches that use functional PCA + existing multivariate causal discovery methods. \n\n**[Q2]** The default choice for hyperparameters that we took are:\n\n* Causal effect size $\\gamma$: $a\\_{\\gamma} = 1, b\\_{\\gamma} = 1$\n* Noise variance $\\sigma\\_{j}$: $a\\_{\\sigma} = 0.01, b\\_{\\sigma} = 0.01$\n* Variance of the Gaussian mixture distribution $\\tau\\_{jkm}$: $a\\_{\\tau} = 1, b\\_{\\tau} = 1$\n* Number of known adaptive spline bases $R$: $R = 15$\n* Number of unknown basis functions $S$: $S  = 10$\n* Number of mixture components $M$: $M = 10$\n* Mixing coefficient parameter $\\pi\\_{jk}$: $\\beta = 1$\n\nMany of those choices are standard in Bayesian statistics, for example, $a\\_{\\sigma} = 0.01, b\\_{\\sigma} = 0.01$ ensures non-informativeness of the prior on noise variance (i.e., the data will dominate the inference of the noise variance). We have observed that they work well in all of our numerical experiments. From our experience, our method is not highly sensitive to the tuning of these hyperparameters. In Section C.4 of the supplementary materials, we have conducted a sensitivity analysis of those hyperparameters. It shows that our method is relatively robust to hyperparameter choices within the tested ranges. The various options of the hyperparameters that we considered in the sensitivity analysis are listed as follows:\n\n* Causal effect size $\\gamma$: $(a\\_{\\gamma}, b\\_{\\gamma}) \\in \\{(0.1, 0.1), (0.1, 1), (5, 1)\\}$\n* Noise variance $\\sigma\\_{j}$: $(a\\_{\\sigma}, b\\_{\\sigma}) \\in \\{(0.1, 0.1), (0.01, 0.01), (0.1, 1)\\}$\n* Variance of the Gaussian mixture distribution $\\tau\\_{jkm}$: $(a\\_{\\tau}, b\\_{\\tau}) \\in \\{(0.1, 0.1), (0.1, 1), (5, 1)\\}$\n* Number of known adaptive spline bases $R$: $R \\in \\{20, 25, 30\\}$\n* Number of unknown basis functions $S$: $S \\in \\{10, 15, 20\\}$\n* Number of mixture components $M$: $M \\in \\{15, 20, 30\\}$\n* Mixing coefficient parameter $\\pi\\_{jk}$: $\\beta \\in \\{0.1, 2, 5\\}$\n\nWe hope this sensitivity analysis would alleviate your concerns.\n\n**[Q3]** Figure 1 was meant to illustrate the two crucial steps in establishing Theorem 2.1, which we hope helps readers understand the high-level idea of our proof. On the left-hand side (LHS) of the diagram, we depict the hypergraph-like structure that emerges when assuming the existence of disjoint cycles, whereas, on the right-hand side (RHS), we offer a magnified view of the hypernodes (nodes containing simple cycles). Our approach to proving causal identifiability progresses from the LHS to the RHS. That is, we first prove the identifiability of the hypergraph-like structure depicted on the LHS of Figure 1, and then we proceed to establish the identifiability of each simple cycle within every hypernode in the hypergraph. We will add these explanations in the caption of Figure 1 in the revision of our paper.\n\n**[Q4]** We have added a pdf containing the results of the real data analysis for the baselines, together with our response. We will add these results to the revision of our paper."}]

Author rebuttal:
We would like to thank all reviewers for their constructive comments.
This is a global response to several points mentioned by different reviewers while reviewing our work. 

1. **Assumptions of the model**: Some of the reviewers pointed out limitations regarding some of the assumptions that we have taken to show causal identifiability of our model in the paper. We acknowledge that the causal sufficiency assumption is limitation of our current work and relaxing it would be an interesting future direction. As for the disjoint cycle assumption, we provide an example below, which indicates that allowing overlapping cycles could result in non-identifiability.  
Consider the following two cyclic graphical models,

**Model 1**:
\begin{aligned}
    Y\_{1} &= 0.95Y\_{2} + \epsilon\_{1} \newline
    Y\_{2} &= Y\_{3} - Y\_{4} + \epsilon\_{2} \newline
    Y\_{3} &= -1.05 Y\_{1} + \epsilon\_{3}\newline
    Y\_{4} &= -0.1Y\_{1} + \epsilon\_{4}
\end{aligned}

**Model 2**:
\begin{aligned}
    Y\_{1} &= -0.95Y\_{3} + \epsilon\_{1} \newline
    Y\_{2} &= 1.05Y\_{1} + \epsilon\_{2} \newline
    Y\_{3} &= Y\_{2} + Y\_{4} + \epsilon\_{3}\newline
    Y\_{4} &= -0.1Y\_{1} + \epsilon\_{4}
\end{aligned}

Both Models 1 and 2 have overlapping cycles. Because independent component analysis (ICA) is identifiable up to the column permutation and scaling of the demixing matrix $(\mathbf{I} - \mathbf{B})^{-1}$, Models 1 and 2 are equivalent as their respective $\mathbf{B}$ matrices are row permutation and scaling of each other.


2. **Validation of our real data analysis results**: Some of the reviewers wanted to understand the validation of our real data analysis results. So we have found out references from existing literature (Hayden et al. (2006), Winterer et al. (2003)) that supports our claims in the real data analysis results and thus have provided validations of our results. We will add those references in the revised version of the paper. 

3. **Comparison with the baselines for real data analysis**: We have attached a pdf showing the brain connectivity of alcoholic and control groups for the baseline methods fPCA-PC and fPCA-CCD. We could not carry out fPCA-LiNGAM because LiNGAM is not applicable to cases where $q > n$ with $q = Kp$ being the total number of extracted basis coefficients across all functions, a problem that we mentioned in the simulation studies as well. The brain connectivity for fLiNG is shown in the paper Zhou et al(2022). While we could not provide a clear image close to what we have in the paper due to time constraints, some of the key points that we can gather from these brain connectivity plots are:
    1. fPCA-PC could find out denser connectivity across the temporal regions of the brain in the alcoholic group compared to that of the control group. 
    2. Apart from finding denser connectivity in the left temporal lobe of brain in the alcoholic group compared to the control group, fPCA-CCD could find some of the two way cycles ($X \rightarrow Y$ and $Y \rightarrow X$) for the nodes that are adjacent to each other. 

4. **Typos**: We have acknowledged the typos that the reviewers have pointed out to us and we would make those corrections in the revised version of the paper. 

**Reference**:

1. Zhou, Fangting, et al. "Functional Bayesian Networks for Discovering Causality from Multivariate Functional Data." arXiv preprint arXiv:2210.12832 (2022).

2. Hayden, Elizabeth P., et al. "Patterns of regional brain activity in alcohol‐dependent subjects." Alcoholism: Clinical and Experimental Research 30.12 (2006): 1986-1991.

3. Winterer, G., et al. "EEG phenotype in alcoholism: increased coherence in the depressive subtype." Acta Psychiatrica Scandinavica 108.1 (2003): 51-60.

Meta-review: The reviewers appreciated both the quality of the paper and the results, aimed at learning causal models involving cycles. 

Under the usual limiting assumption of causal sufficiency, and further assuming that the causal cycles are disjoint, the use of the projection on a lower dimensional domain is shown to open a new avenue to deal with causal cycles. The dialog with the reviewers convincingly clarified the merits of the approach compared to the state of the art (NoTEARS and DAG-GNN), discussed the hyper-parameter adjustment, and discussed the relevance of the results on EEG data (though these may violate the second assumption).
Decision: Accept (poster)