CausNet : Generational orderings based search for optimal Bayesian networks via dynamic programming with parent set constraints

Optimal Bayesian Network (BN) Structure Discovery is a method of learning Bayesian networks from data that has applications in wide variety of areas including epidemiology (see e.g. [1], [2], [3], [4]). Disease pathways found using BN leading to a phenotype outcome can lead to better understanding, diagnosis and treatment of a disease. The challenge in finding an optimal BN lies in the super-exponential complexity of the search. Dynamic programming can reduce the complexity to exponential, but still the number of features/nodes being explored remains very small - $10$ to $20$ nodes - and only by restricting the maximum number of parents for each node. To alleviate the challenge of high dimensionality, we implement a dynamic programming algorithm with parent set constraints. We use what we call ‘generational orderings’ based search for optimal networks, which is a novel way to efficiently search the space of possible networks given the possible parent sets.

Most current algorithms do not work for both continuous and categorical nodes, and none works for survival data as the outcome node. We implement support for both continuous and categorical data, and categorical as well as survival outcomes. This is especially useful for disease modeling where mixed data and survival outcomes are common. We also provide options of two common scoring functions and also output multiple best networks when available.

Our main novel contribution aside from providing software is the revision of the Silander algorithm 3 [5] to incorporate possible parent sets for a much more efficient way to explore the search space as compared to the original approach based on lexicographical ordering. In doing so, we cover the entire constrained search space without searching through networks that don’t conform to the parent set constraints.

In this section, we review the necessary background in Bayesian networks and the Bayesian network structure discovery problem (for more background on these topics see, for example, [6], [7]).

A Bayesian network (BN) is a probabilistic graphical model that consists of a labeled directed acyclic graph (DAG) in which the vertices $V=\{v_{1},...,v_{p}\}$ correspond to random variables and the edges represent conditional dependence of one random variable on another. Each vertex $v_{i}$ is labeled with a conditional probability distribution $P(v_{i}|parents(v_{i}))$ that specifies the dependence of the variable $v_{i}$ on its set of parents $parents(v_{i})$ in the DAG. A BN can also be viewed as a factorized representation of the joint probability distribution over the random variables and as an encoding of conditional dependence and independence assumptions.

A Bayesian network $G$ can be described as a vector $G=(G_{1},...,G_{p})$ of parent sets: $G_{i}$ is the subset of $V$ from which there are directed edges to $v_{i}$. Any $G$ that is a DAG corresponds to an ordering of the nodes, given by the ordered set $\{v_{\sigma_{i}}\},i\in\{1,2,..,p\}$, where $\sigma$ is a permutation of $[p]$ - the ordered set of first $p$ natural numbers, with $\sigma(i)=\sigma_{i}$. A BN is said to consistent with an ordering $\{v_{\sigma_{i}}\}$ if parents of $v_{\sigma_{i}}$ are a subset of $\{v_{\sigma_{j}}\}$, if $j<i$, i.e. $G_{\sigma_{j}}\subseteq\{v_{\sigma_{j}}\},j<i$.

One of the main methods for BN structure learning from data uses a scoring function that assigns a real value to the quality of $G$ given the data. For finding a best network structure, we maximize this score over the space of possible networks. Note that we can have multiple best networks with the same score. Scoring functions balance goodness of fit to the data with a penalty term for model complexity. Some commonly used scoring functions are BIC/MDL [8], BDeu [9], and BGe [10] [11]. We use BIC (Bayesian information criterion) and BGe (Bayesian Gaussian equivalent) scoring functions as two options for using Causnet. BIC is a log-likelihood (LL) score where the overfitting is avoided by using a penalty term for the number of parameters in the model, specifically $p\ln(n)$, where $n$ is the sample size. The BGe score is the posterior probability of the model hypothesis that the true distribution of the set of variables is faithful to the DAG model, meaning that it satisfies all the conditional independencies encoded by the DAG, and is proportional to the marginal likelihood and the graphical prior[10] [11].

CausNet uses the dynamic programming approach to finding a best Bayesian network structure for a given dataset. The idea of using dynamic programming for exact Bayesian network structure discovery was first proposed by Koivisto & Sood [12], [13]. Recent work using dynamic programming, includes that by Silander and Myllymäki in [5] and by Singh and Andrew in [14].

We closely follow the algorithm proposed by Silander and Myllymäki (SM algorithm henceforth) and make heuristic modifications focusing on finding sparse networks and disease pathways. This is achieved by dimensionality reduction with parent set identification, and by restricting the search to the space of ‘generational’ orderings rather than lexicographical ordering as in the original SM algorithm.

Finding a best Bayesian network structure is NP-hard [15]. The number of possible structures for $n$ variables is $\mathcal{O}\left(n!2^{{n\choose 2}}\right)$, making exhaustive search impractical. So the dynamic programming algorithms are feasible only for small data sets; usually less than $20$ variables; the number of variables can be increased somewhat by using small in-degree (maximum number of parents for any node). But bounding the in-degree by a constant $k$ does not help much, a lower bound for number of possible graphs still being at least $n!2^{kn\log n}$ (for large enough $n$) [12]. We introduce parent set identification and ‘generational’ orderings based search to reduce the search space and thus scale up the dynamic programming algorithm for higher number of variables.

The dynamical programming approach uses the key fact about DAGs that every DAG has at least one sink. The problem of finding a best Bayesian network given the data $\cal{D}$ starts with finding a best sink for the whole set of nodes. Then that node is removed and the process is repeated recursively for the remaining set of nodes. The result is an ordering of the nodes $\{v_{\sigma_{i}}\},i\in\{1,2,..,p\}$, from which the DAG can be recovered. Denoting the best sink by $s$, and the score of a best network with nodes $V$ by $bestscore(V)$, and the best score of $s$ with parents in $U$ by $bestScore(s,U)$, the recursion is given by the following relation :

To implement the above recursion, the idea of a local score for a node $v_{i}$ with parents $parents(v_{i})$ is used, which we get using a scoring function. The requirement for a score function $score(G)$ for a network $G$ is that it should be decomposable, meaning that the total score of the network is the sum of scores for each node in the network, and the score of a node depends only on the node and its parents. Formally,

where the local scoring function $localscore(x,y)$ gives the score of $x$ with parents $y$ in the network $G$. In a given set of possible parents $pp_{i}$ for node $v_{i}$, we find the best set of parents $bps_{i}$ which give the best local score for $v_{i}$, so that

Now the best sink $s$ can be found by Eq. 5, and the best score for a best network in $V$ can be found by Eq. 6.

In Figure 1, the subset lattice shows all the paths that need to be searched to find a best network. Observe that each edge in the lattice encodes a sink, so that each path also encodes an ordering on the $4$ variables, e.g. the rightmost path encodes the reverse-ordering $\{1,2,3,4\}$. There are a total of $4!$ paths/orderings to be searched. Now suppose we knew the best score corresponding to each edge in all the paths, meaning the best score for the sink corresponding to that edge with the best parents from the subset at the source of that edge. Then naive depth/width first search would compare the sum of scores along all paths to get the best network. In dynamic programming, we proceed from the top of the lattice. Ignoring the empty set, we start with finding the best sink for each singleton in the first row which trivially is the singleton itself. Next, we find the best sink for the subsets of cardinality $2$ in the second row using the edge best scores. And we continue all the way down. Suppose we get the best sinks sublattice as in Figure 2, then the best network is given by the only fully connected path, and corresponds to the reverse-ordering $\{4,1,2,3\}$.

Now, using the dynamic programming approach of SM algorithm, finding the best Bayesian network structure using basic CausNet, without restricting the search space, has the following five steps.

1. 1.

   Calculate the local scores for all $p2^{p-1}$ different (variable, variable set)-pairs.

2. 2.

   Using the local scores, find best parents for all $p2^{p-1}$ (variable, possible parent set)-pairs.

3. 3.

   Find the best sink for all $2^{p}$ variable sets.

4. 4.

   Using the results from Step 3, find a best ordering of the variables.

5. 5.

   Find a best network using results computed in Steps 2 and 4.

The extensions in this base version of Causnet - possible parent sets identification, phenotype driven search, and search space based on ‘Generational orderings’ reduce the search space.

We compare CausNet with three other methods that have been widely used for optimal Bayesian network identification to infer disease pathways from multiscale genomics data. The first method is Bartlett and Cussens’ GOBNILP [16], an integer learning based method that’s considered state-of-art exact method for finding optimal Bayesian network. The other two methods are BNlearn’s Hill Climbing (HC) and Max-min Hill Climbing (MMHC) ([17], [18]), which are both widely used approximate methods, see e.g. [19], [20]. Hill-Climbing (HC) is a score-based algorithm that uses greedy search on the space of the directed graphs [21]. Max-Min Hill-Climbing (MMHC) is a hybrid algorithm [22] that first learns the undirected skeleton of a graph using a constraint-based algorithm called Max-Min Parents and Children (MMPC); this is followed by the application of a score-based search to orient the edges.

We simulated Bayesian networks by generating an $N$ x $p$ data matrix of continuous Gaussian data. The dependencies are simulated using linear regression with the option to control effect sizes. Some number of the $p$ nodes were designated as sources ($p_{1}$), some intermediate ($p_{2}$), and some sinks ($p_{3}$), the remainder ($p_{0}$) being completely independent. The actual DAGs of the $p_{1}+p_{2}+p_{3}$ nodes vary across replicates. The requirement for being a DAG is implemented using the idea of ordering of vertices. We pick a random ordering of a randomly chosen subset of $p$ vertices. Then enforcing each vertex to have parents only from the set of vertices above itself in the ordering guarantees a DAG.

The False Discovery Rate (FDR) and Hamming Distance are used as the metrics to compare the methods. With $FP$ defined as the number of false positives and $TP$ defined as the number of true positives, FDR is defined as :

Controlling for the false discovery rate (FDR) is a way to identify as many significant features (edges in case of BNs) as possible while incurring a relatively low proportion of false positives. This is especially useful metric for high dimensional data and for network analysis ([23]).

The Hamming distance between two labeled graphs $G_{1}$ and $G_{2}$ is given by $|\{\left(\left(e\in E(G_{1})\&\left(e\not\in E(G_{2})\right)\right)or\left((e\not\in E(G_{1})\&e\in E(G_{2})\right)\right)\}|$, where $E(G_{i})$ is the edge set of graph $G_{i}$. Simply put, this is the number of addition/deletion operations required to turn the edge set of $G_{1}$ into that of $G_{2}$. The Hamming distance is a measure of structural similarity, and forms a metric on the space of graphs (simple or directed), and gives a good measure of goodness of a predicted graph ([24], [25]). In the context of predicted and the truth graph, with $FP$ defined as the number of false positives and $FN$ as the number of false negatives, Hamming Distance is defined as :

As the first set of simulations using parent set identification using correlation tests, we ran simulations using multiple replicates of networks, the first with $p=10,20,40,50,60,100$, and $N=500,1000,2000$. Figure 4 shows the plot of average FDR for different values of $p$ and $N$, and their linear trend with BIC scoring for CausNet.

We can see that for lower values of $p$, Gobnilp has the lowest values of FDR with CausNet the second best, but CausNet performs the best for higher values of $p$. The results for the BGe scoring for CausNet are similar qualitatively. In the tables 1 and 2, we show the average $FDR$ across the 9 combinations of $N$ and $p$, both with and without taking directionality into account. The results for CausNet with the choice of scoring function—either BGE or BIC—are given. The results show that our method performs very well compared with the methods considered.

For the number of variables up to $40$ (Table 1), on the metric of FDR, CausNet performs second best after Gobnilp for directed graphs and the best for undirected graphs. This is true for both scoring methods - BIC and BGE. For the number of variables between $50$ and $100$ (Table (Table 2)), our method performs the best for both directed and undirected graphs, using either of the two scoring methods, BIC and BGE.

Figure 5 shows the plot of average Hamming Distance for different values of $p$ and $N$, and their linear trend with BIC scoring for CausNet. We can see that for lower values of $p$, Gobnilp has the lowest values of Hamming Distance with CausNet the second best for most part, but CausNet performs the best for higher values of $p$. The results for the BGe scoring for CausNet are the same qualitatively.

We applied our method to recently published data on gene expression in relation to ovarian cancer prognosis among participants in a clinical trial [26]. The aim of this study was to develop a prognostic signature based on gene expression for overall survival (OS) in patients with high-grade serous ovarian cancer (HGSOC). Expression of 513 genes, selected from a meta-analysis of 1455 tumours and other candidates, was measured using NanoString technology from formalin-fixed paraffin-embedded tumour tissue collected from 3769 women with HGSOC from multiple studies. Elastic net regularization for survival analysis was applied to develop a prognostic model for 5-year OS, trained on 2702 tumours from 15 studies and evaluated on an independent set of 1067 tumours from six studies. Results of this study showed that expression levels of 276 genes were associated with OS (false discovery rate $<0.05$) in covariate-adjusted single-gene analyses.

We applied our method CausNet to this gene expression dataset of $513$ genes and survival outcome. For dimensionality reduction and parent set identification, we used three-level phenotype driven search. For the disease node ‘Status’, we carried out Cox propotional hazard regression for each gene separately, adjusted for age, stage, and stratified site. Then we computed analysis of variance tables for the fitted model objects and created a list of $p$-values based on $\chi^{2}$ distribution. We then adjusted these values using the Benjamini & Hochberg (BH) method. Then we choose $5$ genes - “ZFHX4” , “TIMP3”, “COL5A2”, “FBN1”, and “COL3A1” with the most significant p-values as possible parents of the disease node. At the next level, we used correlation between these and rest of the genes to pick their possible parent sets. After that, once again, parent sets were identified for these possible grand-parents of the disease node. These three levels of ancestors of the disease node picked $16$ genes in all with non-null parent sets. This reduced dataset and list of possible parents was used by CausNet with the BIC score to find a best network. Using the in-degree of $2$, we identified a best network as shown in Figure 10.

We implemented a dynamic programming based Optimal Bayesian Network (BN) Structure Discovery algorithm with parent set identification with ‘generational orderings’ based search for optimal networks, which is a novel way to efficiently search the space of possible networks given the possible parent set.

Our main novel contribution aside from providing software is the revision of the SM algorithm 3 [5] to incorporate possible parent sets and ‘generational orderings’ based search for a more efficient way to explore the search space as compared to the original approach based on lexicographical ordering. In doing so, we cover the entire constrained search space without searching through networks that don’t conform to the parent set constraints. While the basic algorithm can be applied to any dataset from any domain, the phenotype based algorithm is particularly suitable for disease outcome modeling.

The simulation results show that our algorithm performs very well when compared with three state-of-art algorithms that are widely used currently. The parent set constraints reduce the search space , and the runtime significantly, still giving better results than these algorithms, especially for number of variables greater than $60$.

The application to the recently published ovarian cancer gene-expression data with survival outcomes showed our algorithm’s usefulness in disease modeling. It yielded a sparse network of $6$ nodes leading to the disease outcome from gene expression data of $513$ genes in just a few minutes on a basic personal computer. This disease pathway can be used to gain useful insights for designing further medical and biological experiments.

Other features of the algorithm include specifiable parameters - correlation, FDR cutoffs, and in-degree - which can be tuned according to the application domain. Availability of two scoring option - BIC and Bge - and implementation of survival outcomes and mixed data types makes our algorithm very suitable for many types of biomedical data, e.g. GWAS and omics data to find disease pathways.

Funded by a Program Project Grant from the National Cancer Institute, P01 CA58959.

The authors would like to thank Professor Duncan C. Thomas for helpful inputs during the development of this work.

The CausNet software package in R is available at https://github.com/nand1155/CausNet.