Witold Wiecek

Frédéric Y Bois

Ghislaine Gayraud

Structure learning of Bayesian networks involving cyclic structures

Keywords: Bayesian inference, Bayesian networks, network inference, structure learning

published or not. The documents may come   L'archive ouverte pluridisciplinaire

Introduction

Large-scale gene expression studies have invigorated interest in exploratory methods for evaluating patterns of association between random variables. The large number of random variables potentially considered and relatively small data sets challenge known structure learning approaches both conceptually and computationally. Graphical models are often used for to represent network structureand for statistical inference [START_REF] Steffen | Graphical Models[END_REF]. They depict the random variables of interest as nodes in a graph, and conditional independence statements about them by presence or absence of graph edges. We focus on 1 Bayesian networks (BNs) which represent probability distributions by means of directed acyclic graphs (DAGs) and are popular in learning structure of biological networks [e.g. 15, 16]. Under additional assumptions, described for example by [23, Chapter 1], directed edges of BNs can correspond to causal relationships between nodes. In BN models the joint probability distribution can be factorized between nodes and evaluated easily.

Yet, there are important cases in biology and other domains of study where we do want to consider cyclic structures such as feedback loops, common in gene transcription regulation networks and other biological networks [2]. In such cases DAGs cannot be used directly, as they do not offer a sound representation of an essential feature of the networks under analysis. Dynamic Bayesian networks (DBNs) offer an alternative, by unrolling cycles in time. However, they multiply the number of nodes by the number of observation times and require dense and extensive data series, as discussed by [START_REF] Ghahramani | Learning dynamic Bayesian networks[END_REF].

We present a different approach to modelling cyclic structures within a graph,. We contract such structures within the graph to derive an associated acyclic graph. The contracted structures are treated as multidimensional random variables. For inference on graph structure, we use a scorebased implementation in the linear Gaussian case. Our approach is fully Bayesian, with scores being Bayesian prior predictive densities. Our procedure uses the factorization property of BNs and, to our knowledge, is novel. We implemented it in Graph sampler, an efficient C language software for simulated network generation and Bayesian inference on network structures.

The informative priors we use, including on cyclical structures, imply that scores between two graphs may differ even if those graphs entail the same conditional independencies. In this sense our approach is an extension to previously proposed Bayesian approaches to network inference, such as the approach by [START_REF] Mukherjee | Network inference using informative priors[END_REF], where the focus is on working with the full posterior distribution.

This paper is organized in two parts. First, Section 2 presents the statistical model, broken down between graph theory background, graph priors (including hyperparameters) and derivation of likelihood (graph score). Then, Section 3 presents choice of hyperparameters and examples of applications in graphs which involve cyclic structures, including computational benefit to Markov Chain Monte Carlo algorithms.

Statistical Model

Methods for learning structure of a Bayesian network (presence or absence of edges between fixed nodes) can be categorized as either test-based methods for conditional independence or score-based methods. The latter tend to give more accurate results, according to [START_REF] Acid | A comparison of learning algorithms for Bayesian networks: a case study based on data from an emergency medical service[END_REF] and [8], but their major disadvantage is the computational cost: since the number of possible graphs to consider grows superexponentially with their number of nodes, exact inference on structure is a hard problem.

In what follows we use a score-based method in a Bayesian framework. For any directed graph G (not necessarily acyclic) we define graph's score s(G|D) conditionally on observed data D. It is proportional to weight of evidence s(D|G) and prior distribution over space of graphs p(G). We define score for any directed graph, not necessarily acyclic.

Term s(D|G) is obtained from the graph's prior predictive density (p.p.d.), that is the data likelihood integrated over all of the model parameters. Such an approach is computationally more efficient than calculating the full posterior function, and parameter values are not needed to make inferences about structure. However, for brevity we still refer to this quantity as posterior, even though some parameters have been averaged out. In the data model, we make use of prior conjugacy, which helps quickly evaluate the p.p.d. A Metropolis-Hastings Markov Chain Monte Carlo (MCMC) algorithm can then be used to sample random graphs from their posterior distribution, cf. [START_REF] Yu | Advances to Bayesian network inference for generating causal networks from observational biological data[END_REF], [27], [9].

As we will show later in this Section, graphs which imply same conditional independencies may not have the same score. For acyclic graphs, this is due to use of informative priors. When cyclic structures are allowed, this inequality may also arise by choice of hyperpriors which can promote or penalise the occurence of cycles. Before proceeding with a description of the statistical model, we introduce graph theory definitions on which the statistical model depends. The rest of this section describes the priors and data likelihood we use.

Graph model

In what follows, we assume that G = (V, A) is a directed graph (with vertex set V and set of directed edges A), of a given order N = |V |. We do not require for G to be acyclic, but edges from a node to itself (autocycles) are not allowed. We use terms graph and network interchangeably.

We say that a graph is strongly connected if for every pair of vertices there exist paths in each direction between the two. A strongly connected compo-nent (SCC) of a graph is a maximal subgraph that is strongly connected. By definition, every cycle is a strongly connected (although not maximal) subgraph. Not all SCCs are cycles, however; e.g. a "flat eight" graph of three nodes A → B → C → B → A is strongly connected but B is traversed twice to get from A to C, hence it is not a cycle. We call single node components ordinary. For each graph we can create a partition of its vertices into sets of strongly connected components. We denote such partition by SCC(G). It can be performed in linear time, as first proposed by [24]. Since most of the strongly connected components which we encounter in structure learning of biological networks are graph cycles, we will also interchangeably use term "cyclic structures" throughout the paper. For any directed graph G, we can create an associated condensed graph, G c , by repeatedly contracting edges (replacing a pair of vertices connected by an edge by a single vertex, retaining all directed edges) in each strongly connected component until each component corresponds to a single vertex. By construction such graph is acyclic. An illustration of the process is provided in Figure 1.

q q q q q A B C D E Graph G q q A {B, C, D} E Condensed G
For DAGs, a Markov equivalence class of graphs is a set of DAGs that have the same skeleton (set of edges without regards to direction) and v-structures (sets consisting of a child and its two parents that are not themselves connected), as shown by [25] in the context of causal inference. Members of the equivalence class encode same conditional independencies. Various algorithms have been proposed to learn Markov equivalence classes of causal graphs, e.g. by [START_REF] Maxwell | Optimal Structure Identification with Greedy Search[END_REF].

When cyclic structures are present and we are working with condensed graphs, we assume that the conditional independencies implied by the Bayesian network are only the ones that are implied by the condensed graph G C . As G C is a DAG, we can take advantage of the Markov property and factorise the joint probability distribution over nodes of the condensed graph G c . We then model the initial G by considering multivariate normal distributions on each of the strongly connected components.

Priors on graph structure

All possible graphs of a certain size are not equally plausible a priori and we should consider that prior knowledge on graph structure in our inferences. Distributions and parameterizations of the following priors which have been previously described by [START_REF] Mukherjee | Network inference using informative priors[END_REF] and by [9]:

• Bernoulli priors on existence of individual (directed) edges, specified by providing a square matrix of edge probabilities. Gene association studies can provide this type of prior knowledge.

• Concordance prior between the graph adjacency matrix and an edge requirement matrix, where each edge is classified as desired, not desired or no preference. This penalizes networks too different from a canonical one (although, tuning this pseudo-prior is not very easy).

• Degree prior on the distribution of node degrees d(v) in the graph, using a power law with parameter γ. The degree distribution of many physical networks appear to follow approximately such a power law [START_REF] Barabási | Emergence of Scaling in Random Networks[END_REF].

• Edge count prior on expected graph size.

• Motif prior on count of triangular feed forward and feedback loops in the network, as discussed in [START_REF] Bois | Probabilistic generation of random networks taking into account information motifs occurrence[END_REF].

All the above priors are specified on the input graph G (and not G c ). To work with cyclic structures, we introduce two additional structure priors on strongly connected components:

Prior on number of strongly connected components. Consider partition of graph G into SCC(G), with ordinary nodes discarded. We define a prior on the number of strongly connected components (of at least two nodes)

|SCC(G)| via a Poisson distribution: |SCC(G)| ∼ P oisson(λ SCC ), with λ SCC > 0.
Prior on the size of strongly connected components. We also define a prior p SCC on the size of all components present in G, using a power law:

p SCC ∝ S∈SCC(G) (|S|) -γ SCC , with γ SCC > 0.
The global prior probability of graph, p(G), is then proportional as a product of all the priors specified.

Data likelihood

We now turn to calculating the (integrated) data likelihood, i.e. the p.p.d. s(D|G). We want to define that score for any directed graph G, including graphs with cycles. Our aim is to fall back on the key feature of Bayesian networks (the factorisation of its associated probability distribution), and we achieve that by condensing G into the acyclic G c .

Let D = x = (x 1 , ..., x N ) denote the observed data on N nodes, where x i is a n-dimensional vector, with n the number of data points per node. When we condense G, we bifurcate its nodes (and corresponding x's) into nodes obtained by contraction (corresponding to strongly connected components of at least two nodes) and ordinary nodes, for which no contraction was needed (corresponding to strongly connected components of single node). Thus we represent x as:

x = (x 1 D , x 2 D , ..., x N 1 D , x 1 L , x 2 L , ..., x N 2 L )
where N 1 denotes the number of ordinary nodes and thus various x i D 's are just relabelled x i 's from the original data set. N 2 represents the number of strongly connected components of size at least 2. For j = 1, 2, ..., N 2 , we denote x j L = (x j 1 , ..., x jm ) a data set of the 'm' node members of the j-th component. Given this partitioning of the graph data, for any given graph G the score can be factorized into a product over ordinary nodes and strongly connected components:

s(D|G) = N 1 i=1 s(x i D |P a(x i D )). N 2 j=1 s(x j L |P a(x j L ))
where P a(•) denotes the parent nodes of "(

•)" in G c . When P a(x i ) = ∅, s(x i |P a(x i )) reduces to s(x i ).
The remainder of this section describes how to obtain the terms s(x i D |P a(x i D )) and s(x j L |P a(x j L )) under a linear Gaussian model.

Integrated likelihood for contracted nodes

Let us consider a strongly connected component of m nodes j 1 , j 2 , j 3 , . . . , j m . As defined above, P a(x j L ) is the set of its parents in the condensed graph G c , i.e. ∪ i∈1,2,...,m P a(x j i ).

We model the distribution of x j L |P a(x j L ) using a linear multivariate Gaussian model; hence setting Y = x j L , the model can be expressed as

Y = Xθ + (1)
where Y is a matrix of dimension (n × m), X is the design matrix of size (n × k) with ones in its first column and P a(x j L ) in the remaining columns so that k = dim(P a(x j L )) + 1, θ is the matrix involving the coefficient terms of dimension (k × m), while is a n × m dimensional matrix = ( 1 , 2 , .., n ) where all i 's are independent and identically distributed according to multivariate Gaussian distribution N m (0, Σ). Under this model, the likelihood L is multivariate normal and can be expressed as:

L(Y, X, θ, Σ) = 1 (2π) nm/2 |Σ| -n/2 (2) exp -tr[Σ -1 (Y -Xθ) t (Y -Xθ)] 2
where tr(•) denotes the trace of •.

For θ and Σ, we consider independent priors, i.e., p(θ, Σ) = p(θ)p(Σ). In order to have an analytically explicit form of the p.p.d. we define an improper (locally uniform) prior on θ, p(θ) ∝ constant. For the prior distribution of Σ, we use an m-dimensional inverse Wishart distribution with hyperparameters (κ, q), denoted W -1 m (κ, q):

p(Σ) = |κ| q 2 |Σ| -(q+m+1) 2 2 mq 2 Γ m ( q 2 ) exp - 1 2 tr(κΣ -1 ) ,
where κ is a positive definite scale matrix and scalar q (degrees of freedom) is strictly positive; Γ m corresponds to the multivariate gamma function. We refer to this prior as constant-Gamma, to distinguish it from other possible models outlined below. The least square estimate θ for the matrix θ, and the sample co-variance matrix A 0 are given by: θ

= (X t X) -1 X t Y A 0 = (Y -X θ) t (Y -X θ) ( 3
)
where θi is the least square estimate for the i -th column of θ. To define θ, we need (X t X) -1 to exist and thus have the constraint k ≤ n.

The joint posterior distribution of both (θ, Σ)|(Y, X) can therefore be expressed as

P (θ, Σ|Y, X) ∝ C 0 • P (θ|Σ, Y, X) • P (Σ|Y, X) (4)
where C 0 , P (θ|Σ, Y, X) and P (Σ|Y, X) are of the form:

C 0 =(2π) m(k-n) 2 2 m(p-q) 2 Γ m ( p 2 ) Γ m ( q 2 ) |κ| q 2 |X t X| -m 2 |κ + A 0 | -p 2 , (5)
P (θ|Σ, Y, X) = |X t X| m 2 |Σ| -k/2 (2π) mk 2 exp -1 2 trΣ -1 (θ -θ) t X t X(θ -θ) , P (Σ|Y, X) = |Σ| -(p+m+1)/2 2 mp 2 Γ m ( p 2 ) exp - 1 2 trΣ -1 (κ + A 0 ) |κ + A 0 | p 2 ,
where p = q + n -k.

When θ and Σ are marginalized out from (4), we obtain(5), which corresponds to the prior predictive distribution of the strongly connected component under the constant-Gamma model, i.e.,

C 0 = s(x j L |P a(x j L )).

Integrated likelihood for ordinary nodes

In the case of an uncondensed acyclic graph, possible forms of integrated likelihood s(x i D |P a(x D i )) have been discussed previously by [9]. With a univariate linear regression model on x i 's parents, using a classical Normal-Gamma conjugate prior (inverse Gamma on the scale and conditional normal on the mean), integrating out these parameters leads to a multivariate Student's t distribution. Zellner and Dirichlet integrated likelihoods are other possible choices and also described therein, together with the choice of likelihood parameters' hyperpriors.

In our case where cyclic structures are allowed, we treat the ordinary nodes as 1-dimensional special cases of the constant-Gamma prior, owing to the fact that the inverse Wishart distribution with parameters q, κ is the multivariate version of the inverse Gamma distribution with parameters (q/2, κ/2). We show the equality of the two prior predictive distributions in Supplement S1.

Model properties and Applications

In this section we will discuss three topics: how choice of hyperparameters impacts graph score in SCC cases; how inference on structure is accomplished with MCMC algorithm; some examples of applications of our approach. The examples will illustrate role that priors and SCCs play in both learning network structure and the computational aspect of inference.

3.1 Likelihood equivalence in constant-Gamma case [START_REF] Geiger | Parameter priors for directed acyclic graphical models and the characterization of several probability distributions[END_REF] discuss conditions under which graphs in an equivalence class will have the same likelihood. A normal model with inverse Wishart prior is such a case. This notion of equivalence can be extended to marginalized likelihoods. [START_REF] Heckerman | Learning Bayesian networks: The combination of knowledge and statistical data[END_REF] present an additional assumption sufficient for marginal likelihood equivalence. It requires that the Jacobian of the one-to-one mapping between two parameters sets associated with two distribution equivalent graphs exists and the priors of the two parameters sets must be equal after applying the change of variables formula.

However, in our case this property no longer holds since we consider independent prior on each parameter set attached to a single term in the likelihood factorization. For the case of a two-node graph, we derive the p.p.d. explicitly in the Supplement S2 and show how s(D|A → B) = s(D|B → A).

In the case of Markov-equivalent DAGs the differences in graph score are due to sampling variance and tend to 0 with increasing n. The score of SCCs will differ from equivalent DAGs due to marginalisation of the likelihood and the difference doesn't tend to 0 with growing sample, but rather depends on the choice of hyperparameters of the inverse Wishart prior distribution of Σ.

Hyperparameters needed to evaluate the graph score under the constant-Gamma model are the scale matrix κ and degree of freedom q. We use inverse Wishart prior as it is conjugate and allows us to obtain the prior predictive distribution analytically. Additionally, this prior, when informative, can be interpreted in terms of equivalent sample size.

If X ∼ W -1 (κ, q) then E(X ij ) = κ ij q -m -1 , when q > m + 1, (6)
Var(X ij ) = (q -m + 1)κ 2 ij + (q -m -1)κ ii κ jj (q -m)(q -m -1) 2 (q -m -3)
, when q > m + 3.

However, there are known issues with inverse Wishart prior: it implies relationships between variances and covariances and uses a single parameter (q) to describe precision on all parameters. When q > 1 the prior may be biased when the true variance is low, even with large sample sizes, as discussed by [11]; see also [START_REF]Bois' work was funded by the Horizon 2020 project[END_REF] for a simulation study.

In inference on variance-covariance matrices, it is typical to assume κ = I m (identity matrix of order m) and q = m + 1. In such case prior marginal distributions of correlations are uniform on (-1, 1). However, in case of covariance in SCC data this goes against our intuition: typically a priori we assume that nodes of an SCCs are going to be strongly correlated; exactly how strongly depends on context and objectives of analysis. By default we propose to set q = m + 1 and κ to 1 on diagonal elements and to 0.5 on offdiagonal. This creates a monotonic prior on correlation and ensures higher score for SCC than all DAG graphs when the true correlation is higher than 85%-90%. This choice is explored and explained below.

Let us define d as the difference in log scores between a m-node graph where all nodes form an SCC (G SCC ) and a DAG with no conditional independencies (G DAG ):

d = log s(G SCC |D) -log s(G DAG |D),
conditional on the same data D. Positive d's indicate that SCC scores higher than the DAG.

We now briefly explore the behaviour of d in SCCs of different sizes by means of simulated data and show that it is predictable in ways that may be useful in practical applications. For all the examples presented in this section we generate n = 10, 000 draws of data of m nodes from multivariate normal distribution with means 0 with each node having fixed variance σ 2 and same correlation ρ with all other nodes. We start with σ 2 = 1 and vary ρ between 0 and 1. We compare the SCC against one DAG only as at large n all DAGs in the equivalence class will have almost identical scores.

As mentioned, it is typical to set q = m + 1 and κ = I m . Since q is responsible for the precision of the prior and can be interpreted in terms of sample size equivalence, setting a low q is a good default choice. However, even with the default value q = m + 1, κ has an impact on choice between SCCs and DAGs. Setting κ = I m in this situation leads to uniform priors on correlation. This can be desirable in inference on unknown variancecovariance matrix, but our goal in this case is to distinguish cyclic structures from DAG parts of the network: therefore a sensible choice would be to have a prior that "favours" loops when correlation is higher. We illustrate behaviour of d as a function of the off-diagonal elements of κ in two panels of Figure 2. With κ = I m the sign of d is not consistent, but when the off-diagonal elements of κ are set to 0.5 everywhere, d is positive when the true correlation in data exceeds 85%-90% threshold. Therefore we use 0.5 as the default choice of prior. Difference grows larger as m increases. The behaviour of d is also sensitive to variance of random variables. If κ is misspecified, the SCCs are always preferred for low variances and DAGs are always preferred for high variances. However, standardising the inputs can solve this problem, as will scaling κ by sampling variances of each node. This can be done automatically in software implementations and in both cases will "bring back" the behaviour of d to exactly what is seen in Figure 2.

As indicated by Equation 6, we can put a prior on correlation between two elements to any mean ρ by setting off-diagonal elements of κ to (q -m -1)ρ, and any variance by adjusting q. In practice we can use this to manipulate the sign of d, thus allowing us to choose the level of correlation at which SCCs will be chosen over DAGs different from the 85%-90% threshold. This is illustrated in Figure 3.

MCMC algorithm for inference

The computer code needed to perform all of the examples has been implemented in the latest version of the graph sampler software for MCMC inference on graphs, previously introduced by [START_REF] Bois | Probabilistic generation of random networks taking into account information motifs occurrence[END_REF]. Written in ANSI-standard C language, the full software is freely available at www.nongnu.org/graphsampler under the terms and conditions of the GNU General Public License, as published by the Free Software Foundation.

Graph sampler uses Metropolis-Hastings algorithm to sample graphs from Cor(X i ,X j ) d / n q q = m + 2 q = m + 10 q = m + 100 For DAGs, score of proposal is then evaluated by calculating difference in scores on the child node in the proposed addition or deletion. In all cases convergence to the target distribution can be checked by calculating Gelman-Rubin statistic [START_REF] Gelman | Inference from iterative simulation using multiple sequences[END_REF] on chains of graph adjacency matrices. Convergence check function is included as part of the software. When cyclical structures are allowed, the algorithm is modified to take into account situations where the condensed graph changes, (i.e. SCCs are created or deleted. Multiple nodes are affected in such situations and need to have their scores recalculated. We devised an additional decision rule to only condense graph (using Tarjan algorithm) when necessary and recalculate likelihood on the minimal set of nodes that may be affected by additions and deletions. It is presented in Supplement S3.

The MCMC approach yields a set of s graphs sampled from the posterior distribution. We represent them by their adjacency matrices A (1) , . . . , A (s) . Such a sample can be used to approximate posterior probabilities of occurence of edges or motifs. For example, the probability of an edge from i to j, p ij , is obtained by calculating pij = A ij /s. However, such probabilities have to be treated with caution when cyclic structures are allowed. Depending on the objectives of analysis, we can either be interested in pij defined as above or the probability of i and j being part of the same SCC (p SCC ij ) or of i being parent of j, but not in the same SCC (p DAG ij ).

MCMC convergence in SCC setting

In MCMC algorithms for graph inference which operate only by adding or removing edges at each step, reversing the direction of an existing edge can be difficult. It requires two operations: a deletion followed by an addition. The first step will often (e.g., in situations where two nodes are highly correlated) have an extremely low probability. Using tempered MCMC methods can solve this problem [see 5] but requires fine tuning of the tempering algorithm. Using SCCs provides a simpler solution: an addition (creating an SCC) followed by a deletion. Thus for some problems, allowing cyclic structures can be helpful even if we know that the true network is acyclic as it can avoid traversing these "probability wells".

(1)

q q q A B C s = -483 (2) q q q A B C s = -478.8 (3) q q q A B C s = -477.3 (4a) q q q A B C s = -593.8 (4b) q q q A B C s = -451 (5) q q q A B C s = -437.7
Figure 4: "Inverting a v-structure" in five steps. For each graph value s is s(D|G). Leftmost graph is the starting graph for MCMC and the true generative mechanism is the rightmost. While deleting some edges such as CA can be "easy" in terms of log likelihood difference (steps 2 and 3), deletion of the CB edge (step 4a) is difficult as it leads into a "probability well". Adding the edge BC (4b), which results in an SCC, offers a way to reach (5) while avoiding the well.

We illustrate this with a simple example of inverting a v-structure. Assume X A ∼ N (0, 1), X B ∼ N (0, 1) and X C = X A -5X B + with ∼ N (0, 1). We draw 100 realisations of each random variable. Assume the MCMC sampler starts from a graph model A ← C → B with score of -472.5. Assuming that we are working with DAGs only, any path to the true generating graph requires removal of CB edge. This is shown in Figure 4.

The well is a score difference of around 100 (therefore on average we will need e 100 Metropolis-Hastings proposals to remove CB). Using SCCs easily circumvents this by creating an SCC involving the B ↔ C SCC before deleting CB.

Linear model with additive noise

Our first example is straightforward, but designed to be difficult to correctly estimate. For this, we slightly expanded the graph from Figure 1 by adding node P , a parent to A, and Q, a child to E. We assumed a linear relationship and generated 100 draws for each node as follows: for j-th node, i-th generated value x j (i) = P a x j (i) + ij , where P a x i are data for the parent node of X j (for node P we set the mean to zero) and ij are i.i.d. with N (0, 5) for all i and j. For the SCC (nodes B, C and D) we used multivariate Gaussian distribution with same means and variances (equal to 5) and pairwise correlations of 0.9. This way, all of generated data was very highly correlated, making it difficult to distinguish between different graphs using likelihood alone, even under the correct assumption about data generating mechanism being a Gaussian linear additive noise model.

First, in the absence of prior information (Figure 5A) we did not succeed in retrieving the data generating graph and many superfluous edges were found. (Although we usually prefer to work with edge probabilities, for clarity of presentation we only show the best scoring graph here.) Including an informative prior on the out-degree (power law with γ = 3) and size of SCCs (no larger than 3) enabled us to detect the SCC and the undirected edges correctly (Figure 5B). Lastly, adding information on the first cause, i.e., enforcing (through a Bernoulli prior) P a(P ) = ∅, allowed us to retrieve the data generating graph (Figure 5C). q q q q q q q

P A B C D E Q A q q q q q q q P A B C D E Q B q q q q q q q P A B C D E Q C Figure 5
: Graph structure inference when likelihood is the same as the generative model. Highest scoring graphs for: (A) A model without informative priors. (B) A model with prior on out-degree and SCC size. (C) A model with additional prior on the "first cause" (no edges allowed into node P ).

The last two steps illustrate two difficulties with learning network structures. First comes the problem of detecting dependencies from data, which can be helped by putting a strong prior on the types of structures expected to occur in the graph (in this case degrees and SCC sizes). Even if we suc-ceed in this, we are still left with multiple candidate graphs: in this case the condensed graph is a path from P to Q which is Markov-equivalent to a path from Q to P . Only the addition of a prior on whether P or Q is a probable cause can help us retrieve the true network.

As discussed, the choice between SCCs and DAGs is highly sensitive to correlation. We repeated the simulation using the same data-generating mechanism, but with a pairwise correlation of the SCC nodes equal to 0.5 instead of 0.9, DAGs were then preferred when using uninformative structural priors and the best graph resembled Figure 5(C) but without an SCC. Setting q = m + 11 and scaling the "default" κ appropriately is enough to bring up an SCC again.

SCCs detection in a 50-node linear model

In the second simulated study, we generated batches of 100 DAGs of 50 nodes by randomly permuting nodes and drawing each edge e ij with probability of occurrence at 5% if j > i (to avoid SCCs). We then created two threenode SCCs in each graph by adding all possible edges between two groups of three randomly selected nodes. Example of such graphs are presented in Figure 6. Each graph was then used as data generation mechanism for 100 data values for each node, according to a normal linear model with regression coefficients set to 1. That is, for j-th node, the i-th generated value x j (i) = X k ∈P a(X j ) x k (i) + ij , with ij ∼ N (0, 1) i.i.d. for all i and j. For SCCs the distribution was multivariate normal, with correlation between any two nodes fixed at 0.5 or 0.9, to benchmark performance in two different cases.

For inference, we compared four prior assumptions on the inverse Wishart parameters by varying q and setting κ to the desired correlation and scaling it appropriately (see Equation 6). We selected a "default" uninformative prior q = m + 1 assuming a within-SCC correlation of 0.5; a q = m + 11 prior (equivalent to 10 data points) assuming a correlation of 0.5; a m + 11 prior assuming a correlation of 0.9; and a m + 101 prior (equivalent to 100 data points) assuming a correlation of 0.5. A summary of these combinations is given together with results in Table 1.

In each case, we set a Bernoulli prior on the probability of edge occurrence p(e ij ) = 0.025 when i = j (instead of 5%, as in graph generation half of offdiagonal edges were not allowed) and constrained the size of SCCs to be at most three, but imposed no more priors.

For each combination of "true" correlation and prior assumptions we generated data and ran MCMC inference 100 times. For each graph, we used 20 millions MCMC iterations, discarding first ten million, to infer on its q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q Figure 6: Examples of three graphs (out of 800) from which data was generated for this performance benchmark. All graphs were randomly generated with the same settings. White nodes belong to (three-node) SCCs. structure. We only assessed MCMC convergence on a few selected graphs, but assumed that such run length was adequate given the simple nature of the problem and that we are interested in relative, not absolute, performance. For each of these runs, the probability of occurrence of SCCs was calculated from a sample of 100 adjacency matrices drawn from the MCMC chain. We report the area under the receiver operating characteristic curve (AUROC). It is the same as described in [START_REF] Marbach | Revealing strengths and weaknesses of methods for gene network inference[END_REF] -briefly, the k possible edges in the graph were ordered by probabilities obtained from MCMC and we calculated sensitivity and specificity k times, assuming that 1, 2, . . . , k first edges occur and the rest do not. Note that perfect prediction (AUROC = 1) is impossible in this example, as we calculate our score for directed graphs and not equivalence classes. For SCCs, we only assessed sensitivity (as the AUROC statistic captures overall specificity well), by calculating a probability that the "true" SCCs are present in the MCMC results. Under our definition we needed to "detect" all three nodes of the SCC to count as a success. Table 1 presents results for both AUROC and "SCC sensitivity"; results are averaged over 100 inferences for each row.

Generally, the sensitivity and specificity (AUROC) of the score-based method is good under this simple generative model. However, the detection of SCCs is low with "default" settings, with about 10% success rate. (Note that given that the equivalence class for three-node SCC is of size seven, i.e. SCC and six DAG configurations, so we would expect success rate of about 14% assuming equivalence of scores within class.) We can improve this by introducing informative priors. Generally highly correlated SCCs are easier Table 1: Sensitivity and specificity of the scoring method in retrieving true network structure.

True Correlation Prior AUROC SCC Pr 0.50 q = m + 1; Cor = 0.5 0.93 0.12 0.50 q = m + 11; Cor = 0.5 0.91 0.26 0.50 q = m + 11; Cor = 0.9 0.86 0.01 0.50 q = m + 101; Cor = 0.5 0.91 0.52 0.90 q = m + 1; Cor = 0.5 0.89 0.10 0.90 q = m + 11; Cor = 0.5 0.89 0.40 0.90 q = m + 11; Cor = 0.9 0.87 0.95 0.90 q = m + 101; Cor = 0.5 0.91 0.39 to detect successfully, but using an informative prior helps. Misspecification of prior does not seem to overly impact the overall (AUROC) performance, but does affect the detection of SCCs. That is most salient in the case of a high correlation prior, as illustrated in the bottom left panel of Figure 7.

When assuming a correlation of 0.50, the variability in success rate across graphs is large (with a peak around probability of 50%, corresponding to cases where one SCC has been identified perfectly and the other one not at all), but with a prior on correlation equal to 0.90 the behaviour is completely different.

Discussion and Conclusion

We proposed a model to represent cyclic structures within Bayesian networks. Our model offers an alternative way of describing joint probability distribution and performing network inference without apparent computational drawbacks. In our approach, SCCs are condensed to form multivariate nodes, which are still embedded in an acyclic Bayesian network. We can therefore factorise the likelihood, a key computational advantage of Bayesian networks. We use a score-based approach in a fully Bayesian setting. A posterior sample of graphs is obtained by MCMC sampling. This allows us to integrate prior knowledge on presence of edges, degrees, acyclic motifs, occurrence of SCCs etc. The placement of informative priors on network structure also brings faster convergence of MCMC sampling (if the data are not conflicting with the prior) by putting soft constraints of the size of the set of likely graphs.

The likelihood model we present is an additive linear model with Gaussian noise. Such model allows us to easily compute score by integrating out q = m + 11; Cor = 0.9 q = m + 101; Cor = 0.5 q = m + 1; Cor = 0. parameters. The only (arbitrary) constraint imposed by our Gaussian model is that the number of parents for all members of an SCC has to be less than the number of data points for each node. In the future other models for likelihood or other scoring functions should be explored. We also note that in the present form the model can only account for time as an additional linear term in regression, although a dynamic version of it might be workable. Use of SCCs with discrete random variables should also be explored.

Many alternative methods for characterising dependencies in graphs containing cycles have been proposed, including reciprocal graph models based on work by [START_REF] Koster | Markov properties of nonrecursive causal models[END_REF] [see paper by 22, for recent application] or a heuristic algorithm approach to learning cycles from experimental data by [17]. Our work differs from statistical models for purpose of learning causal relationships in observational data, as under our model Markov-equivalent graphs can have different scores due to choice of priors and hyperpriors relating to SCCs.

Under the proposed model, detection of SCCs is sensitive to the choice of hyperparameters. Informative priors can be used to promote or suppress occurence of SCCs in the posterior. We can choose priors to favour SCCs over DAG structures even when correlation is lower than the threshold visible in Figure 2. This may be useful in applications where we know a priori that cycles are present or simply wish to describe joint distribution differently. However, in our model, a limitation in setting informative priors is tied to the properties of W -1 distribution, where the precision of variances and correlations is governed by a single parameter, q. Here again we must decide between standardisation and flexibility.

Our simulation study with a small linear model also shows the importance of prior choices on the inference. Note that pure likelihood-based inference amounts to placing only an indifferent Bernoulli prior on the adjacency matrix, and would bring the same inefficient inference as in Figure 5B. In the case of a larger network, the inference scales up well, but SCCs typically have a 50% chance to be detected. Finally, besides substantive applications, allowing for SCCs can reduce computation time (by reducing the number of nodes) even for underlying DAGs, and improves convergence by easing edge reversals. The practical impact of those computational benefits should be explored in greater detail in the future.

Figure 1 :

 1 Figure 1: Graph G with 5 nodes where nodes B, C and D form a cycle and are contracted to a single node after condensing the graph. Here SCC(G) = {A, {B, C, D}, E}

Figure 2 :

 2 Figure 2: Difference d between SCC and DAG graph scores (divided by n data) as a function of true correlation between pairs of nodes (x axis) and number of nodes m. The hyperparameters are set to q = m + 1 and κ = I m in the left panel. In the right panel we change κ ij = 0.5 for i = j.

Figure 3 :

 3 Figure 3: Difference d between SCC and DAG graph scores (divided by n data) under informative priors on correlation. κ is scaled according to q to yield desired mean correlation: 0.25 in the left panel, 0.5 in the middle, 0.75 in the right panel.

Figure 7 :

 7 Figure 7: How the probability of detecting true SCCs changes with "true" correlation (differently coloured histogram bars) and four priors (panels of the figure). Peaks at 0.5 correspond to cases where only one of two SCCs has been detected.