Bayesian Analysis of Nonparanormal Graphical Models Using Rank-Likelihood

TL;DR

This paper introduces a Bayesian method for nonparanormal graphical models that uses a rank-likelihood to avoid modeling transformation functions, employing a horseshoe prior and Gibbs sampling for efficient inference.

Contribution

It proposes a novel Bayesian approach utilizing rank-likelihood for nonparanormal models, with theoretical consistency and practical efficiency demonstrated.

Findings

The method achieves accurate graph structure recovery in simulations.

Posterior consistency is established for the precision matrix.

The approach performs well on real data applications.

Abstract

Gaussian graphical models, where it is assumed that the variables of interest jointly follow a multivariate normal distribution with a sparse precision matrix, have been used to study intrinsic dependence among variables, but the normality assumption may be restrictive in many settings. A nonparanormal graphical model is a semiparametric generalization of a Gaussian graphical model for continuous variables where it is assumed that the variables follow a Gaussian graphical model only after some unknown smooth monotone transformation. We consider a Bayesian approach for the nonparanormal graphical model using a rank-likelihood which remains invariant under monotone transformations, thereby avoiding the need to put a prior on the transformation functions. On the underlying precision matrix of the transformed variables, we consider a horseshoe prior on its Cholesky decomposition and use an efficient posterior Gibbs sampling scheme. We present a posterior consistency result for the precision matrix based on the rank-based likelihood. We study the numerical performance of the proposed method through a simulation study and apply it on a real dataset.

Equations

$$Y_{i,1},\dots ,Y_{i,p}=H_1(X_{i,1}),\dots ,H_p(X_{i,p})\sim {\mathrm{N}}_p(0,\text{C})$$

$$D=\{\mathbf{Y}\in \text{mathds}{R}^{n\times p}:y_{i_{r-1},j}<y_{i_r,j}<y_{i_{r+1},j}\}.$$

$$\text{Pr}(\mathbf{Y}\in D|\mathbf{C},H_1,\dots ,H_p)={\int }_{D}p(\mathbf{Y}|\mathbf{C})d\mathbf{Y}=\text{Pr}(\mathbf{Y}\in D|\mathbf{C}).$$

$$Y_d=\sum _{k>d}{\beta }_{kd}{Y}_k+{ϵ}_d,{ϵ}_d\sim \mathrm{N}(0,{\sigma }_d^2),d=1,\dots ,p,$$

$${\mathbf{Y}}_d|({\mathbf{Y}}_{k>d},{\mathbf{\beta }}_{k>d},{\sigma }_d^2)\sim \mathrm{N}({\mathbf{Y}}_{k>d}{\mathbf{\beta }}_{k>d},{\sigma }_d^2\mathbf{I}),$$

$${\beta }_{kd}|({\lambda }_d^2,b_{kd},{\sigma }_d^2)\stackrel{\mathrm{ind}}{\sim }\mathrm{N}(0,\frac{{\sigma }_d^2{b}_{kd}{c}^2{\lambda }_d^2}{p^{2}k}),$$

$${\lambda }_d^2|a_d\sim \mathrm{IG}(\frac{1}{2},\frac{1}{a_d}),$$

$$a_d\sim \mathrm{IG}(\frac{1}{2},1),$$

$$b_{kd}|h_{kd}\stackrel{\mathrm{ind}}{\sim }\mathrm{IG}(\frac{1}{2},\frac{1}{h_{kd}}),$$

$$h_{kd}\sim \mathrm{IG}(\frac{1}{2},1),$$

$${\sigma }_d^2\sim \mathrm{IG}(0.01,0.01).$$

$${\rho }_{kd,m}=-\frac{{\psi }_{kd,m}}{\sqrt{{\psi }_{kd,m}{\psi }_{dd,m}}},$$

$${\varphi }_{kd,m}=-\frac{{\lambda }_{kd,m}}{\sqrt{{\lambda }_{kd,m}{\lambda }_{dd,m}}},$$

$$\frac{{\rho }_{kd,m}}{{\varphi }_{kd,m}}>0.5$$

$$\underset{\mathbf{\Psi }}{\mathrm{minimize}}-n\log \det \mathbf{\Psi }+\text{tr}(\mathbf{\Psi }\text{S}),\text{subject to }\mathcal{C}(\mathbf{\Psi }),$$

$$\underset{n\to \infty }{\lim }\pi (\mathbf{\Psi }\in {\mathcal{U}}_0|\mathbf{R})=1\text{mbox}a.e.[P_{\mathbf{\Psi },H}^{\infty }],$$

$$\mathrm{E}(U_{ik}-\frac{R_{n,ik}}{n+1}{)}^2=\frac{1}{n}\sum _{j=1}^n\mathrm{E}[(U_{i,k}-\frac{j}{n+1}{)}^2|R_{n,ik}=j]=\frac{1}{n}\sum _{j=1}^n\frac{j(n-j+1)}{(n+1{)}^2(n+2)}<\frac{1}{n}.$$

$$U_{ik}=\underset{n^{\prime }\to \infty }{\lim }\frac{R_{n^{\prime },ik}}{n^{\prime }+1}\text{mbox}fori\ge 1,\text{mbox}withprobability1forsomesubsequence\{n^{\prime }\}$$

$$\text{Specificity}=\frac{\text{TN}}{\text{TN}+\text{FP}},\text{Sensitivity}=\frac{\text{TP}}{\text{TP}+\text{FN}},$$

$$\text{MCC}=\frac{\text{TP}\times \text{TN}-\text{FP}\times \text{FN}}{\sqrt{(\text{TP}+\text{FP})(\text{TP}+\text{FN})(\text{TN}+\text{FP})(\text{TN}+\text{FN})}},$$

$$\text{Scaled L1-loss}=\frac{1}{p^2}\sum _k\sum _d\Vert {\widehat{\mathbf{{\rm Y}}}}_{kd}-{\mathbf{{\rm Y}}}_{{\text{true}}_{kd}}\Vert .$$

Full text

††journal: Journal of LaTeX Templates

1 Introduction

Graphical models are useful mathematical tools that model complex dependence relationships between variables. Under the Gaussianity assumption, the graph of relations is completely determined by the zeros in the inverse covariance matrix, also known as the precision matrix. If the $(i,j)$th entry of the precision matrix is zero, then the $i$th and the $j$th variables are conditionally independent given all other variables; see Lauritzen [1996] and Edwards [2000] for more information on the properties of Gaussian graphical models (GGMs).

If the variables are not normally distributed, an adjustment needs to be made if one wants to use the statistical properties of a Gaussian graphical model. The standard procedure is to transform the data, typically by the logarithmic or the square root transformation. With such transformations, one has to check each transformation to see if the data appear close to normal, checking which can be tedious and the method is somewhat ad-hoc. From a modeling perspective, an easier solution is to leave the transformation functions unspecified and estimate them instead. The nonparanormal graphical model [Liu et al., 2009] consists of estimating the transformation functions using a truncated marginal empirical distribution function, and then estimating the precision matrix of the transformed variables using the graphical lasso assuming sparsity. A Bayesian approach for nonparanormal graphical models, developed by Mulgrave and Ghosal [2018a, b], uses a random series B-splines prior to estimate the transformation functions and induces sparsity on the off-diagonal entries of the precision matrix by spike-and-slab or continuous shrinkage priors, either directly or through the Cholesky decomposition.

However, one important requirement of the nonparanormal model is that the transformation functions need to be estimated. In order to avoid estimating the transformation functions, alternative rank-based procedures can be employed to transform the data to normally distributed random variables. Dabrowska and Doksum [1988] study the properties of the likelihood function based on transformed observations in a general class called transformation models. An example is given by the rank-likelihood which is the joint distribution of the ranks of the observations. The rank-likelihood is invariant under monotone transformations of the observations. Thus one can ignore the transformations and focus on the main parameter of interest, the precision matrix, or equivalently, the inverse correlation matrix. The rank-likelihood has been used for semiparametric copula estimation [Hoff, 2007], for Receiver Operating Characteristic (ROC) curve estimation [Gu and Ghosal, 2009, Gu et al., 2014], and for ROC surface estimation [Zhu and Ghosal, 2019]. The use of the rank-likelihood results in a gain in robustness and simplification of estimation. Models which can be represented by transformations on Gaussian graphical models are generally known as Gaussian copula graphical models (GCGMs) and have been explored in the Bayesian literature [Pitt et al., 2006, Dobra and Lenkoski, 2011, Liu et al., 2012, Mohammadi and Wit, 2017]. In the frequentist literature, nonparametric rank-based correlation coefficient estimators have been used to construct GCGMs [Liu et al., 2012, Xue and Zou, 2012]. These models can also address binary or ordinal data, but we do not pursue this direction here.

Once the transformed variables have been obtained, one needs to estimate a sparse precision matrix in order to learn the structure of the graphical model. In the frequentist literature, the popular algorithm to use is the graphical lasso [Friedman et al., 2008]. Numerous algorithms have been proposed to solve this problem [Meinshausen and Buhlmann, 2006, Yuan and Lin, 2007, Friedman et al., 2008, Rothman et al., 2008, Banerjee et al., 2008, d’Aspremont et al., 2008, Lu, 2009, Scheinberg et al., 2010, Witten et al., 2011, Mazumder and Hastie, 2012b]. Alternative penalties used to estimate the sparse precision matrix include adaptive LASSO and SCAD penalties [Fan et al., 2009], and LASSO and grouped LASSO penalties [Friedman et al., 2010]. In the Bayesian literature, graphical models can be learned with the use of continuous shrinkage priors. Priors such as the double exponential [Wang, 2012, Peterson et al., 2013], uniform shrinkage priors [Wang and Pillai, 2013], and spike-and-slab priors [Wang, 2015, Peterson et al., 2016, Li and McCormick, 2017, Li et al., 2017, Mulgrave and Ghosal, 2018a] can characterize zero and non-zero elements of the precision matrix with continuous distributions that have a mass at zero and heavy tails. More recently, horseshoe priors have been employed for sparse precision estimation [Mulgrave and Ghosal, 2018b, Williams et al., 2018]. We estimate a sparse precision matrix using a Cholesky decomposition to naturally incorporate the positive definite matrix, a horseshoe prior to regularize the matrix, and a loss procedure to threshold the matrix. Our methods differ from other Bayesian GCGMs which use the rank-likelihood to transform the random variables. Dobra and Lenkoski [2011] use a G-Wishart prior [Roverato, 2002] and Markov Chain Monte Carlo (MCMC) to estimate the sparse precision matrix to construct a GCGM. Mohammadi et al. [2017] estimate a sparse precision matrix using a G-Wishart prior and birth-and-death MCMC [Mohammadi and Wit, 2015]. Li and McCormick [2017] put a normal spike-and-slab prior on the precision matrix and use an expectation-conditional maximization algorithm for estimation.

The paper is organized as follows. In the next section, we state model and prior. In Section 3, we review the posterior computation. In Section 4, we review the thresholding procedure and in Section 5, we discuss our tuning procedure. We derive a posterior consistency result in Section 6 and in Section 7, we review the results of the simulation study. Lastly, in Section 8, we describe an application with gene expression data.

2 Model and Prior

2.1 Estimation of Transformed Variables

Given a set of observed continuous variables $X_{i,1},\ldots,X_{i,p}$ distributed from unknown marginal distributions, there exist monotone increasing transformation functions $H_{1},\ldots,H_{p}$, such that the distribution of the continuous transformed variables is

[TABLE]

where $\bm{C}$ is the correlation matrix. We make this model identifiable by centering the transformed variables and setting the covariance matrix equal to the correlation matrix. We wish to make inference on $\bm{\Psi}=\textbf{C}^{-1}$ and not on the transformation functions $H_{1},\ldots,H_{p}$. Let $\bm{R}$ represent the array of ranks of $\bm{X}$. Since the transformation functions are increasing, $X_{i_{r-1},j}<X_{i_{r},j}$ implies that $Y_{i_{r-1},j}<Y_{i_{r},j}$, where $i_{r}$ is the index for the position of the $r$th smallest observation of the $j$th component of $X$, such that $i=1,\ldots,n$ and $j=1,\ldots,p$. Then for a given data set $\bm{X}=(\bm{X}_{1},\ldots,\bm{X}_{p})^{\prime}$, the transformed variables $\bm{Y}=(\bm{Y}_{1},\ldots,\bm{Y}_{p})^{\prime}$ must lie in the set

[TABLE]

Then following Dabrowska and Doksum [1988] and the notation of Hoff [2007], we calculate the rank-likelihood as

[TABLE]

Thus, this likelihood depends on the parameter $\bm{C}$ and not on the nuisance functions $H_{1},\ldots,H_{p}$.

Since we wish to make inference on $\bm{\Psi}=\textbf{C}^{-1}$, we sample $\bm{\Psi}$. We reparameterize the model in terms of the non-identifiable inverse covariance matrix $\bm{\Omega}$, but focus our posterior inference on the identifiable inverse correlation matrix. Thus $\bm{\Psi}=\bm{A}\bm{\Omega}\bm{A}$, where $\bm{A}=\textup{diag}(\sigma_{dd})$ and $\sigma_{dd}$ are the square roots of the diagonal elements of $\bm{\Sigma}=\bm{\Omega}^{-1}$. The rank-likelihood is scale invariant so the non-identifiable and identifiable models lead to the same posterior distribution, $\textup{Pr}(\bm{\Psi}|\bm{Y}\in D)\propto p(\bm{\Psi})\times p(\bm{Y}\in D|\bm{\Psi})$.

2.2 Estimation of Inverse Correlation Matrix

We put a horseshoe prior on $\bm{\Omega}$ and sample $\bm{\Omega}$ using a regression-based Cholesky decomposition method discussed in Mulgrave and Ghosal [2018b]. Denote the Cholesky decomposition of $\bm{\Omega}$ as $\bm{\Omega}=\bm{L}\bm{L}^{T}$, where $\bm{L}=(\!(l_{kd})\!)$ is a lower triangular matrix. Define the coefficients $\beta_{kd}=-l_{kd}/l_{dd}$ and the precision as $\phi_{d}=1/\sigma_{d}^{2}=l_{dd}^{2}$. Then the multivariate Gaussian model $\bm{Y}\sim\mathrm{N}_{p}(\mathbf{0},\bm{\Sigma}$), where $\bm{\Sigma}=\bm{\Omega}^{-1},$ leads to the set of independent regression problems,

[TABLE]

where $\beta_{kd}$ are the regression coefficients for $k=d+1,\ldots,p$ and $d=1,\ldots,p$. Denoting $\bm{Y}_{d}$ to be the $d$th column of $\bm{Y}$, $\bm{Y}_{k>d}$ the matrix formed by columns of $\bm{Y}$ greater than $d$, and $\bm{\beta}_{k>d}=(\beta_{d+1},\ldots,\beta_{p})$, we may write the regression relation in the vector form as

[TABLE]

which gives rise to the likelihood.

We use a standard conjugate noninformative prior on the variances with improper density proportional to $\sigma_{d}^{-2}$. We enforce a sparsity constraint along the rows of the lower triangular matrix in order to ensure that the probability that an entry is nonzero (i.e. sparsity) remains roughly the same over different rows. We choose $\rho_{k}=$Prob(nonzero in $k$th row)$={c}/({p\sqrt{k}})$, and tune the value of $c\in\{0.1,1,10,100\}$ to cover a range of four orders of magnitude, i.e. $10^{-1},\;10^{0},\;10^{1},\;10^{2}$; see Mulgrave and Ghosal [2018b] for the more information on the sparsity constraint.

We use a horseshoe prior on $\beta_{kd}$ described in Neville et al. [2014]

[TABLE]

for $k,d=1,\ldots,p$; here IG stands for the inverse gamma distribution.

According to van der Pas et al. [2014], the global scale parameter $\lambda$ is roughly equivalent with the probability of a nonzero element. We enforce the sparsity constraint by replacing $\lambda$ with $(\lambda_{d}c)/(p\sqrt{k})$. Thus, since we are working with the squared parameter, the factor in the variance term for $\beta_{kd}$ is $(\lambda^{2}c^{2})/(p^{2}k)$. The prior on $\bm{\Omega}$ leads to an induced prior on $\bm{\Psi}$.

3 Posterior Computation

We obtain samples of $\textup{Pr}(\bm{\Psi}|\bm{Y}\in D)$ by employing the following Gibbs sampler:

Sample $\bm{Y}\sim\mathrm{N}(\mathbf{0},\bm{\Psi}^{-1}):$

For $d=1,\ldots,p$,

1. (a)

Compute the ranks of $\bm{X}_{d}$, where $r=1,\ldots,n$ represents the ranks.

For $r=1,\ldots,n,$

1. (i)

Compute $y_{\textup{lower}}=y_{i_{r-1},d}$ and $y_{\textup{upper}}=y_{i_{r+1},d}$, where $i$ is the index for the position of the $r$-th rank of $\bm{X}_{d}$. Let $y_{i_{0},d}=-\infty$ and $y_{i_{n+1},d}=\infty$. 2. (ii)

Compute $\mu_{i_{r},d}=-\bm{\psi}_{d,d}^{-1}\bm{\Psi}_{d,\setminus d}\bm{Y}_{i_{r},\setminus d}$ 3. (iii)

$\sigma_{d}^{2}=\bm{\psi}_{d,d}^{-1}$ 4. (iv)

Sample $Y_{i_{r},d}\sim\mathrm{TN}(\mu_{i_{r},d},\sigma_{d}^{2};\;y_{\textup{lower}}<Y_{i_{r},d}<y_{\textup{upper}}),$ where

$\mathrm{TN}(a,b;\;c<x<d)$ denotes a univariate truncated normal distribution with mean $a$, variance $b$, and truncation limits $c$ and $d$. Sampling from the truncated normal distribution is implemented with the fast sampling function trandn in MATLAB that uses minimax tilting [Botev, 2017]. 2. 2.

Sample $\bm{\Omega}$:

For $d=1,\ldots,p-1,$

1. (a)

Sample the variables $\bm{\beta}_{k>d}|(\sigma_{d},\bm{b}_{k>d},\lambda_{d}^{2})\sim\mathrm{N}(\bm{A}^{-1}\bm{Y}_{k>d}^{T}\bm{Y}_{d},\sigma_{d}^{2}\bm{A}^{-1})$, where

$\bm{A}=(\bm{Y}_{k>d}^{\prime}\bm{Y}_{k>d}+\textup{diag}({p^{2}k}/({\lambda_{d}^{2}\bm{b}_{k>d}c^{2}}))).$

Since sampling from this normal distribution can be expensive with large $p$, we used an exact sampling algorithm for Gaussian priors that uses data augmentation [Bhattacharya et al., 2016]:

1. (i)

Sample $t\sim\mathrm{N}(\mathbf{0},\bm{D})$ and $\delta\sim\mathrm{N}(0,I_{n})$, where $\bm{D}=\displaystyle\sigma_{d}^{2}\textup{diag}(\frac{\lambda_{d}^{2}\bm{b}_{k>d}c^{2}}{p^{2}k});$ 2. (ii)

set $v=\bm{\Phi}t+\delta$, where $\bm{\Phi}=\bm{Y}_{k>d}/\sigma_{d};$ 3. (iii)

solve $(\bm{\Phi}\bm{D}\bm{\Phi}^{\prime}+\bm{I}_{n})w=(\alpha-v)$, where $\alpha=\textbf{Y}_{d}/\sigma_{d};$ 4. (iv)

set $\beta=t+\bm{D}\bm{\Phi}^{\prime}w.$ 2. (b)

Sample $\lambda_{d}^{2}\sim\mathrm{IG}(\displaystyle\frac{\#(k>d)}{2}+\frac{1}{2},\frac{1}{2}\bm{\beta}_{k>d}^{\prime}\textup{diag}(\frac{p^{2}k}{\sigma_{d}^{2}\bm{b}_{k>d}c^{2}})\bm{\beta}_{k>d}+\frac{1}{a_{d}}).$ 3. (c)

Sample $a_{d}\sim\mathrm{IG}(1,\lambda_{d}^{-2}+1).$ 4. (d)

Sample $b_{kd}\sim\mathrm{IG}(1,\frac{p^{2}k\beta_{kd}^{2}}{2\sigma_{d}^{2}\lambda_{d}^{2}c^{2}}+\frac{1}{h_{kd}}).$ 5. (e)

Sample $h_{kd}\sim\mathrm{IG}(1,\frac{1}{b_{kd}}+1).$ 6. (f)

Sample $\sigma_{d}^{2}\sim\mathrm{IG}(\frac{n+\#(k>d)}{2}+0.01,\frac{1}{2}\left\|\bm{Y}_{d}-\bm{Y}_{k>d}\bm{\beta}_{k>d}\right\|^{2}$

$+\frac{1}{2}\bm{\beta}_{k>d}^{\prime}\textup{diag}(\frac{p^{2}k}{\lambda_{d}^{2}\bm{b}_{k>d}c^{2}})\bm{\beta}_{k>d}+0.01).$ Sample $\sigma_{p}^{2}\sim\mathrm{IG}({n}/{2}+0.01,1/2\left\|\bm{Y}_{p}\right\|^{2}+0.01).$ 7. (g)

Compute $l_{kd}=-\beta_{kd}/\sigma_{d}\textup{ and }l_{dd}=1/\sigma_{d}.$ 8. (h)

Compute $\bm{\Omega}=\textbf{L}\textbf{L}^{\prime}.$ 3. 3.

Set $\bm{\Psi}=\bm{A}\bm{\Omega}\bm{A}$, where $\bm{A}=\textup{diag}(\sigma_{dd})$ and $\sigma_{dd}$ are the square roots of the diagonal elements of $\bm{\Sigma}=\bm{\Omega}^{-1}$.

These steps are repeated until convergence.

4 Thresholding

4.1 0-1 Loss Procedure

We find the posterior partial correlation using the inverse correlation matrices from the Gibbs sampler of the horseshoe prior (2) and the posterior partial correlation using the standard conjugate Wishart prior. The posterior partial correlation using the matrices from the Gibbs sampler is defined as

[TABLE]

where $\psi_{kd,m}$ is the $m$th sample of $M$ Markov chain Monte Carlo (MCMC) draws after burn-in from the posterior distribution of $\psi_{kd}$, $k,d=1,\ldots,p$, $m=1,\ldots,M$. The posterior partial correlation using the standard conjugate Wishart prior is found by starting with the latent observation, $\bm{Y}_{m},$ which is obtained from the MCMC output. We put a standard Wishart prior on the inverse correlation matrix, $\bm{\Psi}\sim\mathrm{W}_{p}(3,\bm{I})$, where $\bm{I}$ is the identity matrix. By conjugacy, the posterior is $\mathrm{W}_{p}(n+3,(\bm{I}+\bm{S})^{-1})$, where $\bm{S}=\bm{Y}^{\prime}\bm{Y}$. We compute the mean of the posterior distribution given $\bm{Y}$, $\bm{\Lambda}=\mathds{E}(\bm{\Psi}|\bm{Y})=(n+3)(\bm{I}_{p}+\bm{S})^{-1}$. Finally, we find the posterior partial correlation coefficients

[TABLE]

where $\lambda_{kd,m}$ is the $(k,d)$th entry of $\bm{\Lambda}$ at the $m$th MCMC iteration.

We link these two posterior partial correlations for the 0-1 loss method. Our convention is that the event $\{\psi_{kd,m}\neq 0\}$ holds if and only if

[TABLE]

for $k,d=1,\ldots,p$ and $m=1,\ldots,M$. The rationale for this thresholding procedure is that we are comparing the regularized precision matrix to the non-regularized precision matrix from the Wishart prior. If the absolute value of the partial correlation coefficient from the regularized precision matrix is similar in size or larger than the absolute value of the partial correlation coefficient from the Wishart precision matrix, then the edge matrix should have an edge. If the absolute value of the partial correlation coefficient from the regularized precision matrix is much smaller than the absolute value of the coefficient from the Wishart matrix, then the edge matrix should not have an edge. The precision matrix from the Wishart prior serves as a means of comparison to determine whether the element of the regularized precision matrix is truly large or small.

5 Choice of Prior Parameters

For the sparsity constraint on the inverse correlation matrix $\bm{\Psi}$, we need to select the value of the parameter $c$. We solve a convex constrained optimization problem in order to use the Bayesian Information Criterion (BIC), originally described in Dahl et al. [2005, 2008] and developed further in Mulgrave and Ghosal [2018a]. First, we find the Bayes estimate of the inverse correlation matrix, $\hat{\bm{\Psi}}=\mathrm{E}(\bm{\Psi}|\bm{Y})$. We also find the average of the transformed variables, $\bar{\bm{Y}}=M^{-1}\sum_{m=1}^{M}\bm{Y}_{m}$, where $\bm{Y}_{m}$, $m=1,\ldots,M$, are obtained from the MCMC output. Then, using the sum of squares matrix, $\bm{S}=\bar{\bm{Y}}^{\prime}\bar{\bm{Y}}$, we solve the following to obtain the maximum likelihood estimate of the inverse correlation matrix, $\hat{\bm{\Psi}}_{\mathrm{MLE}}$:

[TABLE]

where $\mathcal{C}$ stands for the constraint that an element of $\bm{\Psi}$ is zero if and only if the estimated edge matrix from the MCMC sampler has that element as zero. The estimated edge matrix from the MCMC sampler will be described in more detail in Subsection 7.1. For computational simplicity, in the code, we represent this problem as an unconstrained optimization problem as described in Dahl et al. [2005, 2008].

Lastly, we calculate $\textup{BIC}=-2\ell(\hat{\bm{\Psi}}_{\mathrm{MLE}})+k\log n$, where $k=\#\mathcal{C}(\hat{\bm{\Psi}}_{\mathrm{MLE}})$, the sum of the number of diagonal elements and the number of edges in the estimated edge matrix, and $-\ell(\hat{\bm{\Psi}}_{\mathrm{MLE}})=-n\log\det\hat{\bm{\Psi}}_{\mathrm{MLE}}+\textup{tr}(\hat{\bm{\Psi}}_{\mathrm{MLE}}\textbf{S})$. We select the $c$ that results in the smallest BIC.

6 Posterior Consistency

In this section, we show that in the fixed dimensional setting (i.e. $p$ is fixed), the rank-based posterior distribution of $\bm{\Psi}$ is consistent at its true value $\bm{\Psi}_{0}$ for almost all $\bm{\Psi}_{0}$ with respect to the Lebesgue measure under the only assumption that the prior distribution for $\bm{\Psi}$ has positive density on the space of inverse correlation matrices, i.e., the collection of all $p\times p$ positive definite matrices $\bm{\Psi}$ such that all diagonal elements of $\bm{\Psi}^{-1}$ are $1$.

Theorem 6.1**.**

Assume that $\bm{\Psi}$ has prior density $\pi(\bm{\Psi})>0$ a.e. over the space of positive definite matrices, with respect to the Lebesgue measure $\nu$ and that the dimension $p$ is fixed. Then for $\bm{\Psi}_{0}$ a.e. $[\nu]$, and for any neighborhood $\mathcal{U}_{0}$ of $\bm{\Psi}_{0}$, we have that

[TABLE]

where $[P^{\infty}_{\bm{\Psi},H}]$ denotes the joint distribution of all $\bm{X}$’s and ranks $\bm{R}$ with $\bm{\Psi}_{0}$ as the true value of $\bm{\Psi}$ and $\bm{H}$ denotes the underlying normality restoring transformations.

Proof.

The proof is based on an application of Doob’s Theorem [Ghosal and van der Vaart, 2017], Section 6.2. Doob’s theorem is a very general posterior consistency result, which only requires that in the joint distribution of all parameters and observables, the parameter can be a.s. written as a function of all observables of all stages, and then concludes that posterior consistency holds for almost all parameters with respect to the prior distribution. Here observations are the rank information $\bm{R}_{n}=(\!(R_{n,ik})\!)$ for each variable, where $R_{n,ik}$ stands for the rank of the $i$th observation in the $k$th component, $i=1,\ldots,n$, $k=1,\ldots,p$, $n=1,2,\ldots$. We follow some arguments given in Gu and Ghosal [2009]. Let $U_{ik}=F_{k}(X_{ik})$, the “population quantile” of the $i$th observation regarding the $k$th variable, $i=1,\ldots,n$, $k=1,\ldots,p$. Note that $Y_{ik}=\Phi^{-1}(U_{ik})$. By Theorem a on page 157 of Hájek and Šidák [1967], for any $k=1,\ldots,p$,

[TABLE]

As such, $U_{ik}$ is an in-probability limit of $\mathcal{F}_{n}$-measurable random variables, where $\mathcal{F}_{n}$ is the $\sigma$-field generated by $\{R_{n,ik}:i=1,\ldots,n,\,k=1,\ldots,p\}$. Thus, for any $k=1,\ldots,p$,

[TABLE]

and hence, $U_{ik}$ is an $\mathcal{F}_{\infty}$-measurable random variable, where $\mathcal{F}_{\infty}=\sigma\langle\cup_{1}^{\infty}\mathcal{F}_{n}\rangle$, and so is $Y_{ik}=\Phi^{-1}(U_{ik})$. Therefore it suffices to show that $\bm{\Psi}$ can be written as the almost sure limit of a sequence of functions of $\{Y_{ik}:i=1,\ldots,n,\,k=1,\ldots,p\}$.

Let $C_{jk}$ stand for the $(j,k)$th element of $\bm{C}=\bm{\Psi}^{-1}$. Then clearly $C_{jk}=n^{-1}\lim_{n\to\infty}\sum_{i=1}^{n}Y_{ij}Y_{ik}$ almost surely. Thus $\bm{C}$, and hence $\bm{\Psi}$ is expressible as an almost sure limit of $\{Y_{ik}:i=1,\ldots,n,\,k=1,\ldots,p\}$. Thus, by Doob’s theorem, the consistency of the posterior (4) at $\bm{\Psi}_{0}$ holds a.e. $[\pi]$. However, as the prior density is positive throughout the parameter space, it also follows that the posterior (4) at $\bm{\Psi}_{0}$ holds a.e. $[\nu]$. ∎

Usually, the main criticism against a posterior consistency result obtained by applying Doob’s theorem is that the exceptional set where consistency may fail may be “large” since it only needs to be null with respect to the prior, which is somewhat arbitrary. However, in the present application, since the parameter space for the parameter of interest $\bm{\Psi}$ is finite dimensional, where we have the Lebesgue measure as a benchmark measure, the exceptional set of points where the posterior may be inconsistent is characterized as Lebesgue null, which can be regarded as “small”. It is important to note that the normality restoring transformations are taken to be fixed, and no prior is assigned on them. Since the underlying procedure is unaffected by the transformations, the fact that these transformations are unknown does not matter. Note that the fixed dimensionality of the variables is essential since the posterior consistency result is applicable only if the parameter space is fixed.

7 Simulation Results

We conduct a simulation study to compare the performance of the proposed Bayesian rank-likelihood method, the Bayesian GCGM [Mohammadi et al., 2017], the empirical method in a nonparanormal graphical model [Liu et al., 2009], the Bayesian method in a nonparanormal graphical model [Mulgrave and Ghosal, 2018b], and the empirical method based on GCGM [Liu et al., 2012]. Both the proposed Bayesian rank-likelihood method, indicated as Rank-Likelihood, and the Bayesian method of Mulgrave and Ghosal [2018b], indicated as B-splines, use a horseshoe prior on the Cholesky decomposition of the precision matrix and MCMC estimation. The B-splines method uses a random series B-splines prior to estimate the transformation functions. The Bayesian method based on GCGM, indicated as Bayesian Copula, uses the rank-likelihood to transform the random variables and puts a G-Wishart prior on the inverse correlation matrix, a uniform prior on the graph, and estimates the sparse matrix using a birth-and-death MCMC [Mohammadi and Wit, 2015]. The empirical method in a nonparanormal graphical model, indicated as Truncation, uses a truncated marginal empirical distribution function to transform the variables and the graphical lasso to estimate the sparse precision matrix. The empirical method based on GCGM, indicated as SKEPTIC, uses Spearman’s rho to transform the variables and estimates the sparse precision matrix with the graphical lasso.

The random variables, $Y_{1},\ldots,Y_{p}$, are simulated from a multivariate normal distribution such that $Y_{i1},\ldots,Y_{ip}\stackrel{{\scriptstyle\mathrm{i.i.d.}}}{{\sim}}\mathrm{N}(\bm{\mu},\bm{\Omega}^{-1})$ for $i=1,\ldots,n$. The means $\bm{\mu}$ are selected from an equally spaced grid between 0 and 5 with length $p$. We consider nine different combinations of $n,p,$ and sparsity for $\bm{\Omega}$:

$p=25,\;n=50$, AR(4) model;

• 2.

$p=50$, $n=100$, AR(4) model;

• 3.

$p=100$, $n=500$, AR(4) model;

• 4.

$p=25,\;n=50$, AR(1) model;

• 5.

$p=50$, $n=100$, AR(1) model;

• 6.

$p=100$, $n=500$, AR(1) model;

• 7.

$p=25$, $n=50$, sparsity = $10\%$ non-zero entries in the off-diagonals;

• 8.

$p=50$, $n=100$, sparsity = $5\%$ non-zero entries in the off-diagonals;

• 9.

$p=100$, $n=500$, sparsity = $2\%$ non-zero entries in the off-diagonals;

where the AR(4), AR(1), and star models are described by the relations

AR(4) model: $\omega_{i,i}=1,\;\omega_{i,i-1}=\omega_{i-1,i}=0.2,\;\omega_{i,i-2}=\omega_{i-2,i}=0.2,\;\omega_{i,i-3}=\omega_{i-3,i}=0.2,\;\omega_{i,i-4}=\omega_{i-4,i}=0.1$;

• 2.

AR(1) model: $\omega_{1,1}=1.9608,\;\omega_{i,i-1}=\omega_{i-1,i}=-1.3725$, $\omega_{p,p}=1.9608$

The percent sparsity levels for $\bm{\Omega}$ are computed using lower triangular matrices that have diagonal entries normally distributed with mean 1 and standard deviation 0.1, and non-zero off-diagonal entries normally distributed with mean 0 and standard deviation 1.

The observed variables $\bm{X}=(X_{1},\ldots,X_{p})$ are constructed from the simulated variables $Y_{1},\ldots,Y_{p}$. The functions used to construct the observed variables were four cumulative distribution functions (c.d.f.s): asymmetric Laplace, extreme value, logistic, and stable. We could choose any values of the parameters for the c.d.f.s, but instead of selecting them ourselves, we automatically choose the values of the parameters to be the maximum likelihood estimates with the mle function in MATLAB. The values of the parameters for each of the c.d.f.s are the maximum likelihood estimates for the parameters of the corresponding distributions (asymmetric Laplace, extreme value, logistic, and stable), using the variables $Y_{1},\ldots,Y_{p}$.

We follow the procedure in Mulgrave and Ghosal [2018a] to estimate the transformation functions for the B-splines method. The hyperparameters for the normal prior are chosen to be $\nu=1,\tau=1,$ and $\sigma^{2}=1$. To choose the number of basis functions, we use the Aikaike Information Criterion. Samples from the truncated multivariate normal posterior distributions for the B-spline coefficients are obtained using the exact Hamiltonian Monte Carlo (exact HMC) algorithm [Pakman and Paninski, 2014]. After finding the initial coefficient values $\bm{\theta}_{d}$, we construct initial values for $Y_{d,\textup{initial}}=\sum_{j=1}^{J}\theta_{dj,\textup{initial}}B_{j}(X_{d})$ using the observed variables. These initial values $\bm{Y}_{\textup{initial}}$ are used to find initial values for $\bm{\Sigma},\bm{\mu}$, and $\bm{\Omega}$ for the algorithm.

For the Rank-Likelihood method, we initialize the algorithm using the ranks of the observed variables. First, we create a $p\times n$ matrix of ranks, $\bm{R}=\mbox{rank}(\bm{X})$. Then, we divide the columns by $n+1$ and transform each entry by the standard normal quantile function using the inverse standard normal c.d.f., such that $\bm{Y}_{\textup{initial}}=\bm{\Phi}^{-1}(\bm{R}/(n+1))^{\prime}$. Finally, we take the transpose to obtain an $n\times p$ matrix. We obtain an initial value for the inverse correlation matrix using $\bm{\Psi}_{\textup{initial}}=(\mbox{corr}(\bm{Y}_{\textup{initial}}))^{-1},$ where corr stands for the correlation coefficients of $\bm{Y}_{\textup{initial}}.$

For both the Rank-Likelihood and the B-splines methods, the hyperparameters for the horseshoe prior are initialized with ones. To impose the sparsity constraint on the Cholesky decomposition of the matrices for the B-splines and the Rank-Likelihood methods, we consider four values for tuning: $c\in\{0.1,1,10,100\}$. We select the graphical model with the value of $c$ having the lowest BIC. The 0-1 loss procedure (4.1) is used to threshold the matrices for the B-splines and the Rank-Likelihood methods and construct the corresponding edge matrices. The codes for these methods are written in MATLAB.

For the simulation study, we run 100 replications for each of the nine combinations and assess structure learning and parameter estimation for each replication. We collect $10000$ MCMC samples for inference after discarding a burn-in of $5000$ and we do not apply thinning. The Bayesian Copula method is implemented in the R package BDgraph [Mohammadi and Wit, 2017, 2019] using the option “gcgm”. Bayesian model averaging is used for inverse correlation matrix and graph selection. The default option in the BDgraph package selects the graph formed by links having estimated posterior probabilities greater than $0.5$. The Truncation and SKEPTIC methods are implemented in the R package huge [Zhao et al., 2015] using the “truncation” and “skeptic” options respectively. For both empirical methods, the graphical lasso method is used for the graph estimation based on transformed variables and the default lossless screening method [Witten et al., 2011, Mazumder and Hastie, 2012a] is applied. A sequence of 100 regularization parameters, $\lambda$, is generated starting from $\lambda_{\mbox{max}}$ to $0.01*\lambda_{\mbox{max}}$, in the log-scale. We define $\lambda_{\mbox{max}}=\max(\max(\bm{M}-\mbox{diag}(p)),-\min(\bm{M}-\mbox{diag}(p)))$, where $\bm{M}$ is the correlation matrix from the SKEPTIC method and $\bm{M}$ is the matrix constructed for the Truncation method after scaling the transformed variables and converting it to a correlation matrix. Note that $p$ is the dimension. The method used to select the graphical model along the regularization path is Generalized Stability Approach to Regularization Selection [Müller et al., 2016], or G-StARS, implemented using the R package pulsar. We use upper and lower bounds to reduce the computational burden and the number of random subsamples taken for graph re-estimation was 100. All codes are provided in the Supplementary Material.

7.1 Performance Assessment

For the Rank-Likelihood method, we find the Bayes estimate of the inverse correlation matrix $\hat{\bm{\Upsilon}}=\mathds{E}(\bm{\Psi}|\bm{Y})$ and average it over MCMC iterations. For the B-splines method, the Bayes estimate of the precision matrix is $\hat{\bm{\Upsilon}}=\mathds{E}(\bm{\Omega}|\bm{X})$, using the MCMC samples. The median probability model [Berger and Barbieri, 2004] is used to find the Bayes estimate of the edge matrix for both the Rank-Likelihood and B-splines methods. We find the estimated edge matrix by first using the 0-1 loss procedure (4.1) to threshold the MCMC inverse correlation and precision matrices, and then we take the mean of the thresholded matrices. If the off-diagonal element of the mean is greater than $0.5$, the element is registered as an edge; else it is registered as a no-edge.

We compute the specificity (SP), sensitivity (SE), and Matthews Correlation Coefficient (MCC), previously used for assessing the accuracy of classification procedures [Baldi et al., 2000], to assess the performance of the graphical structure learning. They are defined as follows:

[TABLE]

where TP stands for true positives, TN stands for true negatives, FP stands for false positives, and FN stands for false negatives. Specificity and sensitivity values are between 0 and 1, where 1 is the best value. MCC values are between $-1$ and 1, and 1 is the best value.

We also assess the strength of parameter estimation. We consider the scaled $L_{1}$-loss function, the average absolute distance. The scaled $L_{1}$-loss is defined as

[TABLE]

For the Rank-Likelihood, SKEPTIC, and Bayesian Copula methods, $\hat{\bm{\Upsilon}}$ is the estimated inverse correlation matrix and $\bm{\Upsilon}_{\textup{true}}$ is the true inverse correlation matrix. For the Truncation and B-splines method, $\hat{\bm{\Upsilon}}$ is the estimated inverse covariance matrix and $\bm{\Omega}_{\textup{true}}$ is the true inverse covariance matrix. The results are presented in Figures 1–4.

The Rank-Likelihood method performs consistently better than the B-splines method in terms of structure learning and parameter estimation. In particular, the Rank-Likelihood method appears to be more sensitive and specific to signals than the B-splines method. In addition, the scaled $L_{1}$-loss of the Rank-Likelihood method is significantly better than the scaled $L_{1}$-loss of the B-splines method. The Bayesian Copula method performs similar to or better than SKEPTIC and Truncation methods in terms of structure learning for all models considered. The Bayesian Copula method performs similar to the SKEPTIC and Truncation methods with regard to parameter estimation. The proposed Rank-Likelihood method performs the best at structure learning for all models considered except for the AR(4) model at dimension $p=100$, at which the Bayesian Copula model performs the best. For parameter estimation, the proposed Rank-Likelihood method generally outperforms all competing methods.

Thus overall, compared to competing methods, the Rank-Likelihood method performs nearly the same or better for structure learning and parameter estimation for all models excluding the AR(4) model at dimension $p=100$. However, the Rank-Likelihood model has good performance when considering structure learning and parameter estimation together, compared to the competing models.

8 Real Data Application

We demonstrate the methods on a gene expression data set originally referenced in Stranger et al. [2007] with Gene Expression Omnibus database Series Accession Number GSE6536 19 and its funding is supported in part by the US National Institutes of Health ENDGAME. Data are collected to measure the gene expression in B-lymphocyte cells from inhabitants in Utah with European ancestry. The interest is on the single nucleotide polymorphisms that are found in the 5’ untranslated region of messenger RNA with a minor allele frequency $\geq 0.1$. Following Bhadra and Mallick [2013], of the 47,293 total available probes, we considered the 100 most variable probes that correspond to different Illumina TargetID transcripts. The data for these 100 transcripts are available in the R package BDgraph [Mohammadi and Wit, 2017, 2019]. The data consist of $n=60$ unrelated individuals and $p=100$ transcripts. The variables in the data are continuous but do not appear Gaussian. A Bayesian estimate based on a Gaussian graphical model using a spike-and-slab type prior constructed by Bhadra and Mallick [2013] detected 55 edges.

To construct the graph using our method, we convert the original values to be between 0 and 1 using the affine transform $({x-\min(x_{i})})/({\max(x_{i})-\min(x_{i})})$. We use the identity matrix as the initial matrix for the covariance and inverse covariance matrices for the Rank-Likelihood and B-splines methods. The Rank-Likelihood method results in 252 edges and the B-splines method results in 99 edges. Convergence of the Rank-Likelihood method can be obtained in about 28 minutes and for the B-splines method in about 60 minutes for a given $c$ for these data on a laptop computer with an Intel i7 processor and 24 GB of RAM. Using the same set-up as in the simulation study, the SKEPTIC method results in no edges and the Truncation method results in 363 edges. The Bayesian Copula method results in 834 edges. The proposed Rank-Likelihood and B-splines methods result in the sparsest models, with the B-splines method as the most sparse model. The graphs of the proposed methods are shown in Figure 1. The graphs of the Bayesian Copula and truncation methods are shown in Figure 2. Since the SKEPTIC method resulted in no edges, it is not included in the comparison. Plots are made with the circularGraph function in MATLAB.

SUPPLEMENTAL MATERIALS

GitHub Repository:

The code for the methods described in this paper can be found in the following GitHub repository: https://github.com/jnj2102/RankLikelihood.