Fast Bayesian inference in large Gaussian graphical models

TL;DR

This paper introduces an efficient Bayesian method for inferring large Gaussian graphical models by using closed-form Bayes factors for hypothesis testing, enabling scalable analysis of high-dimensional data.

Contribution

It develops a computationally efficient approach with exact tail probabilities for Bayesian inference in large Gaussian graphical models, addressing high-dimensional challenges.

Findings

Method performs well on simulated data

Successfully applied to large gene expression dataset

Controls error rates effectively in edge selection

Abstract

Despite major methodological developments, Bayesian inference for Gaussian graphical models remains challenging in high dimension due to the tremendous size of the model space. This article proposes a method to infer the marginal and conditional independence structures between variables by multiple testing of hypotheses. Specifically, we introduce closed-form Bayes factors under the Gaussian conjugate model to evaluate the null hypotheses of marginal and conditional independence between variables. Their computation for all pairs of variables is shown to be extremely efficient, thereby allowing us to address large problems with thousands of nodes. Moreover, we derive exact tail probabilities from the null distributions of the Bayes factors. These allow the use of any multiplicity correction procedure to control error rates for incorrect edge inclusion. We demonstrate the proposed approach to graphical model selection on various simulated examples as well as on a large gene expression data set from The Cancer Genome Atlas.

Figures

Tables

Table 1. Table 1: Average and standard deviation (in parenthesis) of areas under the receiver operating characteristic and precision-recall curves over the simulated datasets, as a function of the true graph structure and sample size n 𝑛 n . beam , our method; bdmcmc and rjmcmc methods of Mohammadi and Wit, ( 2015 ) ; saturnin method of Schwaller et al., ( 2017 ) ; AUC ROC ROC {}_{\text{ROC}} , area under the receiver operating characteristic curve; AUC PR PR {}_{\text{PR}} area under the precision-recall curve. Best performances are boldfaced.

$n$	Method	AUC${}_{\text{ROC}}$	AUC${}_{\text{PR}}$	AUC${}_{\text{ROC}}$	AUC${}_{\text{PR}}$
		Band structure		Cluster structure
100	beam	0.89 (0.02)	0.65 (0.03)	0.80 (0.02)	0.54 (0.03)
100	bdmcmc	0.89 (0.03)	0.67 (0.03)	0.79 (0.02)	0.51 (0.04)
100	rjmcmc	0.88 (0.03)	0.63 (0.05)	0.78 (0.03)	0.50 (0.04)
100	saturnin	0.89 (0.02)	0.61 (0.04)	0.77 (0.02)	0.53 (0.04)
50	beam	0.84 (0.03)	0.53 (0.04)	0.73 (0.02)	0.39 (0.04)
50	bdmcmc	0.82 (0.03)	0.51 (0.06)	0.72 (0.03)	0.37 (0.04)
50	rjmcmc	0.81 (0.03)	0.47 (0.05)	0.72 (0.02)	0.35 (0.04)
50	saturnin	0.82 (0.02)	0.44 (0.04)	0.68 (0.02)	0.33 (0.04)
25	beam	0.78 (0.04)	0.39 (0.05)	0.66 (0.03)	0.24 (0.04)
25	bdmcmc	0.75 (0.04)	0.32 (0.05)	0.65 (0.03)	0.23 (0.03)
25	rjmcmc	0.75 (0.04)	0.27 (0.05)	0.64 (0.03)	0.22 (0.03)
25	saturnin	0.73 (0.03)	0.28 (0.05)	0.58 (0.02)	0.15 (0.02)
		Hub structure		Random structure
100	beam	0.88 (0.03)	0.62 (0.03)	0.87 (0.03)	0.65 (0.03)
100	bdmcmc	0.89 (0.02)	0.67 (0.04)	0.86 (0.03)	0.66 (0.03)
100	rjmcmc	0.89 (0.02)	0.65 (0.05)	0.85 (0.03)	0.65 (0.04)
100	saturnin	0.92 (0.01)	0.63 (0.02)	0.86 (0.02)	0.59 (0.02)
50	beam	0.84 (0.03)	0.53 (0.03)	0.83 (0.03)	0.56 (0.04)
50	bdmcmc	0.84 (0.03)	0.52 (0.05)	0.81 (0.03)	0.53 (0.05)
50	rjmcmc	0.84 (0.03)	0.48 (0.06)	0.80 (0.03)	0.49 (0.06)
50	saturnin	0.86 (0.02)	0.48 (0.03)	0.83 (0.02)	0.47 (0.03)
25	beam	0.80 (0.03)	0.42 (0.04)	0.79 (0.03)	0.43 (0.05)
25	bdmcmc	0.79 (0.04)	0.32 (0.05)	0.75 (0.02)	0.33 (0.05)
25	rjmcmc	0.77 (0.04)	0.27 (0.04)	0.74 (0.03)	0.30 (0.05)
25	saturnin	0.80 (0.03)	0.35 (0.04)	0.77 (0.02)	0.35 (0.04)

Table 2. Table 2: Average and standard deviation (in parenthesis) areas under the receiver operating characteristic and precision-recall curves over the simulated datasets, and as a function of the true graph structure and sample size n 𝑛 n . beam , our method; saturnin method of Schwaller et al., ( 2017 ) ; GeneNet method of Schäfer and Strimmer, ( 2005 ) ; FastGGM method of Ren et al., ( 2015 ) ; AUC ROC ROC {}_{\text{ROC}} , area under the receiver operating characteristic curve; AUC PR PR {}_{\text{PR}} area under the precision-recall curve.

$p$	Method	AUC${}_{\text{ROC}}$	AUC${}_{\text{PR}}$	AUC${}_{\text{ROC}}$	AUC${}_{\text{PR}}$
		Band structure		Cluster structure
200	beam	0.88 (0.01)	0.55 (0.02)	0.91 (0.01)	0.58 (0.01)
200	GeneNet	0.89 (0.01)	0.57 (0.02)	0.91 (0.01)	0.59 (0.01)
200	FastGGM	0.87 (0.01)	0.57 (0.02)	0.89 (0.01)	0.60 (0.02)
500	beam	0.91 (0.01)	0.58 (0.01)	0.89 (0.01)	0.50 (0.01)
500	GeneNet	0.91 (0.01)	0.60 (0.01)	0.89 (0.01)	0.52 (0.01)
500	FastGGM	0.90 (0.01)	0.61 (0.01)	0.85 (0.01)	0.49 (0.01)
1000	beam	0.88 (0.01)	0.49 (0.01)	0.90 (0.00)	0.48 (0.01)
1000	GeneNet	0.88 (0.01)	0.49 (0.01)	0.90 (0.00)	0.49 (0.01)
1000	FastGGM	0.87 (0.01)	0.51 (0.01)	0.87 (0.00)	0.48 (0.01)
		Hub structure		Random structure
200	beam	0.90 (0.01)	0.56 (0.01)	0.86 (0.01)	0.43 (0.02)
200	GeneNet	0.85 (0.01)	0.21 (0.03)	0.86 (0.01)	0.47 (0.02)
200	FastGGM	0.87 (0.01)	0.46 (0.02)	0.85 (0.01)	0.47 (0.02)
500	beam	0.92 (0.01)	0.54 (0.01)	0.82 (0.01)	0.35 (0.01)
500	GeneNet	0.90 (0.00)	0.43 (0.01)	0.82 (0.01)	0.34 (0.01)
500	FastGGM	0.88 (0.01)	0.44 (0.01)	0.81 (0.00)	0.34 (0.01)
1000	beam	0.93 (0.00)	0.54 (0.01)	0.77 (0.00)	0.22 (0.01)
1000	GeneNet	0.92 (0.00)	0.49 (0.01)	0.77 (0.00)	0.21 (0.01)
1000	FastGGM	0.89 (0.00)	0.44 (0.01)	0.77 (0.00)	0.22 (0.01)

Equations

$$(2{)}^{-\frac{\alpha d}{2}}{\Gamma }_d^{-1}\left(\frac{\alpha }{2}\right)|A{|}^{-\frac{\alpha }{2}}|\Omega {|}^{\frac{\alpha -d-1}{2}}\exp \left\{-\frac{1}{2}\text{tr}(A^{-1}\Omega )\right\},$$

$$\begin{array}{r}\text{vec}(Y)\mid \Sigma \sim N_{np}(0,\Sigma \otimes I_n)\text{and}\Sigma \mid D,\delta \sim I{W}_p\left(\left(\delta -p-1\right)D,\delta \right)\end{array},$$

$$E\left(\Sigma \mid Y\right)=\left\{\left(\delta -p-1\right)D+S\right\}/(\delta +n-p-1),$$

$$E\left(\Omega \mid Y\right)=(\delta +n){\left\{\left(\delta -p-1\right)D+S\right\}}^{-1}.$$

$$\hat{\delta }=\underset{\delta }{\mathrm{argmax}}p(Y;\delta ),$$

$$\begin{array}{r}p(Y;\delta )={\pi }^{-(np)/2}\frac{{\Gamma }_p\left(\frac{\delta +n}{2}\right)}{{\Gamma }_p\left(\frac{\delta }{2}\right)}\frac{\mid \left(\delta -p-1\right)D{\mid }^{\frac{\delta }{2}}}{\mid \left(\delta -p-1\right)D+S{\mid }^{\frac{\delta +n}{2}}}.\end{array}$$

$${\text{BF}}_{ij}^{\text{C}}=\frac{\int p_1(Y\mid \Sigma )p_1(\Sigma )d\Sigma }{\int p_0(Y\mid {\Sigma }^0)p_0({\Sigma }^0)d{\Sigma }^0},$$

$$p(Y|\Sigma )=p(Y_b|{\Sigma }_{bb})p(Y_a|Y_b,B_{a|b},{\Sigma }_{aa.b}),$$

$${\text{BF}}_{ij}^{\text{C}}=\frac{\iint p_1(Y_a\mid Y_b,B_{a\mid b},{\Sigma }_{aa.b})p_1(B_{a\mid b},{\Sigma }_{aa.b})d{B}_{a\mid b}d{\Sigma }_{aa.b}}{\iint p_0(Y_a\mid Y_b,B_{a\mid b},{\Sigma }_{aa.b}^0)p_0(B_{a\mid b},{\Sigma }_{aa.b}^0)d{B}_{a\mid b}d{\Sigma }_{aa.b}^0}.$$

$$\begin{array}{rl}\text{vec}(Y_a)\mid Y_b,B_{a\mid b},{\Sigma }_{aa.b}& \sim N_{n\times 2}\left(\text{vec}(Y_b{B}_{a\mid b}),{\Sigma }_{aa.b}\otimes I_n\right),\\ \text{vec}(B_{a\mid b})\mid {\Sigma }_{aa.b}& \sim N_{(p-2)\times 2}\left(\text{vec}(F_{a\mid b}),{\Sigma }_{aa.b}\otimes F_{bb}^{-1}\right),\\ {\Sigma }_{aa.b}& \sim I{W}_2\left(F_{aa.b},\delta \right),\end{array}$$

$$\begin{array}{rl}\text{vec}(Y_a)\mid Y_b,B_{a\mid b},{\Sigma }_{aa.b}^0& \sim N_{n\times 2}\left(\text{vec}(Y_b{B}_{a\mid b}),{\Sigma }_{aa.b}^0\otimes I_n\right),\\ \text{vec}(B_{a\mid b})\mid {\Sigma }_{aa.b}^0& \sim N_{(p-2)\times 2}\left(\text{vec}(F_{a\mid b}),{\Sigma }_{aa.b}^0\otimes F_{bb}^{-1}\right),\\ p_0(B_{a\mid b},{\Sigma }_{aa.b}^0)& =p_1(B_{a\mid b},{\Sigma }_{aa.b}\mid {\text{H}}_{0,ij}^{\text{C}})\\ & =\frac{p_1(B_{a\mid b},{\Sigma }_{aa.b},{\text{H}}_{0,ij}^{\text{C}})}{\iint p_1(B_{a\mid b},{\Sigma }_{aa.b},{\text{H}}_{0,ij}^{\text{C}})d{B}_{a\mid b}d{\Sigma }_{aa.b}},\end{array}$$

$${\text{BF}}_{ij}^{\text{C}}=\frac{\Gamma \left(\frac{\delta +n}{2}\right)\Gamma \left(\frac{\delta +n-1}{2}\right){\Gamma }^2\left(\frac{\delta +1}{2}\right)}{\Gamma \left(\frac{\delta }{2}\right)\Gamma \left(\frac{\delta -1}{2}\right){\Gamma }^2\left(\frac{\delta +n+1}{2}\right)}\frac{(1-r_{g_{ij}}^2{)}^{\frac{\delta }{2}}}{(1-r_{q_{ij}}^2{)}^{\frac{\delta +n}{2}}}{\left(\frac{g_{ii}{g}_{jj}}{q_{ii}{q}_{jj}}\right)}^{\frac{1}{2}},$$

$${\text{BF}}_{ij}^{\text{M}}=\frac{\int p_1(Y\mid \Sigma )p_1(\Sigma )d\Sigma }{\int p_0(Y\mid {\Sigma }^0)p_0({\Sigma }^0)d{\Sigma }^0},$$

$$p(Y\mid \Sigma )=p(Y_a\mid {\Sigma }_{aa})p(Y_b\mid Y_a,B_{b\mid a},{\Sigma }_{bb.a}),$$

$$\begin{array}{rl}{\text{BF}}_{ij}^{\text{M}}& =\frac{\int p_1(Y_a\mid {\Sigma }_{aa})p_1({\Sigma }_{aa})d{\Sigma }_{aa}}{\int p_0(Y_a\mid {\Sigma }_{aa}^0)p_0({\Sigma }_{aa}^0)d{\Sigma }_{aa}^0}.\end{array}$$

$$\begin{array}{rl}\text{vec}(Y_a)\mid {\Sigma }_{aa}& \sim N_{n\times 2}(0,{\Sigma }_{aa}\otimes I_n),\\ {\Sigma }_{aa}& \sim I{W}_2(F_{aa},\delta -p+2),\end{array}$$

$$\begin{array}{rl}\text{vec}(Y_a)\mid {\Sigma }_{aa}^0& \sim N_{n\times 2}(0,{\Sigma }_{aa}^0\otimes I_n),\\ p_0({\Sigma }_{aa}^0)& =p_1({\Sigma }_{aa}\mid {\text{H}}_{0,ij}^{\text{M}})=\frac{p_1({\Sigma }_{aa},{\text{H}}_{0,ij}^{\text{M}})}{\int p_1({\Sigma }_{aa},{\text{H}}_{0,ij}^{\text{M}})d{\Sigma }_{aa}}.\end{array}$$

$${\text{BF}}_{ij}^{\text{M}}=\frac{{\Gamma }_2\left(\frac{\delta +n-p+2}{2}\right){\Gamma }^2\left(\frac{\delta -p+3}{2}\right)}{{\Gamma }_2\left(\frac{\delta -p+2}{2}\right){\Gamma }^2\left(\frac{\delta +n-p+3}{2}\right)}\frac{{\left(1-r_{f_{ij}}^2\right)}^{\frac{\delta -p+2}{2}}}{{\left(1-r_{t_{ij}}^2\right)}^{\frac{\delta +n-p+2}{2}}}{\left(\frac{t_{ii}{t}_{jj}}{f_{ii}{f}_{jj}}\right)}^{\frac{1}{2}},$$

$${\text{sBF}}_{ij}^{\text{C}}=\frac{\Gamma \left(\frac{\delta +n}{2}\right)\Gamma \left(\frac{\delta +n-1}{2}\right){\Gamma }^2\left(\frac{\delta +1}{2}\right)}{\Gamma \left(\frac{\delta }{2}\right)\Gamma \left(\frac{\delta -1}{2}\right){\Gamma }^2\left(\frac{\delta +n+1}{2}\right)}\frac{(1-r_{g_{ij}}^2{)}^{\frac{\delta }{2}}}{(1-r_{q_{ij}}^2{)}^{\frac{\delta +n}{2}}},$$

$${\text{sBF}}_{ij}^{\text{M}}=\frac{{\Gamma }_2\left(\frac{\delta +n-p+2}{2}\right){\Gamma }^2\left(\frac{\delta -p+3}{2}\right)}{{\Gamma }_2\left(\frac{\delta -p+2}{2}\right){\Gamma }^2\left(\frac{\delta +n-p+3}{2}\right)}\frac{{\left(1-r_{f_{ij}}^2\right)}^{\frac{\delta -p+2}{2}}}{{\left(1-r_{t_{ij}}^2\right)}^{\frac{\delta +n-p+2}{2}}},$$

$$\Phi =\left(\begin{array}{cc}{\varphi }_1^2& {\varphi }_1{\varphi }_2\phi \\ {\varphi }_1{\varphi }_2\phi & {\varphi }_2^2\end{array}\right)\text{and}\Sigma =\left(\begin{array}{cc}{\sigma }_1^2& {\sigma }_1{\sigma }_2\rho \\ {\sigma }_1{\sigma }_2\rho & {\sigma }_2^2\end{array}\right)$$

$$({\phi }^2\mid \rho =0)\sim Beta(1/2,(d-1)/2).$$

$$\begin{array}{rl}{Y}_a^T{Y}_a-{\bar{B}}_{a\mid b}^T& (Y_b^T{Y}_b+F_{bb}){\bar{B}}_{a\mid b}+F_{ab}{F}_{bb}^{-1}{F}_{ba}=\\ & (Y_a-Y_b{F}_{a|b}{)}^T(I_n+Y_b{F}_{bb}^{-1}{Y}_b^T{)}^{-1}(Y_a-Y_b{F}_{a|b}),\end{array}$$

$$(Y_a-Y_b{F}_{a|b}{)}^T(I_n+Y_b{F}_{bb}^{-1}{Y}_b^T{)}^{-1}(Y_a-Y_b{F}_{a|b})\sim W_2({\Sigma }_{aa.b},n).$$

$$r_{t_{ij}}=\frac{(f_{ii}{f}_{jj}{)}^{1/2}{r}_{f_{ij}}+(s_{ii}{s}_{jj}{)}^{1/2}{r}_{s_{ij}}}{(f_{ii}+s_{ii}{)}^{1/2}(f_{jj}+s_{jj}{)}^{1/2}},$$

$$pr\left\{{\text{sBF}}_{ij}^{\text{M}}>b_1\right\}=pr\left\{r_{s_{ij}}^2>c_1\right\},$$

$$r_{q_{ij}}=\frac{(g_{ii}{g}_{ii}{)}^{1/2}{r}_{g_{ij}}+(z_{ii}{z}_{jj}{)}^{1/2}{r}_{z_{ij}}}{(g_{ii}+z_{ii}{)}^{1/2}(g_{jj}+z_{jj}{)}^{1/2}},$$

$$Pr\left\{{\text{sBF}}_{ij}^{\text{C}}>b_2\right\}=Pr\left\{r_{z_{ij}}^2>c_2\right\},$$

$${\text{BFDR}}_t=E({\text{blfdr}}_{ij}\mid {\text{blfdr}}_{ij}<t)=\frac{{\sum }_{i,j}(1-{\pi }_{ij})I_{{\pi }_{ij}>t}}{{\sum }_{i,j}{I}_{{\pi }_{ij}>t}}.$$

$$\begin{array}{r}\iint p_1(Y_a\mid Y_b,B_{a\mid b},{\Sigma }_{aa.b})p_1(B_{a\mid b},{\Sigma }_{aa.b})d{B}_{a\mid b}d{\Sigma }_{aa.b}=\frac{|F_{bb}||F_{aa.b}{|}^{\frac{\delta }{2}}{\Gamma }_2\left(\frac{\delta +n}{2}\right)}{{\pi }^n{\Gamma }_2\left(\frac{\delta }{2}\right)|T_{bb}||T_{aa.b}{|}^{\frac{\delta +n}{2}}}\end{array}.$$

Full text

1 Introduction

Graphical models provide a natural basis for the statistical description and analysis of interplay between variables. In applications, interest often lies in the bidirected and undirected graphs that respectively describe the marginal and conditional dependence structures among variables (Wang,, 2015). When the joint distribution of the variables is assumed to be Gaussian, these are known to be fully coded in the covariance matrix $\Sigma=\{\sigma_{ij}\}$ and its inverse $\Omega=\{\omega_{ij}\}$ (Dempster,, 1972, Cox and Wermuth,, 1993). Precisely, a pair $(i,j)$ of variables with $1\leq i<j\leq p$, will be marginally independent when $\sigma_{ij}=0$ and conditionally independent (given all the remaining variables) when $\omega_{ij}=0$. The present article treats inference of both types of graphs in context of the Gaussian model when the number of variables $p$ is large.

Despite major methodological developments, Bayesian inference for Gaussian graphical models remains challenging. The standard approach casts the problem as a model selection problem, and first requires specification of prior distributions over all possible graphical models and their parameter spaces. Such specification is not straightforward as it is desirable to favour parsimonious models and address the compatibility of priors across models (Carvalho and Scott,, 2009, Consonni and La Rocca,, 2012). Next, the inference procedure is hindered by the search over a very high-dimensional model space where the number of possible graphical models grows super-exponentially with the number of variables. Full exploration of the model space is, therefore, only possible when the number of variables is very small (say $p\leq 10$). In moderate- and high-dimensional settings where $p$ is in the tens, hundreds or thousands, the model space must generally be searched stochastically (Mohammadi and Wit,, 2015, Wang and Li,, 2012, Lenkoski and Dobra,, 2011, Giudici and Green,, 1999). However, due to the tremendous size of the model space in such settings, it may be difficult (nay impossible) to identify with confidence the graphical model that is best supported by the data. Accordingly, it has become common practice to account for model uncertainty by performing Bayesian model averaging and infer the graphical structure by selecting edges with the highest marginal posterior probabilities, for example by exploiting their connection to a Bayesian version of the false discovery rate (Mitra et al.,, 2013, Peterson et al.,, 2015, Baladandayuthapani et al.,, 2014).

To tackle the difficulties associated with the standard approach this article proposes a method to directly select edges by multiple testing of hypotheses about pairwise (marginal or conditional) independence (Drton and Perlman,, 2007) using closed-form Bayes factors. These are obtained using the conditional approach of Dickey, (1971), in which the prior under the null hypothesis is derived from that of the alternative by conditioning on the null hypothesis. This approach was also adopted by Giudici, (1995) to derive a closed-form Bayes factor for conditional independence. However, the latter relies on elements of the inverse of the sample covariance matrix which is singular when the number of variables is large relative to the sample size. We bypass this issue and also introduce a closed-form Bayes factor for marginal independence. Moreover, we show the consistency of the Bayes factors and derive exact tail probabilities from their null distributions to help address the multiplicity problem and control error rates for incorrect edge inclusion. The proposed procedure, available via the R package beam on the CRAN website (R Core Team,, 2016, http://cran.r-project.org), is shown to be computationally very efficient, addressing problems with thousands of nodes in just a few seconds.

The article is structured as follows. The next section introduces basic notations and the Gaussian conjugate model. In section 3 we present closed-form Bayes factors to evaluate the null hypotheses of marginal and conditional independence between any two variables and study their consistency. Section 4 details graph inference and discuss the multiple testing problem and error control. The performance of the proposed approach is compared to Bayesian and non-Bayesian methods on simulated data in section 5. Section 6 illustrates our method on a large gene expression data set from The Cancer Genome Atlas.

2 Background

2.1 Notation

The following notation will be used throughout this paper. We employ the notation $x\mid\mu,\Sigma\sim N_{p}(\mu,\Sigma)$ to say that the random column vector $x\in\mathbb{R}^{p}$, has a multivariate normal distribution with mean $\mu\in\mathbb{R}^{p}$ and positive definite covariance matrix $\Sigma$. We also write $\Omega\mid A,\alpha\sim W_{d}(A,\alpha)$ to indicate that the $d\times d$ random matrix $\Omega$ with density

[TABLE]

has a Wishart distribution with scale matrix $A$ and degrees of freedom $\alpha>d+1$. Here $\mid\Omega\mid$ represents the determinant of $\Omega$, $\text{tr}(A)$ denotes the trace of matrix $A$ and $\Gamma_{d}\left(x\right)=\pi^{\frac{d(d-1)}{4}}\prod_{i=1}^{d}{\Gamma\left(x+(1-i)/2\right)}$ is the multivariate gamma function. The inverse of $\Omega$ is said to have an Inverse-Wishart distribution with scale matrix $A^{-1}$ and $\alpha$ degrees of freedom. We shall use the notation $\Omega^{-1}\mid A^{-1},\alpha\sim IW_{d}(A^{-1},\alpha)$. A random variable $\beta$ following a beta distribution with shape parameters $b_{1}$ and $b_{2}$ will be denoted by $\beta\sim Beta(b_{1},b_{2})$. We use the operator $vec$ to denote the linear transformation that stacks the columns of a matrix into a vector and $\otimes$ to denote the Kronecker product. We refer to Gupta and Nagar, (2000) for more details on these operators. Last, we shall use the subscripts $aa$, $bb$, $ab$ and $ba$ to refer to the submatrices $\Sigma_{aa}$, $\Sigma_{bb}$, $\Sigma_{ab}$ and $\Sigma_{ba}$ of a $p\times p$ symmetric matrix $\Sigma$ whose block-wise decomposition is implied by a partition of its rows and columns into two disjoint subsets indexed by $a\subset\{1,\ldots,p\}$ and $b=\{1,\ldots,p\}\setminus a$.

2.2 The Gaussian conjugate model

Given an $n\times p$ observation matrix $Y=(Y_{1},\ldots,Y_{p})$ the Gaussian conjugate model is defined by

[TABLE]

with $D$ positive definite, $I_{n}$ the $n$-dimensional identity matrix and $\delta>p+1$. Here, the covariance matrix with kronecker product structure makes explicit the assumption of independence for the rows of $Y$ and the dependence of its columns via the covariance $\Sigma$.

Due to conjugacy, model (1) offers closed-form Bayesian estimators of the covariance matrix $\Sigma$ and its inverse $\Omega=\Sigma^{-1}$. The posterior expectation of $\Sigma$ is

[TABLE]

where $S=Y^{T}Y$, and that of its inverse is

[TABLE]

It is important to note that estimator (2) is a linear shrinkage estimator that is a convex linear combination between the maximum likelihood estimator $\widehat{\Sigma}_{\text{mle}}=n^{-1}S$ of $\Sigma$ and $E\left(\Sigma\right)=D$, with weight $\alpha=(\delta-p-1)/(\delta+n-p-1)\in(0,1)$ (Chen,, 1979, Hannart and Naveau,, 2014). Likewise, estimator (3) is recognized as a ridge-type estimator of the precision matrix (Kubokawa and Srivastava,, 2008, Van Wieringen and Peeters,, 2016). The next proposition presents some properties of these two estimators. All proofs are presented in the Appendix.

Proposition 1**.**

Let estimators (2) and (3) depend on $\delta$ with $D$, $n$ and $p$ fixed, and denote them by $\widehat{\Sigma}_{\delta}$ and $\widehat{\Omega}_{\delta}$, respectively. Then the following properties hold:

1. (i)

$\lim_{\delta\rightarrow\infty}\widehat{\Sigma}_{\delta}=D$** 2. (ii)

$\lim_{\delta\rightarrow\infty}\widehat{\Omega}_{\delta}=D^{-1}$** 3. (iii)

$\lim_{\delta\rightarrow p+1}\widehat{\Sigma}_{\delta}=\widehat{\Sigma}_{\text{mle}}$** 4. (iv)

$\lim_{\delta\rightarrow p+1}\widehat{\Omega}_{\delta}=\left\{(n+p+1)/n\right\}\widehat{\Sigma}_{\text{mle}}^{-1}$, if $n>p$, 5. (v)

$\widehat{\Sigma}_{\delta}$* and $\widehat{\Omega}_{\delta}$ are positive definite*

Additionally, the asymptotic properties of estimators (2) and (3) when $n\rightarrow\infty$ are the same as those of the maximum likelihood estimators $\widehat{\Sigma}_{\text{mle}}$ and $\widehat{\Sigma}_{\text{mle}}^{-1}$ of $\Sigma$ and $\Omega$. Proposition 2 summarizes.

Proposition 2**.**

Let estimator (2) and (3) depend on $n$ with $D$, $\delta$ and $p$ be fixed, and denote them by $\widehat{\Sigma}_{n}$ and $\widehat{\Omega}_{n}$, respectively. Then the following properties hold:

1. (i)

$\lim_{n\rightarrow\infty}\widehat{\Sigma}_{n}=\widehat{\Sigma}_{\text{mle}}$** 2. (ii)

$\lim_{n\rightarrow\infty}\widehat{\Omega}_{n}=\widehat{\Sigma}_{\text{mle}}^{-1}$**

2.3 Choice of hyperparameters

In model 1, the prior matrix $D$ represents the prior expectation of $\Sigma$. It may also be interpreted as the shrinkage target towards which the maximum likelihood estimator of the covariance matrix is shrunk, since the posterior expectation of $\Sigma$ is a linear shrinkage estimator. For these reasons, $D$ can be chosen to encourage estimator (2) to have specific structures (e.g. autoregressives or low-ranks). Ideally, in such cases the matrix $D$ should be parameterised by a low-dimensional vector of hyperparameters that are interpretable and for which prior knowledge exists. As often this knowledge is absent, it is common to choose $D=I_{p}$. Throughout this paper we use $D=I_{p}$ and standardize the $n\times p$ observation matrix $Y$ so that for $1\leq j\leq p$, $Y_{j}^{T}1_{n}=0$ and $Y_{j}^{T}Y_{j}/n=1$, where $1_{n}$ is an $n\times 1$ vector whose elements are all equal to 1.

The other hyperparameter $\delta$ clearly acts as a regularization parameter (see equation (2) and (3)) and its value must therefore be chosen carefully. Following Chen, (1979) and Hannart and Naveau, (2014) we use empirical Bayes and estimate $\delta$ by the value maximizing the marginal (or integrated) likelihood of the model, i.e. by

[TABLE]

where

[TABLE]

The above optimization problem is easily solved because the marginal likelihood is concave in $\delta$. Moreover, remark that $\mid\left(\delta-p-1\right)D\mid=\left(\delta-p-1\right)^{p}\prod_{r=1}^{R}{d_{r}}$ and $\mid\left(\delta-p-1\right)D+S\mid=\left(\prod_{r=1}^{R}{d_{r}}\right)\left(\prod_{l=1}^{L}{\left(\delta-p-1+e_{l}\right)}\right)$, where $d_{r}$ and $e_{l}$ are respectively the r${}^{\text{th}}$ and l${}^{\text{th}}$ largest eigenvalues of $D$ and $D^{-1}S$. Hence, evaluating the objective function for different values of $\delta$ is computationally cheap provided the eigenvalues of $D$ and $D^{-1}S$ have been pre-computed. We are referring the reader to Hannart and Naveau, (2014, Section 2.3.) for the proof that the asymptotic properties of estimator (2) and (3) (Proposition 1) hold when $\delta=\widehat{\delta}$.

3 Bayes factors

3.1 Bayes factor for conditional independence

In this section we derive an analytic expression for the Bayes factor evaluating the null hypothesis of conditional independence between two variables in context of model (1). For ease of notation we define $F=(\delta-p-1)D$ and $T=F+S$. We wish to evaluate the null hypothesis of conditional independence, denoted $\text{H}_{0,ij}^{\text{C}}$, between two coordinates $i$ and $j$, $1\leq i<j\leq p$. We test $\text{H}_{0,ij}^{\text{C}}:\omega_{ij}=0$ against the alternative hypothesis $\text{H}_{1,ij}^{\text{C}}:\omega_{ij}\neq 0$, where $\omega_{ij}$ is the $(i,j)^{\text{th}}$ element of $\Omega$. The Bayes factor evaluating evidence in favour of $\text{H}_{1,ij}^{\text{C}}$ is

[TABLE]

where, by definition, $\Sigma_{0}$ is such that $\omega_{ij}=0$.

Giudici, (1995) showed that (4) could be obtained in closed-form by reparameterising the Gaussian conjugate model and defining a compatible prior under the null hypothesis using the approach of Dickey, (1971). However, the proposed Bayes factor does not exist in high dimensional settings because it depends on elements of $S^{-1}$. This problem is here circumvented by factorising the joint likelihood of the observed data as

[TABLE]

the product of a marginal and conditional likelihood. This factorisation arises from the partition of $Y=\left[Y_{a},Y_{b}\right]$ into two disjoint subsets indexed by $a=\{i,j\}$ and $b=V\setminus a$. The quantity $B_{a|b}=\Sigma_{bb}^{-1}\Sigma_{ba}$ represents the matrix of regression coefficients obtained when regressing the variables indexed by $a$ onto the variables indexed by $b$, whereas $\Sigma_{aa.b}=\Sigma_{aa}-\Sigma_{ab}\Sigma_{bb}^{-1}\Sigma_{ba}$ denotes the residual covariance matrix.

The factorisation of the likelihood allows conveniently to simplify (4). Using the change of variable from $(\Sigma_{aa},\Sigma_{ab},\Sigma_{bb})$ to $(\Sigma_{aa.b},B_{a\mid b},\Sigma_{bb})$ together with the fact that $\Sigma_{bb}$ is independent of $(B_{a\mid b},\Sigma_{aa.b})$, most nuisance parameters are integrated out and equation (4) becomes

[TABLE]

Note that by the standard properties of the multivariate normal and Inverse-Wishart distributions Gupta and Nagar, (2000, Theorems 2.3.12. and 3.3.9.) the densities under the alternative model are

[TABLE]

where $F_{a\mid b}=F_{bb}^{-1}F_{ba}$ and $F_{aa.b}=F_{aa}-F_{ab}F_{bb}^{-1}F_{ba}$. Therefore, the simplification of Bayes factor (4) intuitively tells us that evaluating the conditional independence between any two coordinates within the $p$-dimensional Gaussian conjugate model (1) is equivalent to evaluating the diagonality of the residual covariance matrix in a bivariate response regression model.

To obtain (5) in closed-form we, similarly to Giudici, (1995), define a compatible prior for $(B_{a\mid b},\Sigma_{aa.b})$ under the null hypothesis $\text{H}_{0,ij}^{\text{C}}$ using the conditional approach of Dickey, (1971). Precisely, the prior density under $\text{H}_{0,ij}^{\text{C}}$ is derived from that under $\text{H}_{1,ij}^{\text{C}}$ by conditioning on $\text{H}_{0,ij}^{\text{C}}$. The densities under the null model are therefore

[TABLE]

where $\Sigma_{aa.b}^{0}$ is such that $\omega_{ij}=0$.

We now state the main result of this section.

Lemma 1**.**

Assume (5) holds with densities defined by (6) and (7). Then the Bayes factor in favour of $\text{H}_{1,ij}^{\text{C}}$ is

[TABLE]

with F_{aa.b}=\left[\begin{array}[]{cc}g_{ii}&g_{ij}\\ g_{ij}&g_{jj}\end{array}\right], $r_{g_{ij}}=g_{ij}(g_{ii}g_{jj})^{-1/2}$, T_{aa.b}=\left[\begin{array}[]{cc}q_{ii}&q_{ij}\\ q_{ij}&q_{jj}\end{array}\right] and $r_{q_{ij}}=q_{ij}(q_{ii}q_{jj})^{-1/2}$.

Remark 1**.**

In Lemma 1, the quantities $g_{ii}$ and $q_{ii}$ (resp. $g_{jj}$ and $q_{jj}$) can be thought of representing prior and posterior partial variances for coordinate $i$ (resp. $j$), whereas $r_{g_{ij}}$ and $r_{q_{ij}}$ can be thought of representing prior and posterior partial correlations.

Remark 2**.**

The Bayes factor proposed by Giudici, (1995, lemma 3), in contrast to Lemma 1, defines the quantities $g_{ij}$ and $q_{ij}$ such that the matrices $F_{aa.b}=\{g_{ij}\}$ and $F_{aa.b}+S_{aa.b}=\{q_{ij}\}$, with $S_{aa.b}=S_{aa}-S_{ab}S_{bb}^{-1}S_{ba}$. Note that here $S_{aa.b}$ only exists when $S_{bb}$ is invertible (i.e. when $n$ is large relatively to $p$) whereas $T_{aa.b}=T_{aa}-T_{ab}T_{bb}^{-1}T_{ba}$ defined in Lemma 1 exists even when $p>n$ because $T$ is always positive definite (a consequence of Proposition 1).

Remark 3**.**

Standard matrix algebra (Gupta and Nagar,, 2000, Theorem 1.2.3.v) tells us that $F_{aa.b}=\{(F^{-1})_{aa}\}^{-1}$ and $T_{aa.b}=\{(T^{-1})_{aa}\}^{-1}$. This means that the elements of the $2\times 2$ matrices $F_{aa.b}$ and $T_{aa.b}$ can respectively be obtained from the elements of $F^{-1}$ and $T^{-1}$. The computation of the Bayes factor in Lemma 1 for all pairs of variables $(i,j)$, $1\leq i<j\leq p$, hence boils down to computing $F^{-1}$ and $T^{-1}$.

3.2 Bayes factor for marginal independence

We now derive an analytic expression for the Bayes factor evaluating the null hypothesis of marginal independence between any two variables in context of model (1). We test $\text{H}_{0,ij}^{\text{M}}:\sigma_{ij}=0$ against the alternative hypothesis $\text{H}_{1,ij}^{\text{M}}:\sigma_{ij}\neq 0$, where $\sigma_{ij}$ is the $(i,j)^{\text{th}}$ element of $\Sigma$. The Bayes factor evaluating evidence in favour of $\text{H}_{1,ij}^{\text{M}}$ is

[TABLE]

where now $\Sigma^{0}$ is such that $\sigma_{ij}=0$.

We adopt a similar approach as in section 3.1 to obtain (8) in closed-form. We first write the joint likelihood as

[TABLE]

and make a change of variable from $(\Sigma_{aa},\Sigma_{ab},\Sigma_{bb})$ to $(\Sigma_{aa},B_{b\mid a},\Sigma_{bb.a})$, where $B_{b\mid a}=\Sigma_{aa}^{-1}\Sigma_{ab}$ and $\Sigma_{bb.a}=\Sigma_{bb}-\Sigma_{ba}\Sigma_{aa}^{-1}\Sigma_{ab}$. Then, using the fact that $(B_{b\mid a},\Sigma_{bb.a})$ is independent of $\Sigma_{aa}$ it is easily seen that the Bayes factor (8) simplifies to

[TABLE]

Here, the densities under the alternative model, by properties of the multivariate normal and Inverse-Wishart distributions, are

[TABLE]

whereas the densities under the null model are

[TABLE]

We now state the following lemma.

Lemma 2**.**

Assume (9) holds with densities defined by (10) and (11). Then the Bayes factor in favour of $\text{H}_{1,ij}^{\text{M}}$ is

[TABLE]

with F_{aa}=\left[\begin{array}[]{cc}f_{ii}&f_{ij}\\ f_{ij}&f_{jj}\end{array}\right], $r_{f_{ij}}=f_{ij}(f_{ii}f_{jj})^{-1/2}$, T_{aa}=\left[\begin{array}[]{cc}t_{ii}&t_{ij}\\ t_{ij}&t_{jj}\end{array}\right] and $r_{t_{ij}}=t_{ij}(t_{ii}t_{jj})^{-1/2}$.

Remark 4**.**

In Lemma 2, the quantities $f_{ii}$ and $t_{ii}$ (resp. $f_{jj}$ and $t_{jj}$) can be thought of representing prior and posterior marginal variances for coordinate $i$ (resp. $j$), whereas $r_{f_{ij}}$ and $r_{t_{ij}}$ can be thought of representing prior and posterior marginal correlations.

Remark 5**.**

The computation of the Bayes factor in Lemma 2 for all pairs of variables $(i,j)$, $1\leq i<j\leq p$, boils down to computing $T$.

3.3 Consistency

In this section we consider the selection consistency of the Bayes factors defined in Lemma 1 and 2. A Bayes factor is said to be consistent when $\lim_{n\rightarrow\infty}\text{BF}_{ij}=0$ if $\text{H}_{0,ij}$ is true and $\lim_{n\rightarrow\infty}\text{BF}_{ij}=\infty$ if $\text{H}_{1,ij}$ is true (Fernández et al.,, 2001, Casella et al.,, 2009, Wang and Maruyama,, 2016). In other words, the consistency property means that the true hypothesis will be selected when enough data are provided.

To prove the consistency of the Bayes factors, we make the following assumption.

Assumption 1**.**

The sample correlation matrix has a limit as $n\rightarrow\infty$ that is positive definite.

Assumption 1 also appears in Maruyama and George, (2011). We now state the following result.

Lemma 3**.**

Under Assumption 1 the Bayes factors $\text{BF}_{ij}^{\text{C}}$ and $\text{BF}_{ij}^{\text{M}}$ are consistent in selection.

4 Graph structure recovery

4.1 Inference by multiple testing

We propose to infer the marginal and conditional independence graphs by multiple testing of hypotheses using the Bayes factors introduced in the previous section. Precisely, we propose to infer the edge set $E_{U}=\{(i,j)\mid\omega_{ij}\neq 0\}$ of the undirected graph $U=(V,E_{U})$ with vertex set $V$ by evaluating $\text{H}_{0,ij}^{\text{C}}$ versus $\text{H}_{1,ij}^{\text{C}}$ for $1\leq i<j\leq p$. Similarly, we propose to infer the edge set $E_{B}=\{(i,j)\mid\sigma_{ij}\neq 0\}$ of the bidirected graph $B=(V,E_{B})$ by evaluating $\text{H}_{0,ij}^{\text{M}}$ versus $\text{H}_{1,ij}^{\text{M}}$ for $1\leq i<j\leq p$. On the whole, the approach consists in translating the pattern of rejected hypotheses into a graph (Drton and Perlman,, 2007).

4.2 Scaled Bayes factors

To infer either graph structure it is necessary to compare Bayes factors between all $p(p-1)/2$ pairs of variables. However, the Bayes factors defined in Lemma 1 and 2 are not scale-invariant (due to their last terms) and, hence, comparable between different pairs of variables. In light of this, we define scaled versions of the Bayes factors defined in Lemma 1 and 2 that can more appropriately rank edges of graphs $U$ and $B$. Corollary 1 and 2 summarize.

Corollary 1**.**

The scaled Bayes factor in favour of $\text{H}_{1,ij}^{\text{C}}$ is

[TABLE]

with quantities defined as in Lemma 1.

Corollary 2**.**

The scaled Bayes factor in favour of $\text{H}_{1,ij}^{\text{M}}$ is

[TABLE]

with quantities defined as in Lemma 2.

Remark 6**.**

When the prior matrix $D$ is proportional to $I_{p}$, then $r_{f_{ij}}=0$ and $r_{g_{ij}}=0$, and the orderings provided by the scaled Bayes factors in Corollaries 1 and 2 for all pairs $(i,j)$ are identical to the orderings provided by the squares of the posterior marginal and partial correlations $r_{t_{ij}}$ and $r_{q_{ij}}$, respectively.

4.3 Multiplicity adjustment and error control

To address the multiplicity problem we propose to use the tail or error probability associated with the null distribution of each scaled Bayes factor. The tail probability is closely related to the notion of a P-value: the Bayes factor is treated as a random variable and its distribution, which follows that of the random data, is used to make a probability statement about its observed value. Then, to recover the structure of a graph, the tail probabilities obtained from all $p(p-1)/2$ comparisons are adjusted using standard multiplicity correction procedures to control, say, the family-wise error or false discovery rates (Goeman and Solari,, 2014).

In the following, we study the conditional null distribution of the Bayes factors statistics defined in Corollaries 1 and 2. The conditional null distribution here refers to the distribution that would be obtained by shuffling or permuting labels of the observations (Jiang et al.,, 2017). Under this scheme, we shall define $Pr\left(\text{sBF}_{ij}^{\text{M}}>b_{1}\right)$ and $Pr\left(\text{sBF}_{ij}^{\text{C}}>b_{2}\right)$ the probabilities of observing values for the two scaled Bayes factors that are respectively larger than $b_{1}$ and $b_{2}$. Next, we show that these tail probabilities can be obtained analytically without the need of a permutation algorithm.

Before, we state three results which will be used in our argumentation.

Proposition 3**.**

Suppose $\Phi\sim W_{2}(\Sigma,d)$, where

[TABLE]

are parametrised in terms of their correlations $-1\leq\varphi\leq 1$ and $-1\leq\rho\leq 1$. Then,

[TABLE]

Proposition 4**.**

The following equality holds:

[TABLE]

where $\bar{B}_{a\mid b}=(Y_{b}^{T}Y_{b}+F_{bb})^{-1}(Y_{b}^{T}Y_{a}+F_{ba})$.

Proposition 5**.**

Let $\Sigma_{aa.b}$ be fixed. Then, according to model (6) we have

[TABLE]

We observe that the only term of the Bayes factor for marginal independence (defined in Corollary 1) that depends on the data is

[TABLE]

via the elements of $S_{aa}=Y^{T}_{a}Y_{a}=\{s_{ij}\}$. Here $r_{s_{ij}}=s_{ij}(s_{ii}s_{jj})^{-1/2}$. This means that we can write

[TABLE]

where $c_{1}$ is a quantity that depends on $\{\delta,n,f_{ii},f_{jj},r_{f_{ij}},s_{ii},s_{jj}\}$. Now, according to our model in equation (10) it is easily verified that $S_{aa}\sim W_{2}(\Sigma_{aa},n)$ and, using Proposition 3, we can establish that $r_{s_{ij}}^{2}\mid\text{H}_{0,ij}^{\text{M}}\sim Beta(1/2,(n-1)/2)$. The tail probability of the Bayes factor can therefore be computed exactly using $Beta(1/2,(n-1)/2)$. We remark that the definition of the type I error is here conditioning on $\{\delta,n,f_{ii},f_{jj},r_{f_{ij}},s_{ii},s_{jj}\}$.

A similar argument holds for obtaining the tail probability associated with the Bayes factor for conditional independence defined in Corollary 2. The only term of the Bayes factor that depends on the data is $r_{q_{ij}}=q_{ij}(q_{ii}q_{jj})^{-1/2}$, where, we recall, $q_{ij}$ is such that $T_{aa.b}=\{q_{ij}\}$. Proposition 4 suggests that we can write $T_{aa.b}=F_{aa.b}+Z$, with $Z=(Y_{a}-Y_{b}F_{a|b})^{T}(I_{n}+Y_{b}F_{bb}^{-1}Y_{b}^{T})^{-1}(Y_{a}-Y_{b}F_{a|b})$. As a result,

[TABLE]

where $Z=\{z_{ij}\}$ and $r_{z_{ij}}=z_{ij}(z_{ii}z_{jj})^{-1/2}$. This means that we can write

[TABLE]

where $c_{2}$ is a quantity that depends on $\{\delta,n,g_{ii},g_{jj},r_{g_{ij}},z_{ii},z_{jj}\}$. Propositions 3 and 5 tell us that $Z\sim W_{2}(\Sigma_{aa.b},n)$ and that $r_{z_{ij}}^{2}\mid\text{H}_{0,ij}^{\text{C}}\sim Beta(1/2,(n-1)/2)$. Therefore, the tail probability of the Bayes factor can also be computed exactly using $Beta(1/2,(n-1)/2)$. We remark that the definition of the type I error is conditioning on $\{\delta,n,g_{ii},g_{jj},r_{g_{ij}},z_{ii},z_{jj}\}$.

5 Numerical experiments

5.1 Comparison to Bayesian methods

In this section we compare the performance of our approach with other Bayesian methods. For computational reasons, we consider a moderate-dimensional problem. We generate 50 datasets of size $n\in\{25,50,100\}$ from a multivariate Gaussian distribution with mean vector [math] and $50\times 50$ inverse covariance matrix $\Phi$. The matrix $\Phi$ is a sparse matrix which we generate from a G-Wishart distribution with scale matrix equal to the identity and $b=4$ degrees of freedom (using the function bdgraph.sim of R package BDgraph). Four different graph structures are considered which we illustrate in Figure 1.

We compare our method, implemented in the R package beam, to two sampling-based approaches based on the birth-death and reversible jump Markov chain Monte Carlo algorithms, developed by Mohammadi and Wit, (2015, 2017) and implemented in the R package BDgraph, using 100,000 sweeps and a burn-in period of 50,000 updates. We also consider the method of Schwaller et al., (2017), implemented in the R package saturnin, that offers closed-form inference within the class of tree-structured graphical models. For each method we obtain the marginal posterior probabilities of edge inclusion, either via the sampling algorithm or exactly.

To evaluate the performance of the methods in recovering the different graph structures we report the area under the receiver operating characteristic (ROC) curve which depicts the true positive rate, $TPR=TP/(TP+FN)$, as a function of the false positive rate, $FPR=(FP)/(FP+TP)$, overall possible thresholds on the marginal posterior probabilities of edge inclusion (or tail probabilities in case of our method). Here, the quantities $TP$, $FP$, $FN$ denote the number of true positives, false positives and false negatives, respectively. We also report the area under the precision-recall (PR) curve which depict the precision, $PR=TP/(TP+FP)$, as a function of the true positive rate (also referred to as recall).

Table 1 summarizes simulation results. It shows that our method performs well compared to other Bayesian methods in recovering the different graph structures. For instance, our method often achieves the largest areas under the receiver operating characteristic and precision-recall curves for different graph structures and sample sizes. Moreover, a marked improvement is observed in cases where the sample size is small ($n=25$).

The results also show non-negligible differences in performance between the birth-death and reversible jump Markov chain Monte Carlo algorithms. These differences does not seem dependent of the graph structure or sample size. This suggests that the performance of the sampling-based methods can be affected by the choice of sampling algorithm.

Overall, the simulation results demonstrate that our method can recover various graphical structures at least as accurately as other Bayesian approaches at a very low computation cost (see Figure 2). Our method achieves generally a greater area under the precison-recall curve than others. The present results also confirm that obtained by Schwaller et al., (2017), namely the relative good performance of tree-structured graphical models compared to sampling-based approaches despite stronger restrictions on the class of graphical models. However, the performance of the approach can degrade in somes cases (e.g. cluster structures).

In conclusion, we remark that the marginal posterior probabilities of edge inclusion obtained from each method can in principle be linked to a Bayesian version of the false discovery rate to carry out edge selection with error control (Mitra et al.,, 2013, Baladandayuthapani et al.,, 2014, Peterson et al.,, 2015). To see this, let $\pi_{ij}=Pr(\text{H}_{1,ij}^{\text{C}}\mid Y)$ be the marginal posterior probability of inclusion for edge $(i,j)$, $1<i<j<p$, for a given method. Then, its complement $Pr(\text{H}_{0,ij}^{\text{C}}\mid Y)=1-\pi_{ij}=\text{blfdr}_{ij}$ can indeed be viewed as a Bayesian version of the local false discovery rate (Efron et al.,, 2001) where the conditioning is on the data rather than a statistic. This connection serves to define the following Bayesian version of the false discovery rate (Newton et al.,, 2004):

[TABLE]

The control of the false discovery rate therefore relies heavily on the appropriate calibration of the marginal posterior probabilities of edge inclusion. The prior distribution obviously plays an important role in the quality of such calibration, however, the latter may also be affected by the sampling algorithm. Due to inherent differences between the different Bayesian approaches under comparison in this simulation study, it appears difficult to achieve a fair comparison on the control of the false discovery rate. Such comparison is therefore omitted here.

5.2 Comparison to non-Bayesian methods

The performance of the proposed method is compared in higher dimensional settings to non-Bayesian approaches that carry out graphical model selection via multiple testing. We generate 50 datasets of size $n=100$ from a $p$-dimensional Gaussian distribution mean vector [math] and inverse covariance matrix $\Psi$. Throughout the simulation, we fix the sample size $n=100$ and vary of the dimensionality $p\in\{200,500,1000\}$. We consider four different sparse precision matrices corresponding to different graph structures (similar to those illustrated in Figure 1):

band structure: $\Psi_{p}^{\text{band}}$ is a tridiagonal matrix, 2. 2.

cluster structure: $\Psi_{p}^{\text{cluster}}$ is a block diagonal matrix whose diagonal blocks are sparse matrices of size $20$ where the off-diagonal entries of non-zero with probability 0.1. 3. 3.

hub structure: $\Psi_{p}^{\text{hub}}$ is a block diagonal matrix whose diagonal blocks are sparse matrices of size $20$ where only the off-diagonal entries in the first row and column are non-zero, 4. 4.

random structure: $\Psi_{p}^{\text{random}}$ is obtained by randomly permuting the rows and columns of $\Psi_{p}^{\text{band}}$.

For all precision matrices the non-zero entries are generated independently from a uniform distribution on $[-1,1]$ and positive definiteness is ensured by adding a constant to the diagonal so that the minimum eigenvalue is equal to 0.1.

We compare our method to that of Schäfer and Strimmer, (2005), implemented in the R package GeneNet, that is based on a linear shrinkage estimator of the covariance matrix (Ledoit and Wolf,, 2004) and a mixture model for false discovery rate estimation (Strimmer,, 2008). We also consider the asymptotic normal thresholding method of Ren et al., (2015) that is implemented in the R package FastGGM (Wang et al.,, 2016). For both methods we obtain P-values associated with the estimated partial correlations, whereas for our method we use the tail probabilities associated with the Bayes factor defined in Corollary 1 for all pairs of variables.

As in the previous section, we compare the performance of the methods using the areas under the receiver operating characteristic and precision-recall curves.

Table 2 shows that the proposed method performs well in recovering large graphical structures compared to non-Bayesian methods. It achieves comparable areas under the receiver operating characteristic and precision-recall curves as other methods for different problem sizes. However, in the case of hub structures the proposed method performs better.

Besides recovering accurately the different graphical structures, Figure 3 shows that the proposed method is the fastest. When $p=1000$, the average computational time is less than a second whereas contenders are 5 to 20 times slower.

6 Gene network in Glioblastoma multiform

We illustrate our method on a large gene expression data set on glioblastoma multiforme from The Cancer Genome Atlas. Glioblastoma multiform is an aggressive form of brain tumor in adults associated with poor prognosis. Level 3 normalized gene expression data (Agilent 244K platform) from 532 patients were obtained from The Cancer Genome Atlas Data Portal. The data comprise measurements of 17,814 genes, of which 14,827 can uniquely be identified in the PathwayCommons database. Missing expression values were imputed using the Bioconductor R package impute (function impute.knn() with default parameters) and the data standardized as described in section 2.3. A small subset of the data were analyzed in Leday et al., (2017). Instead, we here characterize globally the conditional independence structure between all 14,827 genes.

Figure 4 displays the log-marginal likelihood of model (1) as a function of the prior parameter $\alpha$ when the prior matrix $T$ equals the identity. Using the empirical Bayes estimate of $\alpha$ we computed the Bayes factors and their associated tail probabilities for all pair of variables. These computations took 90 seconds overall on 3.40GHz Intel Core i7-3770 CPU without parallel schemes, which is remarkable for a graph with a total number of 109,912,551 possible edges.

The conditional independence graph identified by controlling the family-wise error rate at 10% using the conservative Bonferroni procedure consists of 46,071 edges (0.042% of the total number of edges). Edge degree varies from [math] to $127$ with 9,675 genes having nonzero degrees. The degree distribution seems to follow an exponential distribution (see Figure 4), thereby indicating that a relative small number of genes have a large number of links.

Because it is difficult to visualize the identified graph in its entirety, we determine clusters of densely connected edges using the algorithm of Blondel et al., (2008) implemented in the R package igraph (Csardi and Nepusz,, 2006). The algorithm identifies a partition of edges that yield an overall modularity score equal to 0.91. The modularity score measures the quality of a division of a graph into sub-graphs. Its maximal value being 1, the identified partition presents a high modularity and suggests the presence of densely interconnected groups of nodes in the conditional independence graph. To illustrate this, we report two sub-graphs in Figure 5 that have been identified by the clustering algorithm and correspond to the HOXA and PCDHB gene families. The HOX gene family is known to be involved in the development of human cancers Bhatlekar et al., (2014), including Glioblastoma. The HOXA13 gene has for instance been advanced as potential diagnostic marker for Glioblastoma (Duan et al.,, 2015) and the role of HOXA9 gene in cell proliferation, apoptosis and drug resistance are under active research (Costa et al.,, 2010, Gonçalves et al.,, 2016, Bhatlekar et al.,, 2018). On the other hand, the protocadherin beta (PCDHB) gene cluster, whose function is still poorly understood, have been reported to be associated with poor survival and tumour aggressiveness in Neuroblastoma (Banelli and Romani,, 2015, Lau et al.,, 2012), another neural cancer. The particular methylation status of genes in the PCDHB family has been identified as a mechanism of transcriptional deregulation and associated with high-risk neuroblastoma biology (Henrich et al.,, 2016).

7 Future work

We foresee several promising extensions of the proposed approach. The Bayes factors proposed in this paper can be used for differential network analysis in which the goal is to identify edges that are in common or specific to predefined groups of samples. Provided that samples between groups are independent, the Bayes factors can simply be multiplied across groups so as to obtain new Bayes factors that provide evidence towards the presence or absence of a common edge. Being symmetric, the Bayes factors can also be inverted before being multiplied so as to evaluate more complex hypotheses, e.g. edge losses or gains in a two-group comparison. Last, it would be interesting to derive the Bayes factor in a regression framework so as to compare them with that of Zhou and Guan, (2017).

Acknowledgements

This research was supported by the Medical Research Council core funding number MRC_MC_UP_0801/1 and grant number MR/M004421. The authors wish to thank Ilaria Speranza for helpful comments on the manuscript and improving largely the software. The first author also wishes to thank Catalina Vallejos and Leonardo Bottolo for helpful discussions.

Appendix A Proofs

This appendix contains the proofs for Lemmas 1, 2 and 3, as well as for Propositions 1, 2, 3, 4 and  5.

Proof of Proposition 1.

Let $\alpha_{\delta}=(\delta-p-1)/(\delta+n-p-1)\in(0,1)$ depends on $\delta$ with $n$ and $p$ fixed.

1. (i)/(iii)

From $\widehat{\Sigma}_{\delta}=\alpha_{\delta}D+(1-\alpha_{\delta})\widehat{\Sigma}_{\text{mle}}$ and the fact that $\lim_{\delta\rightarrow\infty}\alpha_{\delta}=1$ and that $\lim_{\delta\rightarrow p+1}\alpha_{\delta}=0$, it follows immediately that $\lim_{\delta\rightarrow\infty}\widehat{\Sigma}_{\delta}=D$ and $\lim_{\delta\rightarrow p+1}\widehat{\Sigma}_{\delta}=\widehat{\Sigma}_{\text{mle}}$ . 2. (ii)/(iv)

Rewrite $\widehat{\Omega}_{\delta}=\frac{\delta+n}{\delta+n-p-1}\widehat{\Sigma}_{\delta}^{-1},$ then it clear that $\lim_{\delta\rightarrow\infty}\widehat{\Omega}_{\delta}=D^{-1}$ and that $\lim_{\delta\rightarrow p+1}\widehat{\Omega}_{\delta}=\left\{(n+p+1)/n\right\}\widehat{\Sigma}_{\text{mle}}^{-1}$, if $\widehat{\Sigma}_{\text{mle}}$ is positive definite. 3. (v)

Since by definition $x^{T}Dx>0$ and $x^{T}\widehat{\Sigma}_{\text{mle}}x\geq 0$, $\forall x\in\mathbb{R}^{p}$, it follows that $x^{T}\widehat{\Sigma}_{\delta}x=\alpha_{\delta}x^{T}Dx+(1-\alpha_{\delta})x^{T}\widehat{\Sigma}_{\text{mle}}x>0$. Therefore, $\widehat{\Sigma}_{\delta}$ is positive definite for $0<\alpha_{\delta}\leq 1$ or equivalently $\delta>0$. As a direct consequence, $\widehat{\Omega}_{\delta}=\frac{\delta+n}{\delta+n-p-1}\widehat{\Sigma}_{\delta}^{-1}$ is also found positive definite.

∎

Proof of Proposition 2.

Let $\alpha_{n}=(\delta-p-1)/(\delta+n-p-1)\in(0,1)$ depends on $n$ with $\delta$ and $p$ fixed.

1. (i)

From $\widehat{\Sigma}_{\delta}=\alpha_{n}D+(1-\alpha_{n})\widehat{\Sigma}_{\text{mle}}$ and the fact that $\lim_{n\rightarrow\infty}\alpha_{n}=0$, it follows immediately that $\lim_{n\rightarrow\infty}\widehat{\Sigma}_{n}=\widehat{\Sigma}_{\text{mle}}$. 2. (ii)

Rewriting $\widehat{\Omega}_{n}=\frac{\delta+n}{\delta+n-p-1}\widehat{\Sigma}_{\delta}^{-1}$ it is clear that $\lim_{n\rightarrow\infty}\widehat{\Omega}_{n}=\widehat{\Sigma}_{\text{mle}}^{-1}$.

∎

Proof of Lemma 1.

The numerator of the Bayes factor is

[TABLE]

Under $\text{H}_{0,ij}^{\text{C}}$, the model likelihood is

[TABLE]

and the probability density of $(\Sigma_{aa.b}^{0},B_{a|b})$ is

[TABLE]

where $B_{a|b}^{(l)}$ and $Y_{a}^{(l)}$ represent the $l^{\text{th}}$ column of $B_{a|b}$ and $Y_{a}$. Therefore,

[TABLE]

Combining (12) and (13) we obtain the Bayes factor in Lemma 1. ∎

Proof of Lemma 2.

The numerator of the Bayes factor is

[TABLE]

Under $\text{H}_{0,ij}^{\text{M}}$, the model likelihood is

[TABLE]

and the probability density of $\Sigma_{aa}^{0}$ is

[TABLE]

As result,

[TABLE]

Combining (14) and (15) we obtain the Bayes factor in Lemma 2. ∎

Proof of Lemma 3.

Using Stirling’s formula, the gamma function can asymptotically be approximated by

[TABLE]

for large values of $x$ (Wang and Maruyama,, 2016). This means that for large values of $n$,

[TABLE]

and

[TABLE]

The Bayes factor in Lemma 2 is therefore asymptotically equivalent to

[TABLE]

Now when $\text{H}^{\text{M}}_{0,ij}$ is true, $\lim_{n\rightarrow\infty}r^{2}_{s_{ij}}=0$ because the sample correlation is asymptotically unbiased. Hence, $\lim_{n\rightarrow\infty}(1-r_{s_{ij}}^{2})^{-\frac{\delta+n-p+2}{2}}=1$. Since $\lim_{n\rightarrow\infty}n^{3/2}=\infty$, we can conclude that $\lim_{n\rightarrow\infty}\text{BF}_{ij}^{\text{M}}=0$, which proves the consistency under $\text{H}^{\text{M}}_{0,ij}$. On the other hand, when $\text{H}^{\text{M}}_{1,ij}$ is true, $\lim_{n\rightarrow\infty}(1-r_{s_{ij}}^{2})^{-1}=c_{1}$, $c_{1}>0$. The application of L’Hôpital’s rule twice, by deriving the numerator and denominator twice with respect to $n$, allows us to conclude that $\lim_{n\rightarrow\infty}\text{BF}_{ij}^{\text{M}}=\infty$. This completes our proof for the consistency of $\text{BF}_{ij}^{\text{M}}$.

Using similar arguments we prove the consistency of the Bayes factor in Lemma 1. The latter is asymptotically equivalent to

[TABLE]

where $r_{p_{ij}}$ denotes the sample partial correlation between variables $i$ and $j$. On one hand, when $\text{H}^{\text{C}}_{0,ij}$ is true $\lim_{n\rightarrow\infty}r^{2}_{p_{ij}}=0$ because the sample partial correlation is asymptotically unbiased. This implies that $\lim_{n\rightarrow\infty}(1-r_{p_{ij}}^{2})^{-\frac{\delta+n-p+2}{2}}=1$. And because $\lim_{n\rightarrow\infty}n^{3/2}=\infty$ we can conclude that $\lim_{n\rightarrow\infty}\text{BF}_{ij}^{\text{C}}=0$. On the other hand, when $\text{H}^{\text{C}}_{1,ij}$ is true, $\lim_{n\rightarrow\infty}(1-r_{p_{ij}}^{2})^{-1}=c_{2}$, $c_{2}>0$. Therefore, by applying L’Hôpital’s rule twice, as above, it is found that $\lim_{n\rightarrow\infty}\text{BF}_{ij}^{\text{C}}=\infty$, which completes our proof for the consistency of $\text{BF}_{ij}^{\text{C}}$. ∎

Proof of Proposition 3.

It follows directly from the probability density of $\Phi$ that $p(\varphi\mid\rho=0)\propto(1-\varphi^{2})^{\frac{d-3}{2}}$, which implies that $p(\varphi^{2}\mid\rho=0)\propto\varphi^{-1}(1-\varphi^{2})^{\frac{d-3}{2}}$. ∎

Proof of Proposition 4.

After some algebra we have

[TABLE]

where, we recall, $\bar{B}_{a\mid b}=(Y_{b}^{T}Y_{b}+F_{bb})^{-1}(Y_{b}^{T}Y_{a}+F_{ba})$. Now,

[TABLE]

and, using twice the Sherman-Morrison-Woodbury matrix identity (Gupta and Nagar,, 2000, Theorem 1.2.3.iv), it is found that

[TABLE]

Thus, by plugin in (17) and (18) in (16), we obtain

[TABLE]

which is further reduced to

[TABLE]

using, again, the Sherman-Morrison-Woodbury matrix identity. ∎

Proof of Proposition 5.

Consider the case where $\Sigma_{aa.b}$ is fixed, then the joint density of model (6) is

[TABLE]

Here the last equality is obtained using Proposition 4. As a result,

[TABLE]

and $\text{vec}(Y_{a})\mid Y_{b},\Sigma_{aa.b}\sim N_{n\times 2}\left(\text{vec}(Y_{b}F_{a\mid b}),\Sigma_{aa.b}\otimes\left(I_{n}+Y_{b}F_{bb}^{-1}Y_{b}^{T}\right)\right)$.

Now, by the scaling property of the multivariate Gaussian distribution it is found that

[TABLE]

and it follows that (Gupta and Nagar,, 2000, Theorem 3.2.2.)

[TABLE]

∎