Fitting very large sparse Gaussian graphical models

Gaussian Graphical Models (GGMs) are extensively used in many research areas, such as genomics, proteomics, neuroimaging, and psychology, to study the partial correlation structure of a set of variables. This structure is visualized by drawing an undirected network, in which the variables constitute the nodes and the partial correlations the edges. In many applications, it makes sense to impose sparsity (i.e., some of the partial correlations are forced to zero) as sparsity is theoretically meaningful and/or because it improves the predictive accuracy of the fitted model. However, as we will show by means of extensive simulations, state-of-the-art estimation approaches for imposing sparsity on GGMs, such as the Graphical lasso, ℓ1 regularized nodewise regression, and joint sparse regression, fall short because they often yield too many false positives (i.e., partial correlations that are not properly set to zero). In this paper we present a new estimation approach that allows to control the false positive rate better. Our approach consists of two steps: First, we estimate an undirected network using one of the three state-of-the-art estimation approaches. Second, we try to detect the false positives, by flagging the partial correlations that are smaller in absolute value than a given threshold, which is determined through cross-validation; the flagged correlations are set to zero. Applying this new approach to the same simulated data, shows that it indeed performs better. We also illustrate our approach by using it to estimate (1) a gene regulatory network for breast cancer data, (2) a symptom network of patients with a diagnosis within the nonaffective psychotic spectrum and (3) a symptom network of patients with PTSD.

A sparse precision matrix can be directly translated into a sparse Gaussian graphical model under the assumption that the data follow a joint normal distribution. This neat property makes high-dimensional precision matrix estimation very appealing in many applications. However, in practice we often face nonnormal data, and variable transformation is often used to achieve normality. In this paper we consider the nonparanormal model that assumes that the variables follow a joint normal distribution after a set of unknown monotone transformations. The nonparanormal model is much more flexible than the normal model while retaining the good interpretability of the latter in that each zero entry in the sparse precision matrix of the nonparanormal model corresponds to a pair of conditionally independent variables. In this paper we show that the nonparanormal graphical model can be efficiently estimated by using a rank-based estimation scheme which does not require estimating these unknown transformation functions. In particular, we study the rank-based graphical lasso, the rank-based neighborhood Dantzig selector and the rank-based CLIME. We establish their theoretical properties in the setting where the dimension is nearly exponentially large relative to the sample size. It is shown that the proposed rank-based estimators work as well as their oracle counterparts defined with the oracle data. Furthermore, the theory motivates us to consider the adaptive version of the rank-based neighborhood Dantzig selector and the rank-based CLIME that are shown to enjoy graphical model selection consistency without assuming the irrepresentable condition for the oracle and rank-based graphical lasso. Simulated and real data are used to demonstrate the finite performance of the rank-based estimators.

Graphical models are used to examine conditional independence among random variables. In this chapter we take graphical models for the multivariate complex normal distribution w.r.t. simple undirected graphs into consideration. This is the first published presentation of these models. Graphical models for the multivariate real normal distribution, also called covariance selection models, have already been studied in the literature. The initial work on covariance selection models is done by Dempster (1972) and Wermuth (1976) and a presentation of these models is given in Eriksen (1992). Graphical models for contingency tables are introduced in statistics by Darroch et al. (1980) and further these are well-described in Lauritzen (1989). Graphical association models are treated in general in Whittaker (1990) and Lauritzen (1993). The complex normal graphical models are quite similar to the covariance selection models. We have chosen to develop this chapter without use of exponential families. We study definition of the model, maximum likelihood estimation and hypothesis testing. To verify some of the results in the chapter we use results from mathematical analysis. These can be found in e.g. Rudin (1987). In graphical models one uses the concentration matrix instead of the variance matrix as it is more advantageous. Therefore we define this matrix and derive a relation which is basic for complex normal graphical models. Afterwards we formally define a complex normal graphical model w.r.t. a simple undirected graph. As these models are used to examine conditional independence of selected pairs of variables given the remaining ones we are mainly interested in inference on the concentration matrix. It is possible to base the maximum likelihood estimation of the concentration matrix on a complex random matrix with mean zero.

Identifying context-specific entity networks from aggregated data is an important task, arising often in bioinformatics and neuroimaging applications. Computationally, this task can be formulated as jointly estimating multiple different, but related, sparse undirected graphical models (UGM) from aggregated samples across several contexts. Previous joint-UGM studies have mostly focused on sparse Gaussian graphical models (sGGMs) and can’t identify context-specific edge patterns directly. We, therefore, propose a novel approach, SIMULE (detecting Shared and Individual parts of MULtiple graphs Explicitly) to learn multi-UGM via a constrained ell 1 minimization. SIMULE automatically infers both specific edge patterns that are unique to each context and shared interactions preserved among all the contexts. Through the ell 1 constrained formulation, this problem is cast as multiple independent subtasks of linear programming that can be solved efficiently in parallel. In addition to Gaussian data, SIMULE can also handle multivariate Nonparanormal data that greatly relaxes the normality assumption that many real-world applications do not follow. We provide a novel theoretical proof showing that SIMULE achieves a consistent result at the rate O(log (Kp)/n_{tot}). On multiple synthetic datasets and two biomedical datasets, SIMULE shows significant improvement over state-of-the-art multi-sGGM and single-UGM baselines (SIMULE implementation and the used datasets @https://github.com/QData/SIMULE).

This paper introduces a MATLAB toolbox, MGraph, which applies graphical models as a natural environment to formulate and solve problems in microarray data analysis. MGraph with its graphical interface allows the user to predict genetic regulatory networks by a graphical gaussian model (GGM), and to quantify the effects of different experimental treatment conditions on gene expression profiles by a graphical log-linear model (GLM). The power of graphical models was explored and illustrated through two example applications. First, four MAPK pathways in yeast were meaningfully reconstructed through GGM. Second, GLM was used to quantify the contributions of sex, genotype and age to transcriptional variance in Drosophila melanogaster. This application may provide a valuable aid in the prediction of genetic regulatory networks, as well as in investigations of various experimental conditions that affect global gene expression profiles. The MATLAB program MGraph is freely available at http://www.uio.no/~junbaiw/mgraph/mgraph.html for academics.

Alzheimer's Disease (AD) is the most common neurodegenerative disease in elderly people. Its development has been shown to be closely related to changes in the brain connectivity network and in the brain activation patterns along with structural changes caused by the neurodegenerative process. Methods to infer dependence between brain regions are usually derived from the analysis of covariance between activation levels in the different areas. However, these covariance-based methods are not able to estimate conditional independence between variables to factor out the influence of other regions. Conversely, models based on the inverse covariance, or precision matrix, such as Sparse Gaussian Graphical Models allow revealing conditional independence between regions by estimating the covariance between two variables given the rest as constant. This paper uses Sparse Inverse Covariance Estimation (SICE) methods to learn undirected graphs in order to derive functional and structural connectivity patterns from Fludeoxyglucose (18F-FDG) Position Emission Tomography (PET) data and segmented Magnetic Resonance images (MRI), drawn from the ADNI database, for Control, MCI (Mild Cognitive Impairment Subjects), and AD subjects. Sparse computation fits perfectly here as brain regions usually only interact with a few other areas. The models clearly show different metabolic covariation patters between subject groups, revealing the loss of strong connections in AD and MCI subjects when compared to Controls. Similarly, the variance between GM (Gray Matter) densities of different regions reveals different structural covariation patterns between the different groups. Thus, the different connectivity patterns for controls and AD are used in this paper to select regions of interest in PET and GM images with discriminative power for early AD diagnosis. Finally, functional an structural models are combined to leverage the classification accuracy. The results obtained in this work show the usefulness of the Sparse Gaussian Graphical models to reveal functional and structural connectivity patterns. This information provided by the sparse inverse covariance matrices is not only used in an exploratory way but we also propose a method to use it in a discriminative way. Regression coefficients are used to compute reconstruction errors for the different classes that are then introduced in a SVM for classification. Classification experiments performed using 68 Controls, 70 AD, and 111 MCI images and assessed by cross-validation show the effectiveness of the proposed method.

In chapter one, we proposed a novel multiple phenotype association test methods to test the association between multiple phenotypes and a SNP. Genome-wide association studies (GWAS) have identified many strongly associated genetic variants with phenotypes and have greatly enhanced our understanding of the genetic architecture of complex phenotypes and diseases. In this study, we proposed a new method to construct a Phenotype-Phenotype Network using the sparse Gaussian Graphical Model (sGGM) based on GWAS summary statistics. This approach will isolate the direct relationship between phenotypes, making it easier to identify clusters of phenotypes conditional on other phenotypes that reflect the more meaningful biological or functional connections. Then we applied the community detection method to partition phenotypes into disjoint modules based on the partial correlation matrix of phenotypes and applied multiple phenotype test methods on each of the clusters and combined p-values to get the final P-value. The comprehensive simulation study showed that our method can control the type I error rate effectively and has the highest power compared to the other methods we compared. Application of this method to the GWAS summary statistics of 92 phenotypes from UK-Biobank data has identified higher number of significant SNPs than other methods. Downstream analysis of the significant SNPs shows the functional and biological importance. The complex traits are often the result of the combined effect of multiple genetic variants, each with a small individual effect and these small effects may be missed due to stringent p-value thresholds in GWAS studies. Gene-based association tests combine the effects of multiple variants within a gene and provide more interpretable associations and enable researchers to detect associations that might be missed in single-SNP analyses. Gene-based tests increase statistical power by considering the cumulative effect of several variants in a gene, reducing the multiple testing burden and prioritizing variants based on their likelihood of contributing to the phenotype. The power of the gene-based test depended on the genetic architecture of the trait and the weight of SNPs in a gene. The genetic architecture of the complex trait is not known in advance and the power of the gene-based association tests decreases if the weight is mis-specified. The equal weighting assumption implies that all SNPs within the gene have an equal effect on the trait. On the other hand, minor allele frequency (MAF)-based weights, such as inverse or beta distribution of MAF, upweight rare variants under the assumption that they

The estimation of functional connectivity structure from functional neuroimaging data is an important step toward understanding the mechanisms of various brain diseases and building relevant biomarkers. Yet, such inferences have to deal with the low signal-to-noise ratio and the paucity of the data. With at our disposal a steadily growing olume of publicly available neuroimaging data, it is however possible to improve the estimation procedures involved in connectome mapping. In this work, we propose a novel learning scheme for functional connectivity based on sparse Gaussian graphical models that aims at minimizing the bias induced by the regularization used in the estimation, by carefully separating the estimation of the model support from the coefficients. Moreover, our strategy makes it possible to include new data with a limited computational cost. We illustrate the physiological relevance of he learned prior, that can be identified as a functional connectivity atlas, based on an experiment on 46 subjects of the Human Connectome Dataset.

BackgroundGraphical Gaussian models are popular tools for the estimation of (undirected) gene association networks from microarray data. A key issue when the number of variables greatly exceeds the number of samples is the estimation of the matrix of partial correlations. Since the (Moore-Penrose) inverse of the sample covariance matrix leads to poor estimates in this scenario, standard methods are inappropriate and adequate regularization techniques are needed. Popular approaches include biased estimates of the covariance matrix and high-dimensional regression schemes, such as the Lasso and Partial Least Squares.ResultsIn this article, we investigate a general framework for combining regularized regression methods with the estimation of Graphical Gaussian models. This framework includes various existing methods as well as two new approaches based on ridge regression and adaptive lasso, respectively. These methods are extensively compared both qualitatively and quantitatively within a simulation study and through an application to six diverse real data sets. In addition, all proposed algorithms are implemented in the R package "parcor", available from the R repository CRAN.ConclusionIn our simulation studies, the investigated non-sparse regression methods, i.e. Ridge Regression and Partial Least Squares, exhibit rather conservative behavior when combined with (local) false discovery rate multiple testing in order to decide whether or not an edge is present in the network. For networks with higher densities, the difference in performance of the methods decreases. For sparse networks, we confirm the Lasso's well known tendency towards selecting too many edges, whereas the two-stage adaptive Lasso is an interesting alternative that provides sparser solutions. In our simulations, both sparse and non-sparse methods are able to reconstruct networks with cluster structures. On six real data sets, we also clearly distinguish the results obtained using the non-sparse methods and those obtained using the sparse methods where specification of the regularization parameter automatically means model selection. In five out of six data sets, Partial Least Squares selects very dense networks. Furthermore, for data that violate the assumption of uncorrelated observations (due to replications), the Lasso and the adaptive Lasso yield very complex structures, indicating that they might not be suited under these conditions. The shrinkage approach is more stable than the regression based approaches when using subsampling.

Loss function, unbiasedness, and optimality of Gaussian graphical model selection

Constructing expression networks using transcriptomic data is an effective approach for studying gene regulation. A popular approach for constructing such a network is based on the Gaussian graphical model (GGM), in which an edge between a pair of genes indicates that the expression levels of these two genes are conditionally dependent, given the expression levels of all other genes. However, GGMs are not appropriate for non-Gaussian data, such as those generated in RNA-seq experiments. We propose a novel statistical framework that maximizes a penalized likelihood, in which the observed count data follow a Poisson log-normal distribution. To overcome the computational challenges, we use Laplace's method to approximate the likelihood and its gradients, and apply the alternating directions method of multipliers to find the penalized maximum likelihood estimates. The proposed method is evaluated and compared with GGMs using both simulated and real RNA-seq data. The proposed method shows improved performance in detecting edges that represent covarying pairs of genes, particularly for edges connecting low-abundant genes and edges around regulatory hubs.

The inference and modeling of network‐like structures in genomic data is of prime importance in systems biology. Complex stochastic associations and interdependencies can very generally be described as a graphical model. However, the paucity of available samples in current high‐throughput experiments renders learning graphical models from genome data, such as microarray expression profiles, a challenging and very hard problem. Here we review several recently developed approaches to small‐sample inference of graphical Gaussian modeling and discuss strategies to cope with the high dimensionality of functional genomics data. Particular emphasis is put on regularization methods and an empirical Bayes network inference procedure.

For many complex heterogeneous diseases, the underlying mechanisms contributing to pathogenesis is still unknown. Although prediction models for complex heterogeneous diseases have been trained on clinical features and imaging data, they may not reveal mechanisms contributing to pathogenesis. The development of human diseases often involves a highly interactive and dynamic system of molecular layers. Networks are a useful framework for representing cross-talk relationships of molecules in biological and disease pathways. Understanding complex biological networks is a fundamental and challenging issue in computational and systems biology. Constructing disease-specific networks is still lacking in many diseases even with current advances. It necessitates the development of methods to integrate protein and genetic interaction information from knowledgebases with omics data to reconstruct networks and develop predictive models to effectively integrate different omics data and network information while being interpretable. Biological networks are often inferred through Gaussian Graphical Models (GGMs) using gene or protein expression data only. GGMs identify conditional dependence by estimating a precision matrix between genes or proteins. However, conventional GGM approaches often ignore prior knowledge about protein-protein interactions (PPI). Recently, several groups have extended GGM to incorporate PPI information. However, these methods are either computationally intractable for large-scale data, or disregard weights in the PPI networks. To address these shortcomings, we developed an augmented high-dimensional graphical lasso (AhGlasso) method to incorporate edge weights in known protein-protein interactions with omics data for global network learning. This new method outperforms weighted graphical lasso-based algorithms with respect to computational time in simulated large-scale data settings while achieving better or comparable prediction accuracy of node connections. Using proteomic data from a study on chronic obstructive pulmonary disease (COPD), we demonstrated that AhGlasso improves protein network inference compared to other approaches that incorporate PPI information. Recently Deep Learning including convolution Graph Neural Network (ConvGNN) has shown great potential for disease classification by integrating a known PPI network from existing databases and a single omics data set. We extended the existing ConvGNN framework to incorporate PPI in addition to multi-omics data to develop a prediction model for COPD. This model achieved significantly better accuracy compared to several conventional classification methods. We also demonstrated that the updated COPD-associated network with AhGlasso further improves the prediction accuracy. Moreover, we applied explainable machine learning methods for interpretation and identified significant genes/proteins and important subnetworks in the ConvGNN model for COPD prediction. Our study suggests the benefits of using ConvGNN for the prediction of heterogeneous complex diseases by integrating multi-omics data and interaction networks.

We propose Bayesian methods for Gaussian graphical models that lead to sparse and adaptively shrunk estimators of the precision (inverse covariance) matrix. Our methods are based on lasso-type regularization priors leading to parsimonious parameterization of the precision matrix, which is essential in several applications involving learning relationships among the variables. In this context, we introduce a novel type of selection prior that develops a sparse structure on the precision matrix by making most of the elements exactly zero, in addition to ensuring positive definiteness - thus conducting model selection and estimation simultaneously. More importantly, we extend these methods to analyze clustered data using finite mixtures of Gaussian graphical model and infinite mixtures of Gaussian graphical models. We discuss appropriate posterior simulation schemes to implement posterior inference in the proposed models, including the evaluation of normalizing constants that are functions of parameters of interest, which result from the restriction of positive definiteness on the correlation matrix. We evaluate the operating characteristics of our method via several simulations and demonstrate the application to real data examples in genomics.

Graph learning for latent-variable Gaussian graphical models under laplacian constraints