# Bayesian Estimation of Gaussian Graphical Models with Predictive Covariance
Selection

Donald R.Williams  Juho Piironen  Aki Vehtari  Philippe Rast  Department of
Psychology, University of California, Davis, USA  Department of Computer
Science, Aalto University, Espoo, Finland

###### Abstract

Gaussian graphical models are used for determining conditional relationships
between variables. This is accomplished by identifying off-diagonal elements
in the inverse-covariance matrix that are non-zero. When the ratio of
variables ($p$) to observations ($n$) approaches one, the maximum likelihood
estimator of the covariance matrix becomes unstable and requires shrinkage
estimation. Whereas several classical (frequentist) methods have been
introduced to address this issue, fully Bayesian methods remain relatively
uncommon in practice and methodological literatures. Here we introduce a
Bayesian method for estimating sparse matrices, in which conditional
relationships are determined with projection predictive selection. With this
method, that uses Kullback-Leibler divergence and cross-validation for
neighborhood selection, we reconstruct the inverse-covariance matrix in both
low and high-dimensional settings. Through simulation and applied examples, we
characterized performance compared to several Bayesian methods and the
graphical lasso, in addition to TIGER that similarly estimates the inverse-
covariance matrix with regression. Our results demonstrate that projection
predictive selection not only has superior performance compared to selecting
the most probable model and Bayesian model averaging, particularly for high-
dimensional data, but also compared to the the Bayesian and classical glasso
methods. Further, we show that estimating the inverse-covariance matrix with
multiple regression is often more accurate, with respect to various loss
functions, and efficient than direct estimation. In low-dimensional settings,
we demonstrate that projection predictive selection also provides competitive
performance. We have implemented the projection predictive method for
covariance selection in the R package GGMprojpred .

Bayesian,

Gaussian graphical models,

Neighborhood selection,

Projection predictive selection,

KL-divergence,

Cross-validation,

Sparsity,

###### keywords:

, , , and

##  1 Introduction

Gaussian graphical models (GGMs) are used for covariance selection, in which
conditional dependencies between random variables are characterized (Dempster,
1972; Peng et al., 2009). This is accomplished by identifying off-diagonal
elements in the inverse-covariance matrix (i.e., precision matrix) that are
non-zero. When these covariances are standardized and the sign reversed, they
correspond to partial correlations. Assuming multivariate normality, partial
correlations imply pairwise conditional dependence (i.e., direct effects)
while controlling for all other variables included in the model (Baba et al.,
2004; Baba and Sibuya, 2005). This stands in contrast to marginal
correlations, where the associations include indirect and overlapping effects,
thereby limiting their use in assessing unique relationships among a set of
variables (Schäfer and Strimmer, 2005a; Roverato and Castelo, 2017). That GGMs
characterize conditional relationships among random variables has led to
extensive use across the sciences, in both methodological and applied
contexts. For example, GGMs are commonly used to characterize genetic co-
expression networks (Wang, 2016), functional connectivity between brain
regions (Das et al., 2017), and more recently to understand chronic mental
illnesses (McNally et al., 2015).

The partial correlation matrix $\tilde{\textbf{{P}}}=(\tilde{\rho}_{ij})$ can
be obtained by inverting the sample covariance matrix $\Sigma$, where
$\Sigma^{-1}$ is the corresponding precision matrix $\Omega$. Here, with the
covariances $\hat{\omega}_{ij}$ and precisions $\omega_{ii,jj}$ $\in$
$\Omega$, the partial correlations $\rho_{ij}$ are computed as

$\tilde{\rho}_{ij}=\frac{-\omega_{ij}}{\sqrt{\omega_{ii}\omega_{jj}}},\hskip
4.26773pt\text{where}\hskip 4.26773pt\Omega=\begin{bmatrix}\omega_{ii}&\\\
\vdots&\ddots\\\ \omega_{ij}&\cdots&\omega_{jj}\\\ \end{bmatrix}=\Sigma^{-1}.$
(1.1)

In ideal situations, for example when the ratio of variables $p$ to
observations $n$ is sufficiently small ($p/n\ll 1$), the customary maximum
likelihood (ML)
$\frac{1}{N}\sum_{i=1}^{n}(X_{i}-\bar{X})(X_{i}-\bar{X})^{\top}$ can provide
an accurate estimate for $\Sigma$. Hence, in this case, inversion to obtain
$\Omega$ is straightforward. However, in many practical applications, there
are often many variables relative to the number of observations. In these
situations, ML becomes unstable and less accurate as $p\rightarrow n$ which is
further magnified when inverting $\Sigma$ (Won et al., 2009; Ledoit and Wolf,
2004a). This does not only influence accuracy of estimation, but also the
ability to correctly identify non-zero elements. For example, while confidence
intervals for the precision matrix have nominal coverage when $p/n\ll 1$
(Drton and Perlman, 2004; Williams and Rast, 2018), they were shown to have
less than nominal probabilities even when the $p/n$ ratio was below 0.50
(Janková and van de Geer, 2015). These errors are a result of increased
variability in the eigenvalues $\lambda$ (Kuismin and Sillanpää, 2016), from
which the determinant of the matrix is computed:

$\text{det}(\Sigma)=\prod_{i=1}^{p}\lambda_{i},\hskip 1.42271pti\hskip
1.42271pt\in\\{1,...,p\\}.$ (1.2)

Further, since the number of non-zero eigenvalues is min($n$, $p$) (Kuismin
and Sillanpää, 2017), it is clear from Equation 1.2 that $\Sigma$ is singular
det($\Sigma$) = 0 when $n<p$. As a result, in this situation, it is not
possible to compute the ML estimate for $\Omega$ (Hartlap et al., 2007).

To address these issues, in the classical literature (i.e., frequentist), two
primary lines of research have emerged to estimate $\Omega$. The first
surrounds accurate estimation with convex combination estimators (Kuismin and
Sillanpää, 2017), and not necessarily determining non-zero relationships
(Ledoit and Wolf, 2004a, b). Of course, once obtained, a decision rule can be
adopted for identifying non-zero relationships (Huang and Wand, 2013; Kubokawa
and Srivastava, 2008; Schäfer and Strimmer, 2005a). Although this typically
requires some form of null hypothesis testing for the partial correlations,
for example using the approach introduced in Schäfer and Strimmer (2005b),
recently Kuismin and Sillanpää (2016) applied a stepwise strategy directly to
$\omega_{ij}$ $\in$ $\Omega$. The most popular approaches (Wang, 2012; Yuan
and Wang, 2013), however, optimize the penalized maximum likelihood with
respect to a $\ell$-1 norm constraint. This approach not only shrinks
$\omega_{ij}$ towards zero, but can also set the elements to exactly zero
($\omega_{ij}=0$). There are two strategies for utilizing $\ell$-1
regularization (Krämer et al., 2009; Friedman et al., 2008), the first of
which estimates $\Omega$ directly such as the popular graphical lasso (glasso)
method (Friedman et al., 2008), whereas the second avoids the issues of
directly estimating $\Omega$ altogether. Instead, the latter uses a column-
wise multiple regression strategy. Here, assuming the input $n\times p$ matrix
is standardized, each variable $p$ is regressed against the remaining $p-1$
variables. The partial correlation then follows

$\rho_{ij}=\text{sign}(\beta_{i}^{(j)})\sqrt{\beta_{i}^{(j)}\beta_{j}^{(i)}},$
(1.3)

where $\beta_{i}^{(j)}$ is the penalized coefficient for predictor $X_{j}$ and
response $X_{i}$. Several regression approaches, including ridge (Krämer et
al., 2009), partial least squares (Tenenhaus et al., 2008), and lasso have
been used in this context (Meinshausen and Bühlmann, 2006). In this work, we
introduce a Bayesian method that also utilizes a regression based approach for
both accurate estimation and determining conditional relationships
($\omega_{ij}\neq 0$) in high dimensional settings $n<p$.

Our primary focus is on fully Bayesian estimation of GGMs, which remain
relatively uncommon in practice compared to classical methods. While several
classical regression strategies have been employed, only recently was a
Bayesian lasso regression characterized (Li and Zhang, 2017). Perhaps the most
noteworthy Bayesian approaches makes use of the $G$-Wishart distribution,
which is conjugate for $\Omega$ whose covariances are constrained to be zero
(Atay-Kayis and Massam, 2005; Dobra et al., 2010). The earliest approaches
utilized reversible jump MCMC (Green, 1995), for example in Wang (2012) and
Lenkoski (2013), whereas the most recent work use a birth and death MCMC
algorithm (Mohammadi and Wit, 2015b). Both provide a familiar framework for
Bayesian model selection, in that covariance selection uses an indicator
function analogous to customary stochastic search algorithms for general
linear models (George and McCulloch, 1993; O’Hara and Sillanpää, 2009). That
is, the underling conditional independence structure can be obtained with
either Bayesian model averaging or maximum a posterior probability.
Alternative Bayesian approaches, on the other hand, determine non-zero
relationships with credible intervals. Recently, block-wise sampling
algorithms have been developed to estimate Bayesian versions of the graphical
lasso (Wang, 2012; Khondker et al., 2013), in addition to an adaptive version
described in Wang (2012).

Although these Bayesian methods have compared favorably to classical methods,
there are several important limitations that we address in this work. First,
as noted above, the ML estimate and corresponding intervals are adequate when
$p/n\ll 1$. However, most of the simulation studies have almost exclusively
characterized performance in low-dimensional settings (Mohammadi and Wit,
2015a). To our knowledge, performance of $G$-Wishart distributions have not
been characterized in high-dimensional settings. Second, since Bayesian
approaches provide a measure of parameter uncertainty via the posterior
distribution, parameters are never actually zero and covariance selection must
be achieved with an explicit decision rule. In Khondker et al. (2013), for
example, they used used credible intervals that ranged from 10-50 %, none of
which had a clear justification. Further, in high dimensional settings (van
der Pas et al., 2016), credible intervals are unlikely to have nominal
frequentist properties that are known to translate into performance measures
for estimating GGMs (Williams and Rast, 2018, i.e., specificity = 1 -
$\alpha$).

In the present paper, we introduce a novel Bayesian method for estimating
Gaussian graphical models. Our method does not rely on posterior model
probabilities or credible interval width, but we adopt a decision theoretic
perspective. To determine the conditional independence structure, we use
projection predictive selection, described in more detail in section 3, that
allows for variable exclusion ($\omega_{ij}=0$) based on predictive utility
(Goutis, 1998; Dupuis and Robert, 2003). In the context of the general linear
model, this approach was shown to outperform information criteria, cross-
validation scores, or selecting the most probable model (Piironen and Vehtari,
2017a). We also focus on high dimensional settings ($n<p$), where the ML
estimate does not exist, but also examine consistency in typical asymptotic
settings (fixed $p$; $n\rightarrow\infty$). This not only shows the generality
of the proposed method, but also provides a contrast to the known
inconsistencies of glasso (Leppä-aho et al., 2017; Williams et al., 2018;
Kuismin and Sillanpää, 2016).

The remainder of the manuscript is organized as follows. Thus far we have been
exclusively discussing covariance selection in terms of $\Omega$, so we first
introduce the customary notation and terminology for GGMs (section 2), and
then describe the relationship between column-wise multiple regression and the
elements of $\Omega$. Second, we then outline the projection predictive method
in the context of both low-dimensional and high-dimensional GGMs. For the
former, posterior samples are generated with the Bayesian bootstrap, whereas
for the latter we use the horseshoe prior distribution (section 3). In section
4, we use numerical experiments to evaluate our methods performance compared
to $G$-Wishart BMA and MAP estimation, in addition to the Bayesian graphical
lasso. As a point of reference, we also include the classical glasso and an
alternative $\ell_{1}$-regularization method that also constructs $\Omega$
with multiple regression. In the following section, we use each method to
estimate a known genetic regulatory network. We end by discussing limitations
and future directions.

##  2 Gaussian Graphical Models

We now introduce notation and terminology specific to graphical models.
Depending on the field, undirected graphical models refer to covariance
selection models and random Markov fields. Here we adopted the term Gaussian
graphical model, because it is common in the Bayesian literature and provides
an informative description of the method. For example, consider a
p-dimensional random variable $X$ that follows a multivariate normal
distribution

$X=\\{X_{1},...,X_{p}\\}\sim\mathcal{N}(\mu,\Sigma),$ (2.1)

where $\Sigma$ is an unknown, but positive definite covariance matrix and
$\mu$ is a vector of means. The inverse of the covariance matrix is
$\Omega=\Sigma^{-1}$, where the $i$th row and $j$th column is denoted by
$\Omega_{ij}$. Without loss of generality, we assume that all variables have
been standardized to have mean zero and variance one:

$0=\\{\mu_{1},...,\mu_{p}\\}^{\mathsf{T}},\hskip 4.26773pt\text{and}\hskip
4.26773pt(\Sigma_{ii})=1.$ (2.2)

Expanded to multiple observations, with $n$ independent samples, let X be the
$n\times p$ data matrix. The graph is then denoted by $\mathcal{G}=(V,E)$.
Here $V$ = $\\{1,...,p\\}$ is the node set and $E\subset V\times V$ is the
edge set. Thus, while $V$ represents the total number of columns in X, each
dimension of Equation 2.1 is denoted by $\textit{X}_{i}$. The edge set for
$\mathcal{G}$ contains nodes $(X_{i},X_{j})$ that share a conditional
relationship $X_{i}\not\Vbar X_{j}|\boldsymbol{X}_{V-i,-j}$. In contrast,
conditionally independent nodes $X_{i}\Vbar X_{j}|\boldsymbol{X}_{V-i,-j}$ are
not included in $E$ and correspond to the zero elements within $\Omega$
($\omega_{ij}=0$). The maximum number of edges possible is $\frac{V(V-1)}{2}$,
which corresponds to the number of non-zero off-diagonal elements
$\omega_{ij}\in$ $\Omega$.

###  2.1 Regression-Based Interpretation

We estimate $\mathcal{G}$ by fitting node-wise regression models, a procedure
known as neighborhood selection. The random variables can be partitioned into
two groups $X=(X_{i},X_{-i})$, where
$X_{-i}=(X_{1},\ldots,X_{i-1},X_{i+1},\ldots,X_{p})$. Here, assuming $X$
follows a multivariate normal distribution (Equation 1.3), the conditional
distribution of $X_{i}$ given $X_{-i}$ is also normally distributed Anderson
(2003), hence,

$X_{i}|X_{-i}\sim\mathcal{N}(\mu_{i}+\Sigma_{i,-i}\Sigma_{-i,-i}^{-1}(X_{-i}-\mu_{-i}),\Sigma_{ii}-\Sigma_{i,-i}\Sigma_{-i,-i}^{-1}\Sigma_{-i,i}).$
(2.3)

The equivalent regression equation is denoted as

$X_{i}=\alpha_{i}+X_{-i}^{\top}\beta_{(i)}+\varepsilon_{i},$ (2.4)

where $\alpha=\mu_{i}-\Sigma_{i,-i}\Sigma_{-i,-i}^{-1}\mu_{-i}$,
$\beta_{(i)}=\Sigma_{-i,-i}^{-1}\Sigma_{-i,i}$ is a $p-1$ dimensional vector,
and
$\varepsilon_{i}\sim\mathcal{N}(0,\Sigma_{ii}-\Sigma_{i,-i}\Sigma_{-i,-i}^{-1}\Sigma{-i,i})$
is independent of $X_{i}$ (Yuan, 2010; Li et al., 2017). The inverse,
resulting in the $i$th column of $\Omega$, follows

$\displaystyle\omega_{ii}$
$\displaystyle=(\Sigma_{ii}-\Sigma_{i,-i}\Sigma_{-i,-i}^{-1}\Sigma_{-i,i})^{-1},$
(2.5) $\displaystyle=\text{Var}(\varepsilon_{k})^{-1},$
$\displaystyle\Omega_{-i,i}$
$\displaystyle=-(\Sigma_{ii}-\Sigma_{i,-i}\Sigma_{-i,-i}^{-1}\Sigma_{-i,i})^{-1}\Sigma_{-i,-i}^{-1}\Sigma_{-i,i},$
(2.6) $\displaystyle=-\text{Var}(\varepsilon_{k})^{-1}\beta_{(i)},$

where the resulting precision matrix is

$\Omega=\begin{bmatrix}\frac{1}{\text{Var}(\epsilon_{k})}_{11}&\frac{-\beta_{12}}{\text{Var}(\epsilon_{k})}_{11}&\cdots&\frac{-\beta_{1p}}{\text{Var}(\epsilon_{k})}_{11}\\\
\frac{-\beta_{21}}{\text{Var}(\epsilon_{k})}_{22}&\frac{1}{\text{Var}(\epsilon_{k})}_{22}&\cdots&\frac{-\beta_{2p}}{\text{Var}(\epsilon_{k})}_{22}\\\
\vdots&\vdots&\ddots&\vdots\\\
\frac{-\beta_{p1}}{\text{Var}(\epsilon_{k})}_{pp}&\frac{-\beta_{p2}}{\text{Var}(\epsilon_{k})}_{pp}&\cdots&\frac{1}{\text{Var}(\epsilon_{k})}_{pp}\\\
\end{bmatrix}.$ (2.7)

In the GGM literature, it has been shown that this approach approximates
$\Omega$. However, when using least squares, it should be noted that the
relationship is exact only when $\Sigma$ can be inverted ($p<n$). Detailed
analytic proofs, that we use to derive $\Omega$ in high-dimensional settings
(section 3.3), are provided in Stephens (1998). This is a fundamental problem
and in the following section we present a solution for cases where $p>n$.

##  3 Projection Predictive Selection

Projective model selection is a model simplification technique and
fundamentally a two-step procedure that provides a decision theoretically
correct inference after selection. That is, the full posterior is projected to
the restricted subspace instead of forming the usual posterior given
constraints. Importantly, the projection method is general in that selection
only requires samples from the posterior, thus allowing for flexible model
specification.

###  3.1 Model Fitting

The first step requires the construction of the encompassing reference model
$M_{*}$, defined as

$\mathbf{y}=\boldsymbol{X}\beta\hskip 1.42271pt+\hskip
1.42271pt\varepsilon,\hskip
1.42271pt\varepsilon\sim\mathcal{N}(0,\sigma^{2}),$ (3.1)

where X is a $n\times p-1$ matrix that includes a subset of nodes $X_{-i}$
from $V$ $\\{1,...,p\\}$, $\beta$ is a $p-1$ dimensional vector including
$\beta_{j}^{(i)}=\\{j,...,p-1\\}$, and y is the $i$th node to be predicted.
This model best captures our assumptions and uncertainties related to the
problem, including the choice of prior distributions for the task at hand
(e.g., the dimensions of the data).

Building upon Equation 3.1, for the purpose of regularization, we use the
horseshoe prior distribution

$\displaystyle\beta_{j}|\lambda_{j},\tau$
$\displaystyle\sim\mathcal{N}(0,\lambda_{j}^{2}\tau^{2}),\hskip
4.26773ptj\in\\{1,...,p-1\\},$ (3.2) $\displaystyle\lambda_{j}$
$\displaystyle\sim\text{C}^{+}(0,1),$ $\displaystyle\tau$
$\displaystyle\sim\text{C}^{+}(0,\tau_{0}),$

where $\text{C}^{+}$ is a half-Cauchy distribution for the local and global
hyperparameters denoted with $\lambda$ and $\tau$, respectively. The latter
provides shrinkage to all estimates, whereas the local hyperparameter allows
stronger signals to avoid shrinkage, thus making this choice of prior ideal
for identifying non-zero relationships. It is possible to incorporate prior
expectations in regards to sparsity, denoted $\tau_{0}$, that serves as the
scale for the prior on $\tau$. Described in Piironen and Vehtari (2017a), the
value for $\tau_{0}$ can be computed as

$\tau_{0}=\frac{p_{0}}{D-p_{0}}\frac{\sigma}{\sqrt[]{N}},$ (3.3)

where $p_{0}$ denotes the number of expected edges and $D$ is the number of
variables included in the model ($p-1$). Since priors can have diminishing
returns as $n\rightarrow\infty$, the term $\frac{\sigma}{\sqrt[]{N}}$ allows
for scaling with the data. In practice, while it is recommended to include
relevant information when available, Piironen and Vehtari (2017b) demonstrate
robustness of the results with respect to the specific value. Here we assume a
generic value ($\tau_{0}$ = 1)

###  3.2 Neighborhood Selection

In the second stage, for each of the $p$ regression models, we then determine
which nodes can be removed, or set to zero ($\hat{\beta}_{j}^{(i)}=0$),
without introducing considerable predictive loss compared to the original
reference model $M_{*}$ (Equation 3.1). In the predictive framework, the model
parameters of the reduced models, herein referred to as submodels including
nodes $x$, are determined so that the resulting predictive distribution
remains close to that of $M_{*}$ (Vehtari and Ojanen, 2012). A theoretically
justified and computationally convenient way to do this is by minimizing the
discrepancy between the predictive distributions conditional on the parameter
values, where the discrepancy $\delta$ is defined in terms of Kullback-Leibler
(KL) divergence averaged over the empirical distribution for $X$ (Goutis,
1998; Dupuis and Robert, 2003)

$\displaystyle\theta_{\perp}$
$\displaystyle=\arg\min_{\theta^{\prime}\in\Theta_{\perp}}\delta(\theta,\theta^{\prime})$
$\displaystyle\triangleq\arg\min_{\theta^{\prime}\in\Theta_{\perp}}\frac{1}{n}\sum_{i=1}^{n}\mathrm{KL}{\left(p(\tilde{y}\,|\,\theta,X_{i},M_{*})\,\|\,p(\tilde{y}\,|\,\theta^{\prime},x_{i},M_{\perp})\right)}.$
(3.4)

This defines a parameter mapping $\Theta\mapsto\Theta_{\perp}$ where
$\theta\in\Theta$ denotes the parameters of the reference model and
$\theta_{\perp}\in\Theta_{\perp}$ the parameters of the projected submodel
$M_{\perp}$. Given a set of posterior draws $\\{\theta_{s}\\}_{s=1}^{S}$ for
the reference model, we can then project these individually according to
Equation (3.4) to obtain the projected posterior draws
$\\{\theta_{\perp}^{s}\\}_{s=1}^{S}$ for any submodel $M_{\perp}$. It can be
shown that the projected regression coefficients are obtained analytically by
the least squares solution where the target values are replaced by the fit of
the reference model for a particular posterior draw (see the appendix of
Piironen and Vehtari, 2017a). The projection loss is then defined as the mean
discrepancy over the posterior of the reference model

$\displaystyle
L=\frac{1}{S}\sum_{s=1}^{S}\delta(\theta_{s},\theta_{\perp}^{s}).$ (3.5)

We then seek a reduced model with minimal loss by using a forward search
strategy, in which nodes are sequentially added that minimize $\delta$. A
direct consequence of defining the model parameters of the submodels as
projections of the reference model is that the prior is only specified for the
reference model and this information is also then transmitted to the submodels
in the projection.

To assess the accuracy of the reduced models in the second step, we use
approximate leave-one-out (LOO) cross-validation with Pareto smoothed
importance sampling (PSIS), which avoids the repeated fitting of the reference
model, but requires the selection to be performed $n$ times. In (PS)IS-LOO,
the posterior draws can be treated as draws from the LOO posteriors after
weighting them by the importance weights

$w_{is}\propto\frac{1}{p(y_{i}\,|\,\theta_{s})}$ (3.6)

which are then smoothed to stabilize the LOO estimates (Vehtari et al., 2017).
Thus the projection is carried out with the data point $i$ left out, where the
submodel draws $\\{\theta_{s}\\}_{s=1}^{S}$ are weighted with $w_{is}$ when
predicting the left-out point. This gives an estimate of the predictive
accuracy $\hat{u}_{k}$ with some standard error $s_{k}$ for a given model
complexity $k$ (number of nodes, for instance). As an accuracy measure, we use
log predictive density (LPD).

Our decision rule for achieving sparsity ($\hat{\beta}_{i}^{(j)}=0$) is based
on $\hat{u}_{k}$. We select the simplest model that has an acceptable
difference $\Delta u>0$ relative to the reference model $\hat{u}_{*}$ with
some level of confidence $\alpha$. That is, we choose the simplest model
satisfying

$\Pr{\left[u_{*}-u_{k}\leq\Delta u\right]}\geq 1-\alpha,$ (3.7)

which is estimated using the Bayesian bootstrap. The choices for $\Delta u$
and $\alpha$ are made on subjective grounds depending on how much one is
willing to sacrifice predictive accuracy for making the model simpler. In
other words, sparsifying $\mathcal{G}$ can be viewed as a trade-off between
node-wise predictive loss and parsimony. A reasonable but generic choice for
$\Delta u$ would be 5 or 10 percent of the difference
$\hat{u}_{*}-\hat{u}_{0}$, where $\hat{u}_{0}$ denotes the accuracy estimate
for the simplest possible (null) model (Piironen and Vehtari, 2017a).

###  3.3 Symmetrization

For each node-wise regression, after the selection is completed, the projected
estimates (including zero and non-zero elements) are placed into a $p$
$\times$ $p$ matrix $\Lambda_{i,j\neq i}$. In low-dimensional settings this
matrix would be symmetric, with all estimates having the same sign, but in
high-dimensional settings this is not guaranteed. We use the “or-rule”,
$\left(\hat{\beta}_{j}^{(i)}\vee\hat{\beta}_{i}^{(j)}\right)\neq 0$, to
achieve symmetric non-zero entries in $\Lambda$, in which the neighborhoods
are re-projected to include the inconsistent estimate(s).

####  3.3.1 Partial Correlation Matrix

We then use the approach, described in Krämer et al. (2009), to compute the
partial correlation matrix

$\displaystyle\tilde{\textbf{{P}}}_{i,-i}$
$\displaystyle=[\hat{\rho}_{ij}],\hskip 2.84544ptj\in\\{1,...,p-1\\},$ (3.8)
$\displaystyle\hat{\rho}_{ij}$
$\displaystyle=\text{sign}\big{(}\Lambda\big{)}\text{min}\Big{\\{}1,\sqrt[]{\Lambda_{ij}\hskip
1.42271pt\circ\hskip 1.42271pt\Lambda_{ij}^{\top}}\Big{\\}},$ $\displaystyle
if\hskip 2.84544pt\text{sign}(\Lambda_{ij})=\text{sign}(\Lambda_{ji}),$
$\displaystyle\text{and 0 otherwise},$

where $\circ$ denotes the Hadamard product between the corresponding matrix
elements. This approach ensures that the partial correlations are bounded [-1,
1], and are well defined in high-dimensional settings (Krämer et al., 2009).

####  3.3.2 Precision Matrix

Although the proposed method, with Equation 3.8, provides an estimate for the
partial correlations, it is also important to consider an estimate for
$\Omega$. In low-dimensional settings, following Equation 2.7, this is
straight forward. We thus follow Yuan (2010) and Liu and Wang (2017), each of
which used node-wise regressions to construct $\Omega$ in high-dimensional
settings.

The diagonal and off-diagonal elements, from Equation 2.7, are defined as

$\omega_{ii}=\frac{1}{\text{Var}(e_{k})_{ii}}\hskip 2.84544pt\text{and}\hskip
2.84544pt\omega_{ij}=\frac{-\beta_{ij}}{\text{Var}(e_{k})_{ii}}$ (3.9)

where $\text{Var}(e_{i})$ is the residual variance. With many variables, due
to over-fitting, $\text{Var}(e_{i})$ can become very small, thereby produce
non-optimal estimates for $\omega_{ii}$ and $\omega_{ij}$. We thus computed
the projected residual variance, that included only the selected nodes for
each neighborhood, as

$\displaystyle\hat{\Omega}_{ii}$
$\displaystyle=\frac{1}{\text{Var}(e_{k})_{ii}},$ (3.10)
$\displaystyle\text{Var}(e_{k})_{ii}$ $\displaystyle=y-\Lambda_{i,-i}X_{-i},$

where $y=X_{i}$ (Equation 2.4), $X_{-i}$ is the model matrix excluding the
$i$th variable, and $\Lambda_{i,-i}$ is the matrix containing the projected
estimates. We then use Equation 3.10, for all $i\neq j$, to compute the off-
diagonal elements

$\hat{\Omega}_{ij}=-\frac{1}{2}\Bigg{(}\frac{\Lambda_{i,i\neq
j}}{\text{Var}(e_{k})_{ii}}+\frac{\Lambda_{j,j\neq
i}}{\text{Var}(e_{k})_{jj}}\Bigg{)},$ (3.11)

where symmetry is obtained by averaging the values. While there are
alternative strategies, such as
$\tilde{\Omega}_{ij}\leftarrow\text{min}\\{\hat{\Omega}_{ij},\hat{\Omega}_{ji}\\}$,
both were shown to provide an adequate estimate for $\Omega$ (Cai et al.,
2011). Whereas this procedure ensures $\Omega$ is asymptotically positive
definite (Fan et al., 2016), a correction for finite sample positive
definiteness may be needed. We iteratively find the nearest positive definite
matrix, as in Higham (2002) and Houduo and Defeng (2006), but with respect to
the estimated graphical structure, thereby ensuring the edge set remains
intact.

Figure 1: Graphical structures used in the low-dimensional simulations
(section 4.1.2). The precision matrices were standardized, and sign reversed,
resulting in the partial correlation matrix. In the random graph, there was 30
% probability of sharing a connection.

##  4 Numerical Experiments

In the following, we present a series of experiments in which our aim is to
demonstrate the generality of projection predictive covariance selection. We
first examine model selection consistency in low-dimensional settings, and for
high-dimensional settings we also examine various loss functions for $\Omega$.
For both dimensions of data, we include the same estimation methods that were
specifically chosen to provide an informative contrast to the projection
method. The first Bayesian approach, that uses a birth-death MCMC sampler,
allowed for obtaining the BMA and MAP graphical structures (Mohammadi and Wit,
2015b, R package: BDgraph). The second Bayesian approach, the graphical lasso
(BGL), allowed for comparing the projection method to using credible intervals
as the decision rule (Wang, 2012; Khondker et al., 2013, low-dimensional: 95
%; high-dimensional: 50 %). For the classical methods, we included the
customary graphical lasso with EBIC (Friedman et al., 2008; Foygel and Drton,
2010, R package: glasso), in addition to the TIGER method that estimates
$\Omega$ with $\ell_{1}$-regularized regression (Liu and Wang, 2017, R
package: flare). For the former we sequenced through 50 potential values for
lambda (0.01 to 0.50) and assumed 0.5 for the gamma parameter of EBIC, whereas
for the TIGER method lambda was computed as $\sqrt{\text{log}(p)/n}$ (Liu and
Wang, 2017).

For the proposed method, outlined in section 3, we chose generic values for
$\tau_{0}$ = 1 (Equation 3.2), $\Delta u=10\%$ (Equation 3.7), and
$\alpha=90\%$ (Equation 3.7). This choice of $\tau_{0}$ follows common
practice, whereas the values for $\Delta u$ and $\alpha$ were shown to be
reasonable choices for variable selection (Piironen and Vehtari, 2017a). The
models were fitted with the R package GGMprojpred that estimates high-
dimensional regression models with horseshoe and uses projpred for the
projection (Paasiniemi et al., 2017). For both the projection method and the
Bayesian glasso 1,000 samples were drawn from the posterior (excluding a 1,000
iteration warm-up). We used the default settings, as implemented in each
package, to fit the remaining models. All computations were done in R version
3.4.2 (R Core Team, 2016).

Figure 2: Simulation results for the the low-dimensional simulations (section
4.1.2).

###  4.1 Low-dimensional

We first examine low-dimensional settings, in the customary asymptotic setting
(fixed $p$ = 20), to evaluate consistency as $n\in\\{100,250,500,1000\\}$
increased. We considered two graphical structures that are displayed in Figure
1. The first followed an AR-2 process, in which $\omega_{ii}=1$,
$\omega_{i,i-1}=\omega_{i-1,i}=0.5$, and $\omega_{i,i-2}=\omega_{i-2,i}=0.5$,
and $\omega=0$ otherwise, whereas the second followed a random structure with
a 30 % probability of two nodes sharing a connection. For the latter, the
precision matrix was generated from a Wishart distribution $\Omega$ $\sim$
$W_{G}$($df$ = 20, $I_{p}$) with 20 degrees of freedom. For each condition,
the data were generated from a multivariate normal distribution, with means of
zero, and standardized covariance matrices.

Performance was assessed with two loss functions, KL-divergence and quadratic
loss, that are defined as

$\displaystyle\text{KL}(\Omega,\hat{\Omega})$
$\displaystyle=\text{tr}(\Omega^{-1}\hat{\Omega})-\text{log}(|\Omega^{-1}\hat{\Omega}|)-p,$
(4.1) $\displaystyle\text{QL}(\Omega,\hat{\Omega})$
$\displaystyle=\text{tr}(\Omega^{-1}\hat{\Omega}-I_{p})^{2}.$

We included two measures to assess the accuracy of the graphical structure
(edge set identification), including Specificity and the Matthew’s correlation
coefficient

$\mathrm{SP}=\frac{\mathrm{TN}}{\mathrm{TN}+\mathrm{FP}}\hskip
5.69046pt\text{and}\hskip
5.69046pt\mathrm{MCC}=\frac{\mathrm{TP}\times\mathrm{TN}-\mathrm{FP}\times\mathrm{FN}}{\sqrt[]{(\mathrm{TP}+\mathrm{FP})(\mathrm{TP}+\mathrm{FN})(\mathrm{TN}+\mathrm{FP})(\mathrm{TN}+\mathrm{FN})}},$
(4.2)

where $\mathrm{TP}$, $\mathrm{FP}$, $\mathrm{TN}$, and $\mathrm{FN}$ are the
number of true positives, false positives, true negatives, and false
negatives. It should be noted that 1 - SP is the false positive rate, and that
the MCC denotes the the correlation between binary variables.

####  4.1.1 Results: Risk Estimates

The results are displayed in Figure 2. While all methods were consistent, in
that the risk decreased with larger sample sizes, there were some notable
differences. The projection method, for both loss functions, provided better
performance than BMA. However, for KL-divergence, the Bayesian glasso often
had superior performance compared to the projection method. On the other hand,
for both graphical structures, the proposed method always had the lowest risk
as measured with quadratic loss. For the losses, while performance was
competitive with the classical methods, it is important to also consider
specificity. For the AR-2 graphical structure, GL performed well with respect
to KL-divergence, but specificity consistently declined with larger sample
sizes. Importantly, the projection method often had superior performance
compared to TIGER, which highlights that our approach for constructing
$\Omega$, described in section 3.3.2, worked well compared to using
$\ell_{1}$-regularized regression to estimate $\Omega$. Further, and
interestingly, these results make it clear that directly estimating $\Omega$
does not necessarily provide a more accurate estimate.

####  4.1.2 Results: Edge Set Identification

For edge set identification, it should be noted that the projection method
generally had the lowest MCC, except when compared to glasso. However, these
results also show that specificity did not decline with larger sample sizes,
and that the false positive rate (1- SP) was consistently lower than the
classical methods. In fact, specificity was mostly uninfluenced by the sample
size, whereas performance varied for the other methods. The alternative
Bayesian methods performed well in these performance measures, and in
particular BMA and the Bayesian glasso. For the latter, specificity was
consistently around 95 % which suggest the intervals were calibrated with
respect to covering zero, thereby suggesting minimal shrinkage. This makes
sense, given the dimensions of these data, in that 1 - SP denotes the $\alpha$
level.

###  4.2 High-dimensional

We now examine performance in high-dimensional settings, where $n$ was fixed
to 50 and $p\in\\{25,50,100,150\\}$ increased. We considered four additional
graphical structures (Figure: 1):

  1. 1)

AR-1: $\Omega^{-1}_{ij}$ = $0.7^{|i-j|}$, and $\Omega^{-1}_{ii}$ = 1.

  2. 2)

Cluster: A graph that contains clusters of connections, each of which are
randomly structured graphs. The number maximum number of clusters is
max{2,[$p$/20]}, with $\Omega\sim W_{g}(3,I_{p})$.

  3. 3)

Random: 90 % of the off-diagonal elements set to zero, and $\Omega\sim
W_{g}(3,I_{p})$

  4. 4)

Scale-free: A graph generated with the B-A algorithm (Albert and Barabasi,
2001), with $p$ \- 1 edges, and $\Omega\sim W_{g}(3,I_{p})$.

For each condition, 50 datasets were generated from a multivariate normal
distribution, with means of zero, and standardized covariance matrices. It has
been noted that methods may perform well in some respects, but not others. We
thus considered two more loss functions, the first of which was the
$L_{2}$-loss defined as
$L_{2}(\Omega,\hat{\Omega})=||\Omega-\hat{\Omega}||_{F}$. The second loss was
mean squared error in which the risk was assessed for the partial correlation
matrix. Additionally, we computed sensitivity and the F1-score

$\mathrm{SN}=\frac{\mathrm{TP}}{\mathrm{TP}+\mathrm{FN}}\hskip
4.26773pt\text{and}\hskip
4.26773ptF1-\text{score}=\frac{\text{2}\mathrm{TP}}{\text{2}\mathrm{TP}+\mathrm{FP}+\mathrm{FN}}.$
(4.3)

The F1-score can range between 0 and 1, with 1 denoting perfect
identification.

Table 1: Speed of estimation (n = 50) averaged across graphical structures. The parentheses include the standard deviations. Only the results for BMA are displayed (the MAP times are identical). BGL: Bayesian glasso; GL: glasso; TIGER: tuning insensitive graph estimation and regression; BMA: Bayesian model averaging; Projection: projection predictive selection. $p$ | BGL | GL | TIGER | BMA | Projection |
---|---|---|---|---|---|---
25 | 1.72 (0.10) | 0.04 (0.01) | 0.19 (0.02) | 3.45 (0.22) | 13.70 (0.24) |
50 | 7.84 (0.07) | 0.55 (0.28) | 0.29 (0.05) | 14.40 (0.64) | 23.30 (0.55) |
100 | 52.60 (0.46) | 4.98 (2.81) | 0.72 (0.20) | 104.00 (5.85) | 40.00 (1.22) |
150 | 129.00 (0.96) | 23.40 (10.90) | 1.50 (0.32) | 345.00 (7.58) | 60.90 (1.13) |

####  4.2.1 Results: Speed of Estimation

We first considered the speed of estimation, with the times averaged across
graphical structures. These comparisons were non-trivial. To ensure a fair
comparison, we implemented each method in the original source code language,
but with R serving as a front end. For example, the Bayesian glasso was fitted
in Matlab, but was called from the Julia language which was in turn called
from R. Not surprisingly the classical methods were much faster, and in
particular, the TIGER method that does not perform model selection. On the
other hand, glasso had notable variation in the timings which was likely a
result of increased estimation time for lower values of the tuning parameter.

Among the Bayesian methods, the projection method scaled with $p$, which is
important for high-dimensional applications. While the projection method was
slowest in the lowest dimension ($p=25$: 13.70 seconds), the time for the
highest-dimension setting was the fastest ($p=150$: 60.90). Indeed, for
$p=150$, the projection method was more than 5 times faster than Bayesian
model averaging with a birth-death MCMC algorithm.

Figure 3: Example graphical structures ($p$ = 50) used in the high-dimensional
simulations (section 1). The precision matrices were standardized, and sign
reversed, resulting in the partial correlation matrix. In the random graph,
for each covariance, there was a 10 % probability of a conditional
relationship.

####  4.2.2 Results: Risk Estimates

These simulation results are presented in Figure 4. We only displayed the BMA
results, since they were very similar to the MAP estimates. Most noticeably,
the performance of BMA was consistently poor in the high-dimensional setting.
In almost all conditions, the risk was highest for the BMA precision matrices
($\Omega$), in addition to the partial correlation matrices (P). Specifically,
for P ($p$ = 150), the MSE was substantially higher than the projection
method. For the AR-1 structure, that is, the risk was 4 times greater than for
the proposed Bayesian method. Compared to the Bayesian glasso, the projection
method also had superior performance in most conditions, particularly in the
high-dimensional settings. Of course, it is possible that the decision rule
(50 % credible) could be adjusted to improve performance. We discuss this
further in the application section, in addition to the discussion.

The projection method not only provided superior performance compared to the
Bayesian methods, but also compared to the classical methods. The TIGER
method, which also constructs $\Omega$ from regression, provides an excellent
comparison to our method. The risk estimates were strikingly similar for P. In
fact, these results make it clear that direct estimation of $\Omega$ is not
necessarily more accurate. There were some notable difference compared to
TIGER, in that the projection method excelled when the measure was either KL-
divergence or the quadratic loss. Further, for the scale-free structure, we
observed that the TIGER method had increased risk for $\Omega$, in addition to
substantial variability. On the other hand, across all conditions, the
projection method either had the lowest risk or was competitive, thereby
demonstrating the accuracy of our method for constructing $\Omega$ (section:
2.7).

While glasso is one of the most popular methods, it is also the case that more
recent methods regularly show superior performance. We thus included glasso as
a point of reference, and note the projection method generally had lower risk.
In particular, the risk estimates for P was often much higher for the glasso
method. The most competitive performance was for quadratic loss, although the
projection method typically had superior performance.

Figure 4: Risk estimates for the high-dimensional simulations (section: 1).
GL: glasso; TIGER: tuning insensitive graph estimation and regression; BMA:
Bayesian model averaging; MAP: maximum a posteriori; Projection: projection
predictive selection.

####  4.2.3 Results: Edge Set Identification

The simulation results are presented in Figure 5. For F1-scores and the MCC,
the projection method often had the highest scores. In particular, for the
AR-1 and scale-free graphical structures. Interestingly, the performances
scores were very similar compared to TIGER which also uses a regression based
approach. As such, like the risk estimates, the results indicated that direct
estimation of $\Omega$ is not necessarily more accurate.

Compared to BMA and MAP estimates, for the high-dimensional settings, the
projection method always had superior performance. Notably, the BMA and MAP
graphs consistently had F1-scores and MCC’s of basically zero ($p$ = 150).
This is important, because novel estimation methods are necessary for high-
dimensional settings, and the BMA and MAP based methods seem most suited for
low-dimensional settings. Further, specificity was low for both BMA and MAP
(thus the false positive rate was inflated compared to the projection method),
which is especially important for high-dimensional graphs in which there are
potentially thousands of covariances to consider. While the results for the
Bayesian glasso depend on the chosen credible interval (50 %), the general
pattern of results would apply to all decision rules. That is, while
specificity approached 90 % for the high-dimensional settings, it is clear the
decision rule would need to change for each dimension of data. That is,
specificity consistently diminished with lower dimensions. On the other hand,
across most conditions, the projection method had the most consistent
performance for specificity, thereby suggesting the generality of our method.
However, it should be noted that sensitivity to detect true edges was also
lower than the other methods.

These results also make it clear that the projection method compares favorably
to glasso. Most notably, the glasso method appears to be erratic (e.g., the
scale-free structure) and the most variable. This may be due to performing
model selection with EBIC, without ensuring that the assumptions of glasso
were satisfied. Within the classical methods, it appears that the TIGER method
is a more viable approach for graph estimation.

Figure 5: Edge identification for the high-dimensional simulations (section
1). GL: glasso; TIGER: tuning insensitive graph estimation and regression;
BMA: Bayesian model averaging; MAP: maximum a posteriori; Projection:
projection predictive selection. Figure 6: Sachs protein-signaling network.
BGL: Bayesian glasso; GL: glasso; TIGER: tuning insensitive graph estimation
and regression; BMA: Bayesian model averaging; MAP: maximum a posteriori;
Projection: projection predictive selection.

##  5 Applications

###  5.1 Low-dimensional

We now compare methods on the so-called Sach’s network ($p=11$ and $n=7,466$).
These data were first described in Sachs et al. (2005), where they used
Bayesian methods to suggest causal protein-signaling between phosphorylated
proteins and phospholipids. The identified edges were then confirmed with
experimental manipulations, thereby providing a unique opportunity to evaluate
a known graphical structure (i.e., a “ground truth”). Because this network is
known, it has been used extensively to evaluate the performance of novel GGM
methods (Friedman et al., 2008; Khondker et al., 2013) or decision rules for
existing methods (Kuismin and Sillanpää, 2017). Based on the present
simulations, we used 95 % credible intervals for the Bayesian glasso.

The estimated graphical structures are displayed in Figure 6, and the
performance scores are reported in Table 2. There were clear differences
between projection selection and the other methods. Most notably, the
alternative methods all estimated very dense graphs in which specificity did
not exceed 0.31 (FRP ¿ 0.69). On the other hand, the projection method had the
lowest false positive rate but also had the lowest sensitivity. From the
estimated graphs, it is clear that the projection method appears to only
include edges with strong signals (based on the thickness of the connections).
While the F1-scores and the MCC were lowest for the projection method, it was
also the case that it was the only method that estimated more true positives
than negatives.

Of course, each method has flexibility in model specification, which can
possibly improve estimation. We explored this possibility with all three
methods. We first examined several $\gamma$ parameters for EBIC model
selection with glasso. With a large value ($\gamma=50$), we were able to
decrease the false positive rate from 0.78 to 0.53. These results may reflect
that the $\ell_{1}$-regularized estimates approach maximum likelihood when the
sample size becomes large (Kuismin and Sillanpää, 2017). For the BMA and MAP
networks, we made several adjustments to the default settings including (1)
more shrinkage, to an identity matrix, by increasing the degrees of freedom
parameter $r$ of the Wishart prior distribution (default = 3); and (2) a
larger cut-off for the posterior inclusion probability (PIP; default = 0.5).
With $r=100$, the false positive rate decreased by only 0.05 which suggests
diminishing returns from the prior distribution in large sample settings. We
then set the PIP cut-off to 0.99, and this further reduced the false positive
rate from 0.72 to 0.58. We also increased the number of iterations, allowing
for more samples from the most probable graphs, but this did not have much of
an effect on the false positive rate.

Table 2: Performance measures for Sach’s protein-signaling network. BGL: Bayesian glasso; GL: glasso; TIGER: tuning insensitive graph estimation and regression; BMA: Bayesian model averaging; MAP: maximum a posteriori; Projection: projection predictive selection. FPR: False positive rate. | BGL | GL | TIGER | BMA | MAP | Projection
---|---|---|---|---|---|---
Sensitivity | 0.95 | 1.00 | 0.95 | 0.95 | 0.95 | 0.37
Specificity | 0.22 | 0.22 | 0.25 | 0.31 | 0.28 | 0.83
F1-score | 0.55 | 0.58 | 0.56 | 0.58 | 0.57 | 0.44
MCC | 0.22 | 0.30 | 0.24 | 0.29 | 0.27 | 0.23
FPR | 0.78 | 0.78 | 0.75 | 0.69 | 0.72 | 0.17
Figure 7: CEU genetic network. BGL: Bayesian glasso; GL: glasso; TIGER: tuning
insensitive graph estimation and regression; BMA: Bayesian model averaging;
MAP: Maximum a posteriori; Projection: projection predictive selection.

###  5.2 High-dimensional

We now analyze a high-dimensional dataset, which includes 100 genes ($p$) from
60 ($n$) unrelated individuals of Northern and Western European ancestry from
Utah (CEU). This data is freely available (GGMprojpred), and thoroughly
described in (Bhadra and Mallick, 2013). Deciding upon the appropriate
credible interval, for the Bayesian glasso, was non-trivial. We ultimately
chose a 50 % interval, which was also used in Khondker et al. (2013) and Li et
al. (2017), but note this can be adjusted to reduce or increase the number of
edges.

The estimated graphical structures are displayed in Figure 7. The results
correspond closely to the simulation results. That is, BMA (843 edges), MAP
(982 edges), and the Bayesian glasso estimated very dense graphs (1,020
edges). We were able to reduce the number of connections, by adjusting the
default settings, but the resulting graphs for BMA and MAP remained much
denser than the projection method (145 edges). In fact, the proposed Bayesian
method provided the sparsest solution, although the glasso (157 edges) and
TIGER (242 edges) methods were similar. In high-dimensional settings, with
potentially thousands of partial correlations, this provides a more
interpretable graph.

##  6 Conclusion

We introduced a Bayesian method for estimating Gaussian graphical models with
projection predictive selection (section 3). For determining conditional
relationships, the proposed method is framed in terms of tolerable predictive
loss from a reference model. This decision theoretic framework stands in
contrast to posterior model probabilities or credible interval exclusion of
zero. In contrast to customary neighborhood selection approaches, that
estimate the partial correlation or adjacency matrix, our method provides an
estimate of the inverse-covariance matrix. With numerical experiments, we
demonstrated that projection predictive covariance selection often outperforms
classical and Bayesian methods, both in terms of accurate estimation and edge
set identification. For specificity (and thus false positives) in particular,
projection selection consistently had the best performance and the scores were
largely independent of the simulation condition. Further, in the application
section, we showed that our approach similarly produced the sparsest graphical
structures, thereby providing a more interpretable network in high-dimensional
settings.

In the classical literature, substantial effort has gone into developing
regression based approaches for covariance selection (Meinshausen and
Bühlmann, 2006; Krämer et al., 2009; Hornik, 2009). These methods often
outperform direct estimation, and have less restrictive assumptions (Buhlmann
and van de Geer, 2011). In high-dimensional settings, the projection method
had consistently lower risk estimates than BMA and MAP. In particular, the
projection method excelled with respect to the log-score (KL-divergence)
Hoeting et al. (1999). In other words, while BMA is known to minimize the log-
score, these simulations showed that projection selection is more accurate in
high-dimensional settings. Direct estimation was also not faster than fitting
several regression models, which further highlights an additional advantage of
the proposed method. Further, while projective covariance selection finds a
trade-off between sparsity and accuracy, it should be noted that sensitivity
was consistently lower than the alternative methods. This was not unexpected,
for instance, if there are correlating features that carry similar
information, then the projection tends to select only one (or a few) of these
and ignore the rest because this does not introduce considerable loss in
predictive accuracy. For this reason, the projection method might not be ideal
when the goal is to find all the true edges, but rather when we want the
simplest possible graph that has almost as low KL/QL compared to a denser
graphical structure.

Despite these promising findings, none of the methods provided the best
performance for all conditions or performance measures. In section 5, for
example, we actively pursued a sparse graph but discovery (minimizing false
negatives) could be important to consider. We focused on specificity for two
reasons: (1) based on the results, we attempted to outperform projection
selection by adjusting the settings of the alternative methods under
consideration. This strengthened our conclusion that projection selection
estimates the fewest false positives; and (2) the classical methods produced
such dense graphs that it was not reasonable to increase the density. Of
course, while not explored here, the projection method readily allows for
denser graphs by adjusting the decision rule for setting coefficients to zero
(section 3).

There are several important limitations for the our method based on projection
predictive selection. First, our method is best suited for moderately high-
dimensional settings which potentially limits its use in ultra high-
dimensions. Second, our method does not initially guarantee a positive
definite estimate. However, even with a correction, we showed that the
projected precision matrix was often more accurate than direct estimation and
a classical method that similarly estimates $\Omega$ with multiple regression.
Third, our conclusions are restricted to these methods, conditions, graphical
structures, and decision rules. For example, in addition to EBIC based tuning
parameter selection, glasso can be used with the rotation information
criterion, stability approach for regularization selection, and cross-
validation. Finally, we did not compare our method to several approximate
Bayesian approaches (Leppä-aho et al., 2017). Our focus was on fully Bayesian
methods, whereas approximate methods often provide only an adjacency matrix,
thereby limiting the losses that can be evaluated.

The proposed method can be extended in a number of important ways; for
example, to handle mixed graphical models (continuous and binary variables),
longitudinal or time series data, or those data in which normative assumptions
are not warranted. Additionally, the projection can be carried out on the
inverse-covariance matrix directly and alternative samplers can be
implemented, which would improve efficiency and allow for improved scalability
to higher dimensions. These future directions will be addressed in our future
work on this topic.

## Acknowledgements

Research reported in this publication was supported by three funding sources:
(1) The National Academies of Sciences, Engineering, and Medicine FORD
foundation pre-doctoral fellowship to DRW; (2) The National Science Foundation
Graduate Research Fellowship to DRW; and (3) the National Institute On Aging
of the National Institutes of Health under Award Number R01AG050720 to PR. The
content is solely the responsibility of the authors and does not necessarily
represent the official views of the National Academies of Sciences,
Engineering, and Medicine, the National Science Foundation, or the National
Institutes of Health.

## References

  * Albert and Barabasi (2001) Albert, R. and Barabasi, A. L. (2001).  “Statistical mechanics of complex networks.”
  * Anderson (2003) Anderson, T. W. (2003).  An Introduction to Multivariate Statistical Analysis.  London: Wiley-Interscience.
  * Atay-Kayis and Massam (2005) Atay-Kayis, a. and Massam, H. (2005).  “A Monte Carlo method to compute the marginal likelihood in non decomposable graphical Gaussian models.”  Biometrika, 92(2): 317–335.
URL http://www.math.yorku.ca/Who/Faculty/Massam/submitted/km15pt.pdf

  * Baba et al. (2004) Baba, K., Shibata, R., and Sibuya, M. (2004).  “Partial Correlation and Conditional Correlation as Measures of Conditional Independence.”  Australian & New Zealand Journal of Statistics, 46(4): 657–664.
URL http://doi.wiley.com/10.1111/j.1467-842X.2004.00360.x

  * Baba and Sibuya (2005) Baba, K. and Sibuya, M. (2005).  “Equivalence of Partial and Conditional Correlation Coefficients.”  J. Japan Statist. Soc., 35(1): 1–19.
  * Bhadra and Mallick (2013) Bhadra, A. and Mallick, B. K. (2013).  “Joint high-dimensional Bayesian variable and covariance selection with an application to eQTL analysis.”  Biometrics, 69(2): 447–457.
  * Buhlmann and van de Geer (2011) Buhlmann, P. and van de Geer, S. (2011).  Statistics for high-dimensional data: methods, theory and applications.  Berlin: Springer Science & Business Media.
  * Cai et al. (2011) Cai, T., Liu, W., and Luo, X. (2011).  “A constrained ℓ1minimization approach to sparse precision matrix estimation.”  Journal of the American Statistical Association, 106(494): 594–607.
  * Das et al. (2017) Das, A., Sampson, A., Lainscsek, C., Muller, L., Lin, W., Doyle, J., Cash, S., Halgren, E., and Sejnowski, T. (2017).  “Interpretation of the Precision Matrix and Its Application in Estimating Sparse Brain Connectivity during Sleep Spindles from Human Electrocorticography Recordings.”  Neural computation, 29(3): 603–642.
  * Dempster (1972) Dempster, A. (1972).  “Covariance Selection.”  Biometrics, 28(1): 157–175.
URL http://www.jstor.org/stable/2528966?origin=crossref

  * Dobra et al. (2010) Dobra, A., Lenkoski, A., and Rodriguez, A. (2010).  “Bayesian inference for general Gaussian graphical models with application to multivariate lattice data.”  1–30.
URL http://arxiv.org/abs/1005.4094

  * Drton and Perlman (2004) Drton, M. and Perlman, M. D. (2004).  “Model selection for Gaussian concentration graphs.”  Biometrika, 91(3): 591–602.
  * Dupuis and Robert (2003) Dupuis, J. A. and Robert, C. P. (2003).  “Variable selection in qualitative models via an entropic explanatory power.”  Journal of Statistical Planning and Inference, 111(1-2): 77–94.
  * Fan et al. (2016) Fan, J., Liao, Y., and Liu, H. (2016).  “An overview of the estimation of large covariance and precision matrices.”  Econometrics Journal, 19(1): C1–C32.
  * Foygel and Drton (2010) Foygel, R. and Drton, M. (2010).  “Extended Bayesian Information Criteria for Gaussian Graphical Models.”  Advances in Neural Information Processing Systems, 604–612.
  * Friedman et al. (2008) Friedman, J., Hastie, T., and Tibshirani, R. (2008).  “Sparse inverse covariance estimation with the graphical lasso.”  Biostatistics, 9(3): 432–441.
URL https://academic.oup.com/biostatistics/article-
lookup/doi/10.1093/biostatistics/kxm045

  * George and McCulloch (1993) George, E. I. and McCulloch, R. E. (1993).  “Variable selection via Gibbs sampling.”  Journal of the American Statistical Association, 88(August 2015): 881–889.
  * Goutis (1998) Goutis, C. (1998).  “Model choice in generalised linear models: a Bayesian approach via Kullback-Leibler projections.”  Biometrika, 85(1): 29–37.
URL http://biomet.oupjournals.org/cgi/doi/10.1093/biomet/85.1.29

  * Green (1995) Green, P. J. (1995).  “Reversible Jump Markov Chain Monte Carlo Computation and Bayesian Model Determination.”  Biometrika, 82(4): 711–732.
URL http://www.jstor.org/stable/2337340?origin=crossref

  * Hartlap et al. (2007) Hartlap, J., Simon, P., and Schneider, P. (2007).  “Unbiased estimation of the inverse covariance matrix.”  a&a, 464: 399–404.
URL https://www.aanda.org/articles/aa/pdf/2007/10/aa6170-06.pdf

  * Higham (2002) Higham, N. J. (2002).  “Computing the nearest correlation matrix - A problem from finance.”  IMA Journal of Numerical Analysis, 22(3): 329–343.
  * Hoeting et al. (1999) Hoeting, J. A., Madigan, D., Raftery, A. E., and Volinsky, C. T. (1999).  “Bayesian Model Averaging: A Tutorial.”  Statistical Science, 14(4): 382–417.
  * Hornik (2009) Hornik, K. (2009).  “Journal of Statistical Software.”  Journal Of Statistical Software, 31(2): 2008–2010.
  * Houduo and Defeng (2006) Houduo, Q. and Defeng, S. (2006).  “A QUADRATICALLY CONVERGENT NEWTON METHOD FOR COMPUTING THE NEAREST CORRELATION MATRIX.”  SIAM J. MATRIX ANAL. APPL, 28(2): 360–385.
  * Huang and Wand (2013) Huang, A. and Wand, M. P. (2013).  “Simple marginally noninformative prior distributions for covariance matrices.”  Bayesian Analysis, 8(2): 439–452.
URL http://projecteuclid.org/euclid.ba/1369407559

  * Janková and van de Geer (2015) Janková, J. and van de Geer, S. (2015).  “Confidence intervals for high-dimensional inverse covariance estimation.”  Electronic Journal of Statistics, 9: 1205–1229.
  * Khondker et al. (2013) Khondker, Z. S., Zhu, H., Chu, H., Lin, W., and Ibrahim, J. G. (2013).  “The Bayesian Covariance Lasso.”  Statistics and its interface, 6(2): 243–259.
URL http://www.ncbi.nlm.nih.gov/pubmed/24551316

  * Krämer et al. (2009) Krämer, N., Schäfer, J., and Boulesteix, A.-L. (2009).  “Regularized estimation of large-scale gene association networks using graphical Gaussian models.”  BMC Bioinformatics, 10(1): 384.
URL http://www.biomedcentral.com/1471-2105/10/384

  * Kubokawa and Srivastava (2008) Kubokawa, T. and Srivastava, M. S. (2008).  “Estimation of the precision matrix of a singular Wishart distribution and its application in high-dimensional data.”  Journal of Multivariate Analysis, 99(9): 1906–1928.
  * Kuismin and Sillanpää (2017) Kuismin, M. and Sillanpää, M. (2017).  “Estimation of covariance and precision matrix, network structure, and a view toward systems biology.”  Wiley Interdisciplinary Reviews: Computational Statistics, 9(6): 1–13.
  * Kuismin and Sillanpää (2016) Kuismin, M. and Sillanpää, M. J. (2016).  “Use of Wishart prior and simple extensions for sparse precision matrix estimation.”  PLoS ONE, 11(2): e0148171.
URL http://www.pubmedcentral.nih.gov/articlerender.fcgi?artid=PMC4734711

  * Ledoit and Wolf (2004a) Ledoit, O. and Wolf, M. (2004a).  “A well-conditioned estimator for large-dimensional covariance matrices.”  Journal of Multivariate Analysis, 88(2): 365–411.
  * Ledoit and Wolf (2004b) — (2004b).  “Honey, I Shrunk the Sample Covariance Matrix.”  The Journal of Portfolio Management, 30(4): 110–119.
URL http://www.iijournals.com/doi/abs/10.3905/jpm.2004.110

  * Lenkoski (2013) Lenkoski, A. (2013).  “A direct sampler for G-Wishart variates.”  Stat, 2(1): 119–128.
  * Leppä-aho et al. (2017) Leppä-aho, J., Pensar, J., Roos, T., and Corander, J. (2017).  “Learning Gaussian graphical models with fractional marginal pseudo-likelihood.”  International Journal of Approximate Reasoning, 83: 21–42.
URL http://dx.doi.org/10.1016/j.ijar.2017.01.001

  * Li and Zhang (2017) Li, F.-q. and Zhang, X.-s. (2017).  “Bayesian Lasso with neighborhood regression method for Gaussian graphical model.”  Acta Mathematicae Applicatae Sinica, English Series, 33(2): 485–496.
URL http://link.springer.com/10.1007/s10255-017-0676-z

  * Li et al. (2017) Li, Y., Craig, B. A., and Bhadra, A. (2017).  “The Graphical Horseshoe Estimator for Inverse Covariance Matrices.”
URL http://arxiv.org/abs/1707.06661

  * Liu and Wang (2017) Liu, H. and Wang, L. (2017).  “TIGER: A tuning-insensitive approach for optimally estimating gaussian graphical models.”  Electronic Journal of Statistics, 11(1): 241–294.
  * McNally et al. (2015) McNally, R. J., Robinaugh, D. J., Wu, G. W. Y., Wang, L., Deserno, M. K., and Borsboom, D. (2015).  “Mental Disorders as Causal Systems.”  Clinical Psychological Science, 3(6): 836–849.
URL http://journals.sagepub.com/doi/10.1177/2167702614553230

  * Meinshausen and Bühlmann (2006) Meinshausen, N. and Bühlmann, P. (2006).  “High-dimensional graphs and variable selection with the Lasso.”  Annals of Statistics, 34(3): 1436–1462.
  * Mohammadi and Wit (2015a) Mohammadi, A. and Wit, E. C. (2015a).  “Bayesian structure learning in sparse Gaussian graphical models.”  Bayesian Analysis, 10(1): 109–138.
  * Mohammadi and Wit (2015b) — (2015b).  “BDgraph: An R Package for Bayesian Structure Learning in Graphical Models.”
URL http://arxiv.org/abs/1501.05108

  * O’Hara and Sillanpää (2009) O’Hara, R. B. and Sillanpää, M. J. (2009).  “A review of bayesian variable selection methods: What, how and which.”  Bayesian Analysis, 4(1): 85–118.
  * Paasiniemi et al. (2017) Paasiniemi, M., Piironen, J., Gabry, J., and Vehtari, A. (2017).  “projpred (R package version 0.6.0).”
URL https://github.com/stan-dev/projpred

  * Peng et al. (2009) Peng, J., Wang, P., Zhou, N., and Zhu, J. (2009).  “Partial Correlation Estimation by Joint Sparse Regression Models.”  Journal of the American Statistical Association I, 104(486): 735–746.
URL http://arxiv.org/abs/0811.4463

  * Piironen and Vehtari (2017a) Piironen, J. and Vehtari, A. (2017a).  “Comparison of Bayesian predictive methods for model selection.”  Statistics and Computing, 27(3): 711–735.
  * Piironen and Vehtari (2017b) — (2017b).  “Sparsity information and regularization in the horseshoe and other shrinkage priors.”
URL http://arxiv.org/abs/1707.01694

  * R Core Team (2016) R Core Team (2016).  “R: A Language and Environment for Statistical Computing.”
URL https://www.r-project.org/

  * Roverato and Castelo (2017) Roverato, A. and Castelo, R. (2017).  “The networked partial correlation and its application to the analysis of genetic interactions.”  Journal of the Royal Statistical Society. Series C: Applied Statistics, 66(3): 647–665.
  * Sachs et al. (2005) Sachs, K., Perez, O., Pe’er, D., Lauffenburger, D. A., and Nolan, G. P. (2005).  “Causal protein-signaling networks derived from multiparameter single-cell data.”  Science, 308(5721): 523–529.
  * Schäfer and Strimmer (2005a) Schäfer, J. and Strimmer, K. (2005a).  “A Shrinkage Approach to Large-Scale Covariance Matrix Estimation and Implications for Functional Genomics.”  Statistical Applications in Genetics and Molecular Biology, 4(1).
URL
https://www.degruyter.com/view/j/sagmb.2005.4.issue-1/sagmb.2005.4.1.1175/sagmb.2005.4.1.1175.xml

  * Schäfer and Strimmer (2005b) — (2005b).  “An empirical Bayes approach to inferring large-scale gene association networks.”  Bioinformatics, 21(6): 754–764.
URL https://www.aanda.org/articles/aa/pdf/2007/10/aa6170-06.pdf

  * Stephens (1998) Stephens, G. (1998).  “On the Inverse of the Covariance Matrix in Portfolio Analysis.”  The Journal of Finance, 53(5): 1821–1827.
  * Tenenhaus et al. (2008) Tenenhaus, A., Guillemot, V., Gidrol, X., and Frouin, V. (2008).  “Gene Association Networks from Microarray Data using a Regularized Estimation of Partial Correlation based on.”  IEEE Transactions on Computational Biology and Bioinformatics, 7(2): 251–262.
  * van der Pas et al. (2016) van der Pas, S., Szabó, B., and van der Vaart, A. (2016).  “Uncertainty quantification for the horseshoe.”  (4): 1221–1274.
URL http://arxiv.org/abs/1607.01892

  * Vehtari et al. (2017) Vehtari, A., Gelman, A., and Gabry, J. (2017).  “Practical Bayesian model evaluation using leave-one-out cross-validation and WAIC.”  Statistics and Computing, 27(5): 1413–1432.
  * Vehtari and Ojanen (2012) Vehtari, A. and Ojanen, J. (2012).  “A survey of Bayesian predictive methods for model assessment, selection and comparison.”  Statistics Surveys, 6(0): 142–228.
URL http://projecteuclid.org/euclid.ssu/1356628931

  * Wang (2012) Wang, H. (2012).  “Bayesian graphical lasso models and eficient posterior computation.”  Bayesian Analysis, 7(4): 867–886.
URL https://projecteuclid.org/download/pdfview_1/euclid.ba/1354024465

  * Wang (2016) Wang, Z. (2016).  “Transactivation of epidermal growth factor receptor by g protein-coupled receptors: Recent progress, challenges and future research.”
URL
http://www.ncbi.nlm.nih.gov/pubmed/26771606%****␣main.bbl␣Line␣450␣****http://www.pubmedcentral.nih.gov/articlerender.fcgi?artid=PMC4730337

  * Williams and Rast (2018) Williams, D. R. and Rast, P. (2018).  “Back to the basics: Rethinking partial correlation network methodology.”  1–15.
  * Williams et al. (2018) Williams, D. R., Wysocki, A. C., and Rast, P. (2018).  “On Non-Regularized Estimation of Psychological Networks.”
  * Won et al. (2009) Won, J.-h., Lim, J., and Kim, S.-J. (2009).  “Maximum Likelihood Covariance Estimation.”  Asilomar Conference on Signals, Systems, and Computers, 1445–1449.
URL
https://pdfs.semanticscholar.org/f82b/c8eb44f5714a1ffcb8ab25aec401c5ca43c0.pdf

  * Yuan (2010) Yuan, M. (2010).  “High Dimensional Inverse Covariance Matrix Estimation via Linear Programming.”  Journal of Machine Learning Research, 11: 2261–2286.
  * Yuan and Wang (2013) Yuan, T. and Wang, J. (2013).  “A coordinate descent algorithm for sparse positive definite matrix estimation.”  Statistical Analysis and Data Mining, 6(5): n/a–n/a.
URL http://doi.wiley.com/10.1002/sam.11185