Bayesian network in R: Introduction

[This article was first published on Ensemble Blogging, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

Bayesian networks (BNs) are a type of graphical model that encode the conditional probability between different learning variables in a directed acyclic graph. There are benefits to using BNs compared to other unsupervised machine learning techniques. A few of these benefits are:
  1. It is easy to exploit expert knowledge in BN models. 
  2. BN models have been found to be very robust in the sense of i) noisy data, ii) missing data and iii) sparse data. 
  3. Unlike many machine learning models (including Artificial Neural Network), which usually appear as a “black box,” all the parameters in BNs have an understandable semantic interpretation. 
This post is the first in a series of “Bayesian networks in R .” The goal is to study BNs and different available algorithms for building and training, to query a BN and examine how we can use those algorithms in R programming.
The R famous package for BNs is called “bnlearn”. This package contains different algorithms for BN structure learning, parameter learning and inference. In this introduction, we use one of the existing datasets in the package and show how to build a BN, train it and make an inference.

First let’s load the “coronary” dataset.
This data contains the following information:
  • Smoking (smoking): a two-level factor with levels no and yes.
  • M. Work (strenuous mental work): a two-level factor with levels no and yes.
  • P. Work (strenuous physical work): a two-level factor with levels no and yes.
  • Pressure (systolic blood pressure): a two-level factor with levels <140 and >140.
  • Proteins (ratio of beta and alpha lipoproteins): a two-level factor with levels.
  • Family (family anamnesis of coronary heart disease): a two-level factor with levels neg and pos.

Learn structure

The first step in a BN is to create the network. There are couples of algorithms in deriving an optimal BN structure and some of them exist in “bnlearn”. However, for the purpose of this post, we limit ourselves to a “max-min hill climbing” algorithm, which is a greedy algorithm [1].
bn_df <- data.frame(coronary)
res <- hc(bn_df)

The above structure finding creates the following conditional dependency between different variables, and the plot function draws the BN as shown below:

The causality between some nodes is intuitive; however, some relations extracted from data does not seem to be correct. For example, it does not make sense to have Family as a variable condition on M.Work. Therefore, we need to modify the derived structure. Let’s remove the link between M.Work and Family.
res$arcs <- res$arcs[-which((res$arcs[,'from'] == "M..Work" & res$arcs[,'to'] == "Family")),]


After learning the structure, we need to find out the conditional probability tables (CPTs) at each node. The function runs the EM algorithm to learn CPT for different nodes in the above graph.
fittedbn <-, data = bn_df)
For example, let look at what is inside the Protein node.
Protein is conditioned on M.Work and Smoking. Since both of these variables are binary variables (only two values) the CPT table has 2x2=4 entries:


Now, the BN is ready and we can start inferring from the network.
cpquery(fittedbn, event = (Proteins=="
which results in 0.61. Note that although the Proteins variable is conditioned on 2 variables, we did the query based on the available evidence on only one variables. But let make our evidence richer by asking the following: What is the chance that a non-smoker with pressure greater than 140 has a Proteins level less than 3?
cpquery(fittedbn, event = (Proteins=="140" ) )
which results in probability 0.63. 

We can also move in the opposite direction of an arc between two nodes. Let’s see if a person’s Proteins level is greater than 3, then what is the chance that his or her Pressure level is greater than 140?
cpquery(fittedbn, event = (Pressure==">140"), evidence = ( Proteins=="
the answer would be Pressure is greater than 140 with probability 0.41

[1] Tsamardinos, Ioannis, Laura E. Brown, and Constantin F. Aliferis. "The max-min hill-climbing Bayesian network structure learning algorithm." Machine learning 65.1 (2006): 31-78.

To leave a comment for the author, please follow the link and comment on their blog: Ensemble Blogging. offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

Never miss an update!
Subscribe to R-bloggers to receive
e-mails with the latest R posts.
(You will not see this message again.)

Click here to close (This popup will not appear again)