  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
