BDgraph: An R Package for Bayesian Structure Learning in Graphical Models

Graphical models provide powerful tools to uncover complicated patterns in multivariate data and are commonly used in Bayesian statistics and machine learning. In this paper, we introduce an R package BDgraph which performs Bayesian structure learning for general undirected graphical models with either continuous or discrete variables. The package efficiently implements recent improvements in the Bayesian literature. To speed up computations, the computationally intensive tasks have been implemented in C++ and interfaced with R. In addition, the package contains several functions for simulation and visualization, as well as two multivariate datasets taken from the literature and are used to describe the package capabilities. The paper includes a brief overview of the statistical methods which have been implemented in the package. The main body of the paper explains how to use the package. Furthermore, we illustrate the package's functionality in both real and artificial examples, as well as in an extensive simulation study.


Introduction
Graphical models (Lauritzen, 1996) can be used to describe the conditional independence relationships among large numbers of variables. In graphical models, each random variable is associated with a node in a graph and links represent conditional dependency between variables, whereas the absence of a link implies that the variables are independent conditional on the rest of the variables (called the pairwise Markov property).
In recent years, significant progress has been made to design efficient algorithms to discover graph structures from high-dimensional multivariate data (Dobra et al., 2011, Jones et al., 2005, Mohammadi and Wit, 2015b, Mohammadi et al., 2015, Friedman et al., 2008, Meinshausen and Buhlmann, 2006, Murray and Ghahramani, 2004. In this regard, Bayesian approaches provide a principled alternative to various penalized approaches. In this paper, we describe the BDgraph package (Mohammadi and Wit, 2015a) in R (R Core Team, 2013) for Bayesian structure learning in undirected graphical models. The package can deal with Gaussian, non-Gaussian, discrete and mixed data sets. The package includes several functional modules, including data generation for simulation, a search algorithm, a graph estimation routine, a convergence check and a visualization tool; see figure 1.
The package efficiently implements recent improvements in the Bayesian literature, including Mohammadi and Wit (2015b), Mohammadi et al. (2015), Lenkoski (2013), Uhler et al. (2014), Hoff (2007). For a Bayesian framework of Gaussian graphical models, we implement the method developed by Mohammadi and Wit (2015b) and for Gaussian copula graphical models we use the method described by Mohammadi et al. (2015). To make our Bayesian methods computationally feasible for high-dimensional data, we efficiently implement the BDgraph package in C++ linked to R. To make the package easy to use, we use the S3 class. The package is available under the general public license (GPL ≥ 3) from the Comprehensive R Archive Network (CRAN) at http://cran.r-project.org/packages=BDgraph.
In the Bayesian literature, the BDgraph is the only software which is available online for Gaussian graphical models and Gaussian copula graphical models. On the other hand, in frequentest literature, existing packages are huge (Zhao et al., 2014), glasso (Friedman et al., 2014), bnlearn (Scutari, 2010), pcalg (Kalisch et al., 2012) and gRain (Højsgaard, 2012). We compare the performance to several packages.
The article is organized as follows. In Section 2 we illustrate the user interface of BDgraph package. In Section 3 we explain some methodological background of the package. In this regard, in Section 3.1, we briefly explain the Bayesian framework for Gaussian graphical models for continuous data. In Section 3.2 we briefly describe the Bayesian framework in the Gaussian copula graphical models for the data that not follow the Gaussianity assumption, such as non-Gaussian continuous, discrete or mixed data. In Section 4 we describe the main functions implemented in the BDgraph package. In Section 5 by a simple simulation examples, we explain the user interface and the performance of the package and we compare it with huge package. In Section 6, by using the functions implemented in the BDgraph package, we study two real data sets.

User interface
In R environment, we can access and load the BDgraph package by suing the following commands R> install.packages( "BDgraph" ) R> library( "BDgraph" ) Loading the BDgraph package automatically loads igraph (Csardi and Nepusz, 2006) and Matrix (Bates and Maechler, 2014) packages, since BDgraph package depends on these two packages. These packages are available on the Comprehensive R Archive Network (CRAN) at http://CRAN.R-project. org. We use igraph package for graph visualization and Matrix package for memory-optimization using sparse matrix output. In our package, all the functions operate on or return Matrix objects.
To maximize computational speed, we efficiently implement the BDgraph package with C++ code linked to R. For C++ code, we used the highly optimized LAPACK (Anderson et al., 1999) and BLAS (Lawson et al., 1979) linear algebra libraries on systems that provide them. We significantly improve program speed by using these libraries.
We design the BDgraph package to provide a Bayesian framework for highdimensional undirected graph estimation for different types of data sets such as continuous, discrete or mixed data. The package facilitates a flexible pipeline for analysis by three functional modules; see Figure 1. These modules are as follows: Module 1. Data Simulation: Function bdgraph.sim simulates multivariate Gaussian, discrete and mixed data with different undirected graph structures, including "random", "cluster", "hub", "scale-free", "circle"and "fixed"graphs. Users can set the sparsity of the graph structure. Users can generate mixed data, including "count", "ordinal", "binary", "Gaussian"and "non-Gaussian"variables.
Module 2. Method: The function bdgraph provides two estimation methods with two different algorithms for sampling from the posterior distributions: (i) Graph estimation for normally distributed data which is based on the Gaussian graphical models (GGMs) by using the birth-death MCMC sampling algorithm, described in Mohammadi and Wit (2015b).
(ii) Graph estimation for non-Gaussian, discrete, and mixed data, which is based on Gaussian copula graphical models (GCGMs) by using the birth-death MCMC sampling algorithm, described in Mohammadi et al. (2015).
Module 3. Result: This module includes four types of functions as follow (1) Convergence check: The functions plotcoda and traceplot provide several visualization plots to monitor the convergence of the sampling algorithm.
(2) Graph selection: The functions phat, select, prob and summary.bdgraph, provide the posterior link inclusion probabilities, posterior probability of each graph, selecting the best graph, respectively.
(3) Comparison and goodness-of-fit: The functions compare and plotroc provide several comparison measures and ROC plot for model comparison.
(4) Visualization: The plotting functions plot.bdgraph and plot.simulate provide visualizations of the simulated data and estimated graphs.

Methodological background
In Section 3.1, we briefly explain the Gaussian graphical model for normal data. Then we explain the birth-death MCMC algorithm for sampling from the joint posterior distribution in the Gaussian graphical models; for more detail see Mohammadi and Wit (2015b). In Section 3.2, we briefly describe the Gaussian copula graphical model, which can deal with non-Gaussian, discrete or mixed data. Then we explain the Birth-death MCMC algorithm which is designed for the Gaussian copula graphical models; for more detail see Mohammadi et al. (2015).

Bayesian Gaussian graphical models
In graphical models, each variable is associated with a node and conditional dependence relationships among random variables are presented as a graph G = (V, E) in which V = {1, 2, ..., p} specifies a set of nodes and a set of existing links E ⊂ V × V (Lauritzen, 1996). Our focus here is on undirected graphs in where (i, j) ∈ E is equivalent with (j, i) ∈ E. The absence of a link between two nodes specifies the pairwise conditional independence of these two variables given the remaining variables, while a link between two variables determines their conditional dependence.
In the Gaussian graphical models (GGMs), we assume the observed data follow a multivariate Gaussian distribution N p (µ, K −1 ), here we assume µ = 0. Let Z = (Z (1) , ..., Z (n) ) T be the observed data of n independent samples, then the likelihood function is where U = Z Z.
In GGMs, the conditional independence is implied by the form of the precision matrix. Based on pairwise Markov property, variables i and j are conditionally independence given the remaining variables, if and only if K ij = 0. This property implies that the links in graph G = (V, E) correspond with the nonzero elements of precision matrix K, it means E = {(i, j)|K ij = 0}. Given graph G, the precision matrix K is constrained to the cone P G of symmetric positive definite matrices with elements K ij equal to zeros for all (i, j) / ∈ E. We consider the G-Wishart distribution W G (b, D) as a prior distribution for the precision matrix K with density where b > 2 is the degree of freedom, D is a symmetric positive definite matrix, is the normalizing constant with respect to the graph G and 1(x) evaluates to 1 if x holds and to 0 otherwise. The G-Wishart distribution is a well-known prior for the precision matrix, since it represents the conjugate prior for normally distributed data. When G is complete the G-Wishart distribution reduces to the standard Wishart distribution, hence, its normalizing constant has an explicit form (Muirhead, 1982). Also, for decomposable graphs, the normalizing constant has an explicit form (Roverato, 2002). But, for non-decomposable graphs, the normalizing constant does not have an explicit form.
Since the G-Wishart prior is a conjugate prior to the likelihood (1), the posterior distribution of K is

Direct sampler from G-Wishart
Several sampling methods from G-Wishart distribution have been proposed; to review existing methods see Wang and Li (2012). More recently, Lenkoski (2013) have developed an exact sampling algorithm to sample from G−Wishart distribution, which borrows an idea from Hastie et al. (2009). The algorithm is as follows. Let N i ⊂ V be the neighbors set of node i in G. Form Ω Ni and Σ Ni,i and solveβ * i = Ω −1 Ni Σ Ni,i

5:
Formβ i ∈ R p−1 by padding the elements ofβ * i to the appropriate locations and zeros in those locations not connected to i in G 6: Update Ω i,−i and Ω −i,i with Ω −i,−iβi 7: end for 8: until convergence 9: return K = Ω −1 In the BDgraph package, we use Algorithm 1 to sample from the posterior distribution of the precision matrix. Besides, we implement the algorithm for general purposes in our package as a function rgwish; see the R code below for an illustration.

BDMCMC algorithm for GGMs
Consider the joint posterior distribution of K and the graph G given by For the graph-prior, as a non-informative prior, we consider a uniform distribution over all graph space, pr(G) ∝ 1. For the prior distribution of K condition on graph G we consider a G-Wishart W G (b, D).
Here, we consider a computationally efficient birth-death MCMC sampler proposed by Mohammadi and Wit (2015b) for the Gaussian graphical models. The algorithm is based on a continuous time birth-death Markov process in which the algorithm explores over the graph space by adding or removing a link in a birth or death event, respectively.
In the birth-death process, at a particular pair of graph G = (V, E) with precision matrix K, each link dies independently of the rest as a Poisson process with death rate δ e (K). Since the links are independent, the overall death rate is δ(K) = e∈G δ e (K). Birth rates β e (K) for e / ∈ G are defined similarly. Thus the overall birth rate is β(K) = e / ∈G β e (K). Since the birth and death events are independent Poisson processes, the time between two successive events exponentially distributed with mean 1/(β(K) + δ(K)). The time between successive events can be considered as support for any particular instance of graph G. The probability of the death and death events are The birth and death rates of links occur in continuous time with the rates determined by the stationary distribution of the process. The algorithm is designed in such a way that the stationary distribution equals the target joint posterior distribution of the graph and the precision matrix (3).
Mohammadi and Wit (2015b, section 3) prove by considering the birth and death rates as ratios of joint posterior distributions, as below, the birth-death MCMC sampling algorithm converges to the target joint posterior distribution of the graph and the precision matrix, if in which G +e = (V, E ∪ {e}) and K +e ∈ P G +e and similarly G −e = (V, E \ {e}) and K −e ∈ P G −e . Based on the above rates, we determine our BDMCMC algorithm as below. We design our BDMCMC sampling algorithm in such a way that, we sample from (G, K) in each steep of jumping to the new state, e.g. {t 1 , t 2 , ...}) in Algorithm 2 Given a graph G = (V, E) with a precision matrix K for N iteration do 1. Sample from the graph. Based on birth and death process 1.1. Calculate the birth rates by (6) and β(K) = e∈ / ∈G β e (K) 1.2. Calculate the death rates by (7) and δ(K) = e∈G δ e (K) 1.3. Calculate the waiting time by W (K) = 1/(β(K) + δ(K)) 1.4. Simulate the type of jump (birth or death) by (4) and (5) 2. Sample from the precision matrix. By using Algorithm 1. end for Figure 2. For efficient inference, we compute the sample means based on the Rao-Blackwellized estimator (Cappé et al., 2003, subsection 2.5); see e.g. 12.
Note, our aim is to estimate the posterior distribution of the graphs based on the data, P r(G|data). By using the Rao-Blackwellized estimator the estimated posterior distribution of the graphs are proportion to the total waiting times of each graph; see Figure 2.

Gaussian copula graphical models
In practice, we encounter both discrete and continuous variables; Gaussian copula graphical modeling has been proposed to describe dependencies between such heterogeneous variables. Let Y (as an observed data) be a collection of continuous, binary, ordinal or count variables with the marginal distribution F j of Y j and F −1 j its pseudo inverse. Towards constructing a joint distribution of Y, we introduce a multivariate Gaussian latent variable as follows where Γ(K) is the correlation matrix for a given precision matrix K. The joint distribution of Y is given by where C(·) is the Gaussian copula given by If all variables are continuous then the margins are unique; thus zeros in K imply conditional independence, as in Gaussian graphical models (Hoff, 2007). For discrete variables, the margins are not unique but still well-defined (Nelsen, 2007).
In semiparametric copula estimation, the marginals are treated as nuisance parameters and estimated by the rescaled empirical distribution. The joint distribution in (9) is then parametrized only by the correlation matrix of the Gaussian copula. We are interested to infer the underlying graph structure of the observed variables Y implied by the continuous latent variables Z. Since Z are unobservable we follow the idea of Hoff (2007) to associate them to the observed data as follows.
Given the observed data Y from a sample of n observations, we constrain the latent samples Z to belong to the set Following Hoff (2007) we infer on the latent space by substituting the observed data Y with the event D and defined the likelihood as The only part of the observed data likelihood relevant for inference on K is P r(Z ∈ D(Y) | K, G). Thus, the likelihood function is given by (11) where P r(Z | K, G) is defined in (1).

BDMCMC algorithm for GCGMs
The joint posterior distribution of the graph G and precision matrix K for our Gaussian copula graphical model is as follow P r(K, G|Z ∈ D(Y)) ∝ P r(K, G)P r(Z ∈ D(Y)|K, G).
Sampling from this posterior distribution can be done by using the birth-death MCMC algorithm. Mohammadi et al. (2015) have developed and extended the Algorithm 3 Given a graph G = (V, E) with a precision matrix K for N iteration do 1. Sample the latent data. For each r ∈ V and j ∈ {1, ..., n}, we update the latent values from its full conditional distribution as follows truncated to the interval L j r (Z), U j r (Z) in (10). 2. Sample from the graph. Same as Step 1 in the Algorithm 2. 3. Sample from the precision matrix. By using Algorithm 1. end for birth-death MCMC algorithm for more general case of GCGMs. We summarize their algorithm as follows.
In each iteration of the Algorithm 3, first conditional on the observed data (Y) we sample from the latent variables (Z). The other steps are the same as the Algorithm 2.
Remark. For the cases that all variables are continuous, we do not need to sample from latent variables in each iteration of the Algorithm 2, since all margins in Gaussian copula are unique. Therefore, for these cases, first based on Gaussian copula approach, we transfer our non-Gaussian data to Gaussian then we run the Algorithm 2.

The BDgraph environment
The BDgraph package provides a set of comprehensive tools related to Bayesian graphical models; here, we describe the essential functions available in this package.

Description of the bdgraph function
The main function of the package is bdgraph, which includes two Bayesian frameworks (GGMs and GCGMs). This function is based on an underlying sampling engine (Algorithms 2 and 3) takes the model definition and returns a sequence of samples from the joint posterior distribution of the graph and precision matrix (3), given the supplied data. By using the following function bdgraph( data, n = NULL, method = "ggm", iter = 5000, burnin = iter / 2, b = 3, D = NULL, Gstart = "empty" ) we can sample from our target joint posterior distribution. bdgraph returns an object of S3 class type "bdgraph". The functions plot, print and summary are working with the object "bdgraph". The input data can be a matrix or a data.frame ( n × p) or a covariance matrix; it can also be an object of class "simulate", which is output of the function bdgraph.sim. The argument method determines the type of BDMCMC algorithm (Algorithms 2 and 3); option "ggm"implements the BDMCMC algorithm based on the Gaussian graphical models (Algorithm 2). It is designed for the data that follow Gaussianity assumption. Option "gcgm"implements the BDMCMC algorithm based on the Gaussian copula graphical models (Algorithm 3). It is designed for the data that do not follow the Gaussianity assumption such as non-Gaussian continuous, discrete or mixed data. This option can deal with missing data.

Description of the bdgraph.sim function
The function bdgraph.sim is designed to simulate different types of data sets with different graph structures. The function bdgraph.sim( n = 2, p = 10, graph = "random", size = NULL, prob = 0.2, class = NULL, type = "Gaussian", cut = 4, b = 3, D = diag(p), K = NULL, sigma = NULL, mean = 0, vis = FALSE ) can simulate multivariate Gaussian, non-Gaussian, discrete and mixed data with different undirected graph structures, including "random", "cluster", "hub", "scale-free", "circle"and "fixed"graphs. Users can determine the sparsity level of the obtained graph by option prob. By this function, users can generate mixed data from "count", "ordinal", "binary", "Gaussian"and "non-Gaussian"distributions. bdgraph.sim returns an object of the S3 class type "simulate". The plot, print and summary functions are working with this object type.

Description of the plotcoda and traceplot functions
In general, convergence in MCMC approaches can be difficult to evaluate. From a theoretical point of view, the sampling distribution will converge to the target joint posterior distribution as the number of iteration increases to infinity. We normally have little theoretical insight about how quickly convergence "KICKS IN"; therefore we rely on post hoc testing of the sampled output. In general, the sample is divided into two parts: a "burn-in"part of sample and the remainder, in which the chain is considered to have converged sufficiently close to the target posterior distribution. Two questions then arise: How many samples are sufficient? How long should the burn-in period be? The

Description of the phat and select functions
In the BDgraph package, the functions phat and select provide the potential tools to do statistical inference from the samples drawn from the joint posterior distribution. The function phat( output, round = 3 ) provides the estimated posterior link inclusion probabilities for all possible links. These probabilities, for all possible link e = (i, j) in graph G, can be calculated by using the Rao-Blackwellization (Cappé et al., 2003, subsection where N is the number of iteration, I(e ∈ G (t) ) is an indicator function, such that I(e ∈ G (t) ) = 1 if e ∈ G (t) and zero otherwise, and W (K (t) ) are the weights of the graph G (t) with the precision matrix K (t) . The function

select( output, cut = NULL, vis = FALSE )
provides the inferred graph, which is the graph with the highest posterior probability as a default. By option cut, users can select the inferred graph based on the estimated posterior link inclusion probabilities.

Description of the compare and plotroc functions
The function compare and plotroc are designed to check and compare the performance of our inferred graph. This function is particularly useful for simulation studies. With the function compare( G, est, est2 = NULL, colnames = NULL, vis = FALSE ) we can check the performance of our BDMCMC algorithm and compare it with other alternative approaches. In this regard, this function provides several measures such as balanced F -score measure (Baldi et al., 2000) which is defined as follows where TP, FP and FN are the number of true positives, false positives and false negatives, respectively. The F 1 -score lies between 0 and 1, where 1 is for perfect identity and 0 for worst case. The function plotroc( G, prob, prob2 = NULL, cut = 20, smooth = FALSE ) provides a ROC plot for visualization comparison based on the estimated posterior link inclusion probabilities (12).

User interface by toy example
In this section, we illustrate the user interface of the BDgraph package with a simple simulation example. We perform all the computations on an Intel(R) Core(TM) i5 CPU 2.67GHz processor. By using the function bdgraph.sim we simulate 60 observations (n = 60) from a multivariate Gaussian distribution from a "random"graph with 8 nodes (p = 8), as below R> data.sim <-bdgraph.sim( n = 60, p = 8, graph = "scale-free", type = "Gaussian" ) R> round( head( data.sim $ data, 4), 2 )

Running the BDMCMC algorithm
Since the generated data follow Gaussianity assumption, we run the BDMCMC algorithm which is based on the Gaussian graphical models. Therefore, we choose method = "ggm", as follow R> output.ggm <-bdgraph( data = data.sim, method = "ggm", iter = 50000 ) Running this function takes around 5 second which is computationally fast as main sampling part is in C++. Since the function bdgraph returns an object of class S3, users can see the summary result of BDMCMC algorithm as follows The summary result are the adjacency matrix of the graph with the highest posterior probability (selected graph), the estimated posterior probabilities of all possible links (phat) and the estimated precision matrix (Khat). In addition, function summary reports a visualization summary of the results as we can see in Figure 3. In the top-left is the graph with the highest posterior probability. In the top-right is the estimated posterior probabilities of all the graphs which are visited by the BDMCMC algorithm. It indicates our algorithm visits more than 600 different graphs and the posterior probability of the selected graph is around 0.22. In the bottom-left is the estimated posterior probabilities of the size of graphs. It indicates our algorithm visited mainly graphs with size between 6 and 14 links. In the bottom-right is the trace of our algorithm based on the size of the graphs.

Convergence check
In our simulation example, we run the BDMCMC algorithm for 50, 000 iterations, 25, 000 of which as burn-in. To check whether the number of iteration is enough or not, we run R> plotcoda( output.ggm ) The results in Figure 4 indicates that our BDMCMC algorithm, roughly speaking, converges after 20, 000 iterations and that burn-in of 25, 000 is sufficient.

Comparison and goodness-of-fit
The function compare provides several measures for checking the performance of our algorithms and compared with other alternative approaches, with respect to the true graph structure. To check the performance of our both algorithms (Algorithms 2 and 3), we also run the Algorithms 3 with the same condition as bellow R> output.gcgm <-bdgraph( data = data.sim, method = "gcgm", iter = 50000 ) where output.ggm is 50, 000 samples from the joint posterior distribution which are generated based on the Gaussian copula graphical models.
Users can compare the performance of these two algorithms by using the code R> plotroc( data.sim, output.ggm, output.gcgm ) As expected, the result indicates that the "GGM"approach performs slightly better than the "GCGM", since the data is generated from a Gaussian graphical model. For more comparisons of these two approaches see Mohammadi and Wit (2014).

Application to real data sets
In this section, we analyze two data sets from biology and sociology, by using the functions that are available in the BDgraph package. In Section 6.1, we analyze a Labor force survey data set, which are mixed data. In Section 6.2, we analyze Human gene expression data, which are high-dimensional data that do not follow the Gaussianity assumption. Both data sets are available in our package.

Application to Labor force survey data
Hoff (2007)  Since those variables are measured on various scales, the marginal distributions are heterogeneous, which makes the study of their joint distribution very challenging. However, we can apply our Bayesian framework based on the Gaussian copula graphical models to this data set. Thus, we run the function bdgraph with option method = "gcgm". For the prior distributions of the graph and precision matrix, as default of the function bdgraph, we place a uniform distribution as a uninformative prior on the graph and a G-Wishart W G (3, I 7 ) on the precision matrix. We run our function for 50, 000 iterations with 25, 000 sweeps as burn-in.
R> output <-bdgraph( data = surveyData, method = "gcgm", iter = 50000 ) R> summary( output ) $selected_graph income degree children pincome pdegree pchildren age The result of the function summary are the adjacency matrix of the graph with the highest posterior probability (selected graph), estimated posterior probabilities of all possible links (phat) and estimated precision matrix (Khat). Figure 7 shows the selected graph which is a graph with the highest posterior probability. Links in the graph are indicated by signs "+"and "-", which represents a positively and negatively correlated relationship between associated variables, respectively.
The result indicate that education, fertility and age determinate income, since there are highly positively correlated relationships between income and those three variables, with posterior probability equal to one for all. It shows that respondent's education and fertility are negatively correlated with posterior probability equal to 0.67. The respondent's education is certainly related to his parent's education, since there is positively correlated relationship with posterior probability equal to one. Moreover, the result indicate that relationships between income, education and fertility hold across generations.
For this data set, Hoff (2007) estimated the conditional independence between variables by inspecting whether the 95% credible intervals for the associated regression parameters do not contain zero. Our result is the same with Hoff's result except for one link. Our result indicate there is a strong relationship between parent's education (pdegree) and fertility (child), which is not selected by Hoff.  Figure 7: Graph with the highest posterior probability for the labor force survey data based on the output from bdgraph. Sign "+"represents a positively correlated relationship between associated variables and "-"represents a negatively correlated relationship.

Application to Human gene expression
Here, by using the functions that are available in the BDgraph package, we study the structure learning of the sparse graphs applied to the large-scale human gene expression data which was originally described by Stranger et al. (2007). They collected the data to measure gene expression in B-lymphocyte cells from Utah individuals of Northern and Western European ancestry. They consider 60 individuals whose genotypes are available online at ftp://ftp. sanger.ac.uk/pub/genevar. Here the focus is on the 3125 Single Nucleotide Polymorphisms (SNPs) that have been found in the 5' UTR (untranslated region) of mRNA (messenger RNA) with a minor allele frequency ≥ 0.1. Since the UTR (untranslated region) of mRNA (messenger RNA) has been subject to investigation previously, it should have an important role in the regulation of the gene expression. The raw data were background corrected and then quantile normalized across replicates of a single individual and then median normalized across all individuals. Following Bhadra and Mallick (2013), among the 47, 293 total available probes, we consider the 100 most variable probes that correspond to different Illumina TargetID transcripts. The data for these 100 probes are available in our package. To see the data, users can run the code R> data(geneExpression) R> dim( geneExpression ) [1] 60 100 The data consist of only 60 observations (n = 60) across 100 genes (p = 100). This data set is an interesting case study for graph structure learning, as it has been used in (Bhadra and Mallick, 2013, Mohammadi and Wit, 2015b, Gu et al., 2015. In this data set, all the variables are continuous but they do not follow the Gaussianity assumption, as you can see in Figure 8. Thus, we apply the Gaussian copula graphical models. Therefore, we run function bdgraph with option method = "gcgm". For the prior distributions of the graph and precision matrix, as default of the function bdgraph, we place a uniform distribution as a uninformative prior on graph and the G-Wishart W G (3, I 100 ) on the precision matrix. We run our function for 100, 000 iterations with 50, 000 sweeps as burn-in as follows R> output <-bdgraph( data = geneExpression, method = "ggm", iter = 100000 ) R> select( output ) This function takes around 13 hours. The function select gives as a default the graph with the highest posterior probability, which has 991 links. We use the following code to visualize the ones that we believe have posterior probabilities larger than 0.995.
R> select( output, cut = 0.995, vis = TRUE ) By using option vis = TRUE, this function plots the selected graph. Figure  9 shows the selected graph with 266 links, for which the posterior inclusion probabilities (12) are greater than 0.995.
The function phat reports the estimated posterior probabilities of all possible links in the graph. For our data the output of this function is a 100×100 matrix. Figure 10 reports the visualization of that matrix.
Most of the links in our selected graph (graph with the highest posterior probability) conform the results in previous studies. For instance, Bhadra and Mallick (2013) found 54 significant interactions between genes, most of which of this method would be desirable in real applications.