Learning Bayesian Networks with the bnlearn R Package

bnlearn is an R package which includes several algorithms for learning the structure of Bayesian networks with either discrete or continuous variables. Both constraint-based and score-based algorithms are implemented, and can use the functionality provided by the snow package to improve their performance via parallel computing. Several network scores and conditional independence algorithms are available for both the learning algorithms and independent use. Advanced plotting options are provided by the Rgraphviz package.


Introduction
In recent years Bayesian networks have been used in many fields, from On-line Analytical Processing (OLAP) performance enhancement (Margaritis 2003) to medical service performance analysis (Acid et al. 2004), gene expression analysis (Friedman et al. 2000), breast cancer prognosis and epidemiology (Holmes and Jain 2008).
The high dimensionality of the data sets common in these domains have led to the development of several learning algorithms focused on reducing computational complexity while still learning the correct network. Some examples are the Grow-Shrink algorithm in Margaritis (2003), the Incremental Association algorithm and its derivatives in Tsamardinos et al. (2003) and in Yaramakala and Margaritis (2005), the Sparse Candidate algorithm in Friedman et al. (1999), the Optimal Reinsertion in Moore and Wong (2003) and the Greedy Equivalent Search in Chickering (2002).
The aim of the bnlearn package is to provide a free implementation of some of these structure learning algorithms along with the conditional independence tests and network scores used to construct the Bayesian network. Both discrete and continuous data are supported. Furthermore, the learning algorithms can be chosen separately from the statistical criterion they are based on (which is usually not possible in the reference implementation provided by the algorithms' authors), so that the best combination for the data at hand can be used.

Bayesian networks
Bayesian networks are graphical models where nodes represent random variables (the two terms are used interchangeably in this article) and arrows represent probabilistic dependencies between them (Korb and Nicholson 2004).
The graphical structure G = (V, A) of a Bayesian network is a directed acyclic graph (DAG), where V is the node (or vertex ) set and A is the arc (or edge) set. The DAG defines a factorization of the joint probability distribution of V = {X 1 , X 2 , . . . , X v }, often called the global probability distribution, into a set of local probability distributions, one for each variable. The form of the factorization is given by the Markov property of Bayesian networks (Korb and Nicholson 2004, section 2.2.4), which states that every random variable X i directly depends only on its parents Π X i : The correspondence between conditional independence (of the random variables) and graphical separation (of the corresponding nodes of the graph) has been extended to an arbitrary triplet of disjoint subsets of V by Pearl (1988) with the d-separation (from direction-dependent separation). Therefore model selection algorithms first try to learn the graphical structure of the Bayesian network (hence the name of structure learning algorithms) and then estimate the parameters of the local distribution functions conditional on the learned structure. This two-step approach has the advantage that it considers one local distribution function at a time, and it does not require to model the global distribution function explicitly. Another advantage is that learning algorithms are able to scale to fit high-dimensional models without incurring in the so-called curse of dimensionality.
Although there are many possible choices for both the global and the local distribution functions, literature have focused mostly on two cases: • multinomial data (the discrete case): both the global and the local distributions are multinomial, and are represented as probability or contingency tables. This is by far the most common assumption, and the corresponding Bayesian networks are usually referred to as discrete Bayesian networks (or simply as Bayesian networks).
• multivariate normal data (the continuous case): the global distribution is multivariate normal, and the local distributions are normal random variables linked by linear constraints. These Bayesian networks are called Gaussian Bayesian networks in Geiger and Heckerman (1994), Neapolitan (2003) and most recent literature on the subject.
Other distributional assumptions lead to more complex learning algorithms (such as the nonparametric approach proposed by Bach and Jordan (2003)) or present various limitations due to the difficulty of specifying the distribution functions in closed form (such as the approach to learn Bayesian network with mixed variables by Boettcher and Dethlefsen (2003), which does not allow a node associated with a continuous variable to be the parent of a node associated with a discrete variable).

Structure learning algorithms
Bayesian network structure learning algorithms can be grouped in two categories: • constraint-based algorithms: these algorithms learn the network structure by analyzing the probabilistic relations entailed by the Markov property of Bayesian networks with conditional independence tests and then constructing a graph which satisfies the corresponding d-separation statements. The resulting models are often interpreted as causal models even when learned from observational data (Pearl 1988).
• score-based algorithms: these algorithms assign a score to each candidate Bayesian network and try to maximize it with some heuristic search algorithm. Greedy search algorithms (such as hill-climbing or tabu search) are a common choice, but almost any kind of search procedure can be used.
Constraint-based algorithms are all based on the Inductive Causation (IC) algorithm by Verma and Pearl (1991), which provides a theoretical framework for learning the structure causal models. It can be summarized in three steps: 1. first the skeleton of the network (the undirected graph underlying the network structure) is learned. Since an exhaustive search is computationally unfeasible for all but the most simple data sets, all learning algorithms use some kind of optimization such as restricting the search to the Markov blanket of each node (which includes the parents, the children and all the nodes that share a child with that particular node).
2. set all direction of the arcs that are part of a v-structure (a triplet of nodes incident on a converging connection X j → X i ← X k ).
3. set the directions of the other arcs as needed to satisfy the acyclicity constraint.
Score-based algorithms on the other hand are simply applications of various general purpose heuristic search algorithms, such as hill-climbing, tabu search, simulated annealing and various genetic algorithms. The score function is usually score-equivalent (Chickering 1995), so that networks that define the same probability distribution are assigned the same score.

Package implementation
4.1. Structure learning algorithms bnlearn implements the following constraint-based learning algorithms (the respective function names are reported in parenthesis): • Grow-Shrink (gs): based on the Grow-Shrink Markov Blanket, the simplest Markov blanket detection algorithm (Margaritis 2003) used in a structure learning algorithm.
• Incremental Association (iamb): based on the Incremental Association Markov blanket (IAMB) algorithm (Tsamardinos et al. 2003), which is based on a two-phase selection scheme (a forward selection followed by an attempt to remove false positives).
• Fast Incremental Association (fast.iamb): a variant of IAMB which uses speculative stepwise forward selection to reduce the number of conditional independence tests (Yaramakala and Margaritis 2005).
• Interleaved Incremental Association (inter.iamb): another variant of IAMB which uses forward stepwise selection (Tsamardinos et al. 2003) to avoid false positives in the Markov blanket detection phase.
• Max-Min Parents and Children (mmpc): a forward selection technique for neighbourhood detection based on the maximization of the minimum association measure observed with any subset of the nodes selected in the previous iterations (Tsamardinos et al. 2006). It learns the underlying structure of the Bayesian network (all the arcs are undirected, no attempt is made to detect their orientation).
Three implementations are provided for each algorithm: • an optimized implementation (used by default) which uses backtracking to roughly halve the number of independence tests.
• an unoptimized implementation (used when the optimized argument is set to FALSE) which is faithful to the original description of the algorithm. This implementation is particularly useful for comparing the behaviour of different combinations of learning algorithms and statistical tests.
• a parallel implementation. It requires a running cluster set up with the makeCluster function from the snow package (Tierney et al. 2008), which is passed to the function via the cluster argument.
The only available score-based learning algorithm is a Hill-Climbing (hc) greedy search on the space of directed graphs. The optimized implementation (again used by default) uses score caching, score decomposability and score equivalence to reduce the number of duplicated tests (Daly and Shen 2007). Random restarts, a configurable number of perturbing operations and a preseeded initial network structure can be used to avoid poor local maxima (with the restart, perturb and start arguments, respectively).

Conditional independence tests
Several conditional independence tests from information theory and classical statistics are available for use in constraint-based learning algorithms and the ci.test function. In both cases the test to be used is specified with the test argument (the label associated with each test is reported in parenthesis).
Conditional independence tests for discrete data are functions of the conditional probability tables implied by the graphical structure of the network through the observed frequencies {n ijk , i = 1, . . . , R, j = 1, . . . , C, k = 1, . . . , L} for the random variables X and Y and all the configurations of the conditioning variables Z: • mutual information: an information-theoretic distance measure (Kullback 1959), defined as It is proportional to the log-likelihood ratio test G 2 (they differ by a 2n factor, where n is the sample size) and it is related to the deviance of the tested models. Both the asymptotic χ 2 test (mi) and the Monte Carlo permutation test (mc-mi) described in Good (2005) are available.
• Pearson's X 2 : the classical Pearson's X 2 test for contingency tables, Again both the asymptotic χ 2 test (x2) and a Monte Carlo permutation test (mc-x2) from Good (2005) are available.
• fast mutual information (fmi): a variant of the mutual information which is set to zero when there aren't at least five data per parameter, which is the usual threshold for establishing the correctness of the asymptotic χ 2 distribution. This is the same heuristic defined for the Fast-IAMB algorithm in Yaramakala and Margaritis (2005).
• Akaike Information Criterion (aict): an experimental AIC-based independence test, computed comparing the mutual information and the expected information gain. It rejects the null hypothesis if which corresponds to an increase in the AIC score of the network.
In the continuous case conditional independence tests are functions of the partial correlation coefficients ρ XY | Z of X and Y given Z: • linear correlation: the linear correlation coefficient ρ XY | Z . Both the asymptotic Student's t test (cor) and the Monte Carlo permutation test (mc-cor) described in Legendre (2000) are available.
• Fisher's Z : a transformation of the linear correlation coefficient used by commercial software (such as TETRAD) and the pcalg package (Kalisch and Bühlmann 2007), which implements the PC constraint-based learning algorithm (Spirtes et al. 2001). It is defined as Both the asymptotic normal test (zf) and the Monte Carlo permutation test (mc-zf) are available.
• mutual information (mi-g): an information-theoretic distance measure (Kullback 1959), defined as It has the same relationship with the log-likelihood ratio as the corresponding test defined in the discrete case.

Network scores
Several score functions are available for use in the hill-climbing algorithm and the score function. The score to be used is specified with the score argument in hc and with the type argument in the score function (the label associated with each score is reported in parenthesis).
In the discrete case the following score functions are implemented: • the likelihood (lik) and log-likelihood (loglik) scores, which are equivalent to the entropy measure used by Weka (Witten and Frank 2005).
• the Akaike (aic) and Bayesian (bic) Information Criterion scores, defined as The latter is equivalent to the Minimum Description Length described by Rissanen (1978) and used as a Bayesian network score in Lam and Bacchus (1994).
• the logarithm of the Bayesian Dirichlet equivalent score (bde), a score equivalent Dirichlet posterior density (Heckerman et al. 1995).
• the logarithm of the K2 score (k2), another Dirichlet posterior density (Cooper and Herskovits 1992) and originally used in the structure learning algorithm of the same name. Unlike the bde score k2 is not score equivalent.
The only score available for the continuous case is a score equivalent Gaussian posterior density (bge), which follows a Wishart distribution (Geiger and Heckerman 1994).

Arc whitelisting and blacklisting
Prior information on the data, such as the ones elicited from experts in the relevant fields, can be integrated in all learning algorithms by means of the blacklist and whitelist arguments. Both of them accept a set of arcs which is guaranteed to be either present (for the former) or missing (for the latter) from the Bayesian network; any arc whitelisted and blacklisted at the same time is assumed to be whitelisted, and is thus removed from the blacklist.
This combination represents a very flexible way to describe any arbitrary set of assumptions on the data, and is also able to deal with partially directed graphs: • any arc whitelisted in both directions (i.e. both A → B and B → A are whitelisted) is present in the graph, but the choice of its direction is left to the learning algorithm. Therefore one of A → B, B → A and A−B is guaranteed to be in the Bayesian network.
• any arc blacklisted in both directions, as well as the corresponding undirected arc, is never present in the graph. Therefore if both A → B and B → A are blacklisted, also A − B is considered blacklisted.
• any arc whitelisted in one of its possible directions (i.e. A → B is whitelisted, but B → A is not) is guaranteed to be present in the graph in the specified direction. This effectively amounts to blacklisting both the corresponding undirected arc (A − B) and its reverse (B → A).
• any arc blacklisted in one of its possible directions (i.e. A → B is blacklisted, but B → A is not) is never present in the graph. The same holds for A − B, but not for B → A.

A simple example
In this section bnlearn will be used to analyze a small data set, learning.test. It's included in the package itself along with other real word and synthetic data sets, and is used in the example sections throughout the manual pages due to its simple structure.

Loading the package
bnlearn and its dependencies (the utils package, which is bundled with R) are available from CRAN, as are the suggested packages snow and graph

Learning a Bayesian network from data
Once bnlearn is loaded, learning.test itself can be loaded into a data frame of the same name with a call to data.

Constraint−based algorithms
The network structure learned by gs, iamb, fast.iamb and inter.iamb is equivalent to the one learned by hc; changing the arc A − B to either A → B or to B → A results in networks with the same score because of the score equivalence property (which holds for all the implemented score functions with the exception of K2). Therefore if there is any prior information about the relationship between A and B the appropriate direction can be whitelisted (or its reverse can blacklisted, which is equivalent in this case).

Network analysis and manipulation
The structure of a Bayesian network is uniquely specified if its graph is completely directed. In that case it can be represented as a string with the modelstring function

> modelstring(bn.hc) [1] "[A][C][F][B|A][D|A:C][E|B:F]"
whose output is also included in the print method for the objects of class bn. Each node is printed in square brackets along with all its parents (which are reported after a pipe as a colon-separated list), and its position in the string depends on the partial ordering defined by the network structure. The same syntax is used in deal (Boettcher and Dethlefsen 2003), an R package for learning Bayesian networks from mixed data.
Partially directed graphs can be transformed into completely directed ones with the set.arc, drop.arc and reverse.arc functions. For example the direction of the arc A − B in the bn.gs object can be set to A → B, so that the resulting network structure is identical to the one learned by the hill-climbing algorithm.

hc) [1] TRUE
Acyclicity is always preserved, as these commands return an error if the requested changes would result in a cyclic graph.
Further information on the network structure can be extracted from any bn object with the following functions: • whether the network structure is acyclic (acyclic) or completely directed (directed); • the labels of the nodes (nodes), of the root nodes (root.nodes) and of the leaf nodes (leaf.nodes); • the directed arcs (directed.arcs) of the network, the undirected ones (undirected.arcs) or both of them (arcs); • the adjacency matrix (amat) and the number of parameters (nparams) associated with the network structure; • the parents (parents), children (children), Markov blanket (mb), and neighbourhood (nbr) of each node.

Debugging utilities and diagnostics
Many of the functions of the bnlearn package are able to print additional diagnostic messages if called with the debug argument set to TRUE. This is especially useful to study the behaviour of the learning algorithms in specific settings and to investigate anomalies in their results (which may be due to an insufficient sample size for the asymptotic distribution of the tests to be valid, for example). For example the debugging output of the call to gs previously used to produce the bn.gs object reports the exact sequence of conditional independence tests performed by the learning algorithm, along with the effects of the backtracking optimizations (some parts are omitted for brevity). Other functions which provide useful diagnostics include (but are not limited to) compare (which reports the differences between the two network structures with respect to arcs, parents and children for each node), score and nparams (which report the number of parameters and the contribution of each node to the network score, respectively).

Practical examples
6.1. The ALARM network The ALARM ("A Logical Alarm Reduction Mechanism") network from Beinlich et al. (1989) is a Bayesian network designed to provide an alarm message system for patient monitoring. It has been widely used (see for example Tsamardinos et al. (2006) and Friedman et al. (1999)) as a case study to evaluate the performance of new structure learning algorithms.
The alarm data set includes a sample of size 20000 generated from this network, which contains 37 discrete variables (with two to four levels each) and 46 arcs. Every learning algorithm implemented in bnlearn (except mmpc) is capable of recovering the ALARM network to within a few arcs and arc directions (see Figure 2). > alarm.gs <-gs(alarm) > alarm.iamb <-iamb(alarm) > alarm.fast.iamb <-fast.iamb(alarm) > alarm.inter.iamb <-inter.iamb(alarm) > alarm.mmpc <-mmpc(alarm) > alarm.hc <-hc(alarm, score = "bic") The number of conditional independence tests, which provides an implementation independent performance indicator, is similar in all constraint-based algorithms (see Table 1); the same holds for the number of network score comparisons performed by hc, even though the learned network has about ten more arcs than the other ones.
The quality of the learned network improves significantly if the Monte Carlo versions of the tests are used instead of the parametric ones, as probability structure of the ALARM network results in many sparse contingency tables. For example a side by side comparison of the two versions of Pearson's X 2 test shows that the use of the nonparametric tests leads to the correct identification of all but five arcs, instead of the 12 missed with the parametric tests.  > alarm.mc <-gs(alarm, test = "mc-x2", B = 10000) > par(mfrow = c(1,2), omi = rep(0, 4), mar = c(1, 0, 1, 0)) > graphviz.plot(dag, highlight = list(arcs = arcs(alarm.gs))) > graphviz.plot(dag, highlight = list(arcs = arcs(alarm.mc)))

The examination marks data set
The marks data set is a small data set studied in Mardia et al. (1979), Whittaker (1990) and Edwards (2000). It contains five continuous variables, the examination marks for 88 students in five different subjects (mechanics, vectors, algebra, analysis and statistics). The purpose of the analysis was to find a suitable way to combine or average these marks. Since they are obviously correlated, the exact weights they are assigned depend on the estimated dependence structure of the data.

Pearson's Linear Correlation
All learning algorithms result in very similar network structures, which agree up to the direction of the arcs (see Figure 4). In all models the marks for analysis and statistics are conditionally independent from the ones for mechanics and vectors, given algebra. The structure of the graph suggests that the latter is essential in the overall evaluation of the examination.

Other packages for learning Bayesian networks
There exist other packages in R which are able to either learn the structure of a Bayesian network or fit and manipulate its parameters. Some examples are pcalg, which implements the PC algorithm and focuses on the causal interpretation of Bayesian networks; deal, which implements a hill-climbing search for mixed data; and the suite composed by gRbase , gRain (Højsgaard 2010), gRc (Højsgaard and Lauritzen 2008), which implements various exact and approximate inference procedures.
However, none of these packages is as versatile as bnlearn for learning the structure of Bayesian networks. deal and pcalg implement a single learning algorithm, even though are able to handle both discrete and continuous data. Furthermore, the PC algorithm has a poor performance in terms of speed and accuracy compared to newer constraint-based algorithms such as Grow-Shrink and IAMB (Tsamardinos et al. 2003). bnlearn also offers a wider selection of network scores and conditional independence tests; in particular it's the only R package able to learn the structure of Bayesian networks using permutation tests, which are superior to the corresponding asymptotic tests at low sample sizes.

Conclusions
bnlearn is an R package which provides a free implementation of some Bayesian network structure learning algorithms appeared in recent literature, enhanced with algorithmic optimizations and support for parallel computing. Many score functions and conditional independence tests are provided for both independent use and the learning algorithms themselves.
bnlearn is designed to provide the versatility needed to handle experimental data analysis. It handles both discrete and continuous data, and it supports any combination of the implemented learning algorithms and either network scores (for score-based algorithms) or conditional independence tests (for constraints-based algorithms). Furthermore, it simplifies the analysis of the learned networks by providing a single object class (bn) for all the algorithms and a set of utility functions to perform descriptive statistics and basic inference procedures.