Bayesian regularization of Gaussian graphical models with measurement error

1. Introduction

In biomedical settings, it is often of interest to use gene expression microarray data to construct a gene regulatory network for metabolic processes (Segal et al., 2003). In such studies, the gene measurements follow a multivariate Gaussian distribution, where the inverse covariance matrix, known as the precision matrix, characterizes conditional dependence between two genes (called two dimensions, or two features). This is accomplished by noting that if an element of the precision matrix is 0, then the two dimensions are conditionally independent; see Lauritzen (1996) for a review. This setting, often referred to as a Gaussian graphical model, is where our analysis takes place.

Estimating the precision matrix is a difficult task when the number of observations $n$ is often much less than the dimension of the features $d$. A naïve approach is to estimate the precision matrix by the inverse of the empirical covariance matrix; this estimate, however, is known to perform poorly and is ill-posed when $n<d$ (Johnstone et al., 2001). The common approach is to assume that the precision matrix is sparse (Dempster, 1972); that is, we assume the precision matrix’s off-diagonal elements are mostly 0. As a result, most pairs of variables are conditionally independent. The sparsity assumption has led to different lines of research with regularized models to estimate the precision matrix. While some approaches utilize a sparse regression technique that estimates the precision by iteratively regressing each variable on the remaining variables, for instance Khare et al. (2015), we instead focus on the direct likelihood approach. The direct likelihood approach optimizes the full likelihood function with an element-wise penalty on the precision matrix; common examples being graphical lasso (Friedman et al., 2008), CLIME (Cai et al., 2011), and TIGER (Liu et al., 2017). We utilize a recent Bayesian optimization procedure, called BAGUS, that relies on optimization performed by the EM-algorithm, which was shown to have desirable theoretical properties, including consistent estimation of the precision matrix and selection consistency of the conditional pair-wise relationships (Gan et al., 2018).

There are many practical issues associated with Gaussian graphical models, such as hyperparameter tuning (Yuan and Lin, 2007), missing data (Liang et al., 2018), and repeated trials (Tan et al., 2016), which practitioners need to adjust for a successful analysis. We address another practical issue that is often involved with the microarray studies, measurement error. In fact, microarray studies tend to have noisy measurements because of technical variations resulting from sources such as sample preparation, labeling, and hybridization (Zakharkin et al., 2005). In other words, the observations are expression values that have been additionally perturbed with noise from some measurement process. The effects of measurement error on statistical models have been studied extensively for classical settings such as density deconvolution and regression (Carroll et al., 2006), but, to our knowledge, has not yet been well studied in the context of Gaussian graphical models, especially in high dimensional settings.

We propose a Bayesian estimator to correct for measurement error in estimating a sparse precision matrix; our new method extends the optimization procedure of Gan et al. (2018), referred to as BAGUS. While directly incorporating the estimate of the uncontaminated variable is possible, we find the incorporation of the imputation–regularization technique of Liang et al. (2018) to provide more desirable results. Our procedure imputes the mismeasured random variables, then performs BAGUS on this imputation; these steps are iterated for a small number of cycles, requiring more computation but giving better results than the naïve estimator. We prove consistency of the estimated precision matrix with the imputed procedure, and illustrate the performance in a simulation study. Finally, we apply the methodology to a microarray dataset that contains gene expression measurements of favorable histology Wilms tumor (Huang et al., 2009).

2. Contaminated Gaussian graphical models

Given a $d$-dimensional random vector, $\mathit{x}=\{x^1,\dots ,x^d\}$, we are interested in the conditional dependence of two variables $x^i$ and $x^j$, for any pair $(i,j)$ with $1\le i<j\le d$, given all the remaining variables. This conditional dependence structure is usually represented by an undirected graph $G=(V,E)$, where $V=\{1,\dots ,d\}$ is the set of nodes and $E\subseteq V\times V=\{(1,1),(1,2),\dots ,(d,d)\}$ is the set of edges (Lauritzen, 1996). In this representation, the two variables $x^i$ and $x^j$ are conditionally independent if there is no edge between node $i$ and node $j$.

If the vector $\mathit{x}$ follows a multivariate normal distribution with mean $\mathbf{0}$ and covariance matrix ${\mathbf{\Sigma }}_x$, $\mathit{x}\sim N_d(\mathbf{0},{\mathbf{\Sigma }}_x)$, every edge corresponds to a non-zero entry in the precision matrix ${\mathbf{\Omega }}_x={\mathbf{\Sigma }}_x^{-1}$, see Lauritzen (1996). The model in this scenario is often known as a Gaussian graphical model. In the high dimensional setting, the set of edges are usually assumed to be sparse, meaning that only a few pairs $(x^i,x^j)$ are conditionally dependent. In the Gaussian case, this assumption implies only a few off-diagonal entries of ${\mathbf{\Omega }}_x$ are non-zero.

Suppose the data consist of iid observations ${\mathit{w}}_1,\dots ,{\mathit{w}}_n$, where ${\mathit{w}}_i={\mathit{x}}_i+{\mathit{u}}_i,i=1,\dots ,n$ with ${\mathit{x}}_i\sim N_d(\mathbf{0},{\mathbf{\Sigma }}_x)$ and ${\mathit{u}}_i\sim N_d(\mathbf{0},{\mathbf{\Sigma }}_u)$. Here, ${\mathit{w}}_i=(w_i^1,\dots ,w_i^d)$, with the subscript and superscript denoting the observation and components respectively. Denote $\mathit{W}$ as the $n\times d$ matrix of observed data. The model is equivalent to the following hierarchical representation. First, the latent random variables ${\mathit{x}}_i$ are generated from a $N_d(\mathbf{0},{\mathbf{\Sigma }}_x)$ distribution, and when conditioned on ${\mathit{x}}_i$ and ${\mathbf{\Sigma }}_u$, we have ${\mathit{w}}_i|{\mathit{x}}_i,{\mathbf{\Sigma }}_u\sim N_d({\mathit{x}}_i,{\mathbf{\Sigma }}_u)$ for each $i=1,\dots ,n$. This forms an intuitive generative process, where first $\mathit{x}$ is realized, then contaminated by measurement error $\mathit{u}$, and observed finally as $\mathit{w}$. The problem of interest is to estimate the precision matrix ${\mathbf{\Omega }}_x$ in the high dimensional setting $n<d$.

Consider an additive measurement error model where $\mathit{w}=\mathit{x}+\mathit{u}$ and $\mathit{w}$ is the observed data. Denote $\mathit{U}={({\mathit{u}}_1,\dots ,{\mathit{u}}_n)}^T$ as measurement errors that are independent from data $\mathit{X}={({\mathit{x}}_1,\dots ,{\mathit{x}}_n)}^T$. For $i=1,\dots ,n$, the amount of measurement error is drawn from another multivariate normal distribution with mean $\mathbf{0}$ and covariance matrix ${\mathbf{\Sigma }}_{\mathit{u}}$, ${\mathit{u}}_i\sim N_d(\mathbf{0},{\mathbf{\Sigma }}_{\mathit{u}})$. We assume ${\mathbf{\Sigma }}_{\mathit{u}}$ to be diagonal, and hence the amount of measurement error on each variable is uncorrelated. We also assume that ${\mathbf{\Sigma }}_{\mathit{u}}$ is known or estimable from ancillary data, such as replicate measurements. This is a common occurrence, such as in microarray studies where multiple replicates of one subject are collected (Nghiem and Potgieter, 2018), and we illustrate the estimation procedure in Section 5. The contaminated variables $\mathit{w}$ in general have a different conditional dependence structure from that of $\mathit{x}$. Indeed, the covariance and precision matrix of $\mathit{w}$ is given by ${\mathbf{\Sigma }}_w={\mathbf{\Sigma }}_{\mathit{x}}+{\mathbf{\Sigma }}_u$and (1)${\mathbf{\Omega }}_w={\mathbf{\Sigma }}_w^{-1}={\left.{\mathbf{\Sigma }}_{\mathit{x}}+{\mathbf{\Sigma }}_u\right)}^{-1}={\mathbf{\Omega }}_x-{\mathbf{\Omega }}_x{(\mathit{I}+{\mathbf{\Sigma }}_u{\mathbf{\Omega }}_x)}^{-1}{\mathbf{\Sigma }}_u{\mathbf{\Omega }}_x,$respectively; here, $\mathit{I}$ denotes the $d\times d$ identity matrix. Eq. (1) follows from the Kailath variant formula in Petersen et al. (2008). Furthermore, Eq. (1) suggests that ${\mathbf{\Omega }}_w$ and ${\mathbf{\Omega }}_x$ are equal if the product ${\mathbf{\Omega }}_x{(\mathit{I}+{\mathbf{\Sigma }}_u{\mathbf{\Omega }}_x)}^{-1}{\mathbf{\Sigma }}_u{\mathbf{\Omega }}_x$ is equal to a zero matrix. This is generally not the case when the matrix ${\mathbf{\Sigma }}_u$ is not zero.

When no measurement error is present, i.e. the ${\mathit{x}}_i$ are directly observed, the sample covariance matrix $\mathit{S}=n^{-1}{\sum }_{i=1}^n({\mathit{x}}_i-\bar{\mathit{x}}){({\mathit{x}}_i-\bar{\mathit{x}})}^{\top }$, with $\bar{\mathit{x}}$ being the sample mean, is a consistent estimator for ${\mathbf{\Sigma }}_x$. However it has the rank of at most $n<d$, so it is not invertible to estimate ${\mathbf{\Omega }}_x$. When measurement error is present, we assume the covariance matrix of measurement error ${\mathbf{\Sigma }}_u$ is known or estimable from replicates. A naïve approach is first to estimate ${\mathbf{\Sigma }}_x$ by ${\tilde{\mathbf{\Sigma }}}_x={\mathit{S}}_w-{\mathbf{\Sigma }}_u$, where ${\mathit{S}}_w$ denotes the sample covariance from contaminated data $\mathit{W}$, and then to invert ${\tilde{\mathbf{\Sigma }}}_x$ to estimate ${\mathbf{\Omega }}_x$. The main issue with this approach is that ${\tilde{\mathbf{\Sigma }}}_x$ is generally not positive definite. This implies its inverse is also not positive definite, which is necessary to find a consistent estimate ${\mathbf{\Omega }}_x$. Hence, a correction procedure to estimate ${\mathbf{\Omega }}_{\mathit{x}}$ need not rely upon the sample covariance matrix ${\tilde{\mathbf{\Sigma }}}_x$ directly. Furthermore, the procedure should also be able to incorporate sparsity constraints to recover the graphical model structure. These requirements are addressed by the procedure described in the next section.

3. The IRO-BAGUS algorithm

In a recent work, Liang et al. (2018) develop a methodology to efficiently handle high dimensional problems with missing data. Their solution is an EM-algorithm variant which alternates between two steps, the imputation step and regularized optimization step; we refer to their algorithm as the IRO algorithm. Denote the missing data as $Y$, and observed data as $X$. Also denote the desired parameter to be estimated by $\theta$, and begin with some initial guess ${\theta }^{\left(0\right)}$. During the $t$th iteration, the IRO algorithm generates Y from the distribution given by the current estimate of $\theta$, i.e. $Y\sim \pi \left(Y\right|X,{\theta }^{(t-1)})$. Then, using $X$ and $Y$, maximizes $\theta$, under regulation, using the full likelihood. Liang et al. (2018) show that this procedure results in a consistent estimate of ${\theta }^{\left(t\right)}$, and results in a Markov chain with stationary distribution.

We make use of this framework for our current problem pertaining to mismeasured observations instead of missing values. The problems are naturally related in the sense that both are generating values of the true process from some estimated underlying distribution. We return to the hierarchical structure of the problem, i.e. $\mathit{w}|\mathit{x},{\mathbf{\Sigma }}_u\sim N_d(\mathit{x},{\mathbf{\Sigma }}_{\mathit{u}})$ and $\mathit{x}\sim N_d(\mathbf{0},{\mathbf{\Sigma }}_{\mathit{x}})$. The IRO algorithm proceeds iteratively between the two following steps:

• •

  Imputation step: At iteration $t$, draw ${\mathit{X}}^{\left(t\right)}=({\mathit{x}}_1^{\left(t\right)},\dots ,{\mathit{x}}_d^{\left(t\right)})$ from the posterior full conditional of $\mathit{X}$, using the current estimate of ${\mathbf{\Omega }}_{\mathit{x}}^{(t-1)}$. Specifically, for $i=1,\dots ,n$, draw ${\mathit{x}}_i^{\left(t\right)}|\mathit{w},{\mathbf{\Omega }}_{\mathit{x}}^{(\mathit{t}-\mathbf{1})}\sim N_d({\mathbf{\Lambda }}^{-\mathbf{1}}{\mathbf{\Omega }}_{\mathit{u}}{\mathit{w}}_i,{\mathbf{\Lambda }}^{-1})$ where $\mathbf{\Lambda }=({\mathbf{\Omega }}_{\mathit{x}}^{(t-1)}+{\mathbf{\Omega }}_{\mathit{u}})$. Note that the posterior distribution of ${\mathit{x}}_i$ depends only on ${\mathit{w}}_i$ due to independence. This allows for easy generation of data from the true underlying distribution.

• •

  Regularization Step: Apply a regularization to the generated ${\mathit{X}}^{\left(t\right)}$ and obtain a new estimate of ${\mathbf{\Omega }}_{\mathit{x}}^{\left(t\right)}$.

In this work, the regularization step is carried out based on a recent Bayesian methodology, called BAGUS. Hence, the whole algorithm is referred to as the IRO-BAGUS algorithm. We briefly note that the algorithm proposed relies on the, typically, unknown measurement error covariance ${\mathbf{\Sigma }}_u$. Assume replicates are observed and an independence assumption is made between the variables, then the diagonal of ${\mathbf{\Sigma }}_u$ is estimable, as elaborated in Appendix. The next Sections 3.1 The spike-and-slab lasso prior specification, 3.2 The full model, 3.3 Variable selection outline prior specifications, the full model, and variable selection for BAGUS. After that, Section 3.4 discusses consistency of the IRO-BAGUS estimate.

4. Simulation study

4.1. Simulation setup

We investigate the performance of our methodology under several different settings. For each setting we generate ${\mathit{x}}_i$ following a $d$-variate Gaussian distribution with mean $\mathbf{0}$ and precision matrix ${\mathbf{\Omega }}_x$ according to some graphical structure; we refer to this as the true data. Then, the contaminated observations ${\mathit{w}}_i$ were generated from ${\mathit{w}}_i={\mathit{x}}_i+{\mathit{u}}_i$, where ${\mathit{u}}_i\sim N_d(\mathbf{0},{\mathbf{\Sigma }}_u),i=1,\dots ,n$. The measurement error covariance matrix ${\mathbf{\Sigma }}_u$ is assumed to be a diagonal matrix, with element ${\left.{\mathbf{\Sigma }}_u\right]}_{ii}=\gamma {\left.{\mathbf{\Sigma }}_x\right]}_{ii}$, where ${\left.{\mathbf{\Sigma }}_x\right]}_{ii}$ is the variance of the dimension $x^i$. In other words, the constant $\gamma$ controls the noise-to-signal ratio on each variable. For the purposes of simulation, we assume the amount of measurement error to be known.

To generate the true data we use the huge package (Zhao et al., 2012). We inspect two different types of graphs, referred to as hub, and random; we expand on these below where ${\omega }_{ij}$ denotes the $(i,j)$ element of ${\mathbf{\Omega }}_x$.

• 1.

  Hub: The $d$ variables are divided into $K=d∕20$ non-overlapping groups, each group having 20 elements. The $k$th group has a “center” at ${\mathit{X}}_{20(k-1)+1}$, which has connection to every other element in that group, $k=1,\dots ,K$. In other words, ${\omega }_{ij}={\omega }_{ji}=1$ if and only if $i=20(k-1)+1,j=i+m,1\le m\le 20,k=1,\dots ,K$. With that structure, the graph only has $d-K$ edges out of $d(d-1)∕2$ possible edges.

• 2.

  Random: For $1\le i<j\le d$, ${\omega }_{ij}=1$ with probability $\frac{3}{d}$, 0 otherwise.

We illustrate the structures in Fig. 1.

Each model was generated with $n=100$ observations. We inspect each model for $d=\{100,200\}$ and $\gamma =\{0.1,0.25,0.5\}$. The amount of correction-imputations was set to be 50, with the first 20% discarded as burn-in; we note that we inspected 25 and 100 imputations with the same percentage of burn-in samples with minimal differences in output. Each setting was replicated 50 times, and the final results are the average of these replicates. We assume the measurement error to be known for the correction procedure. Hyperparameter tuning was done as described in Section 3.5. Because measurement error is often ignored in the context of GGMs, our simulations also provide perspective onto the negative effect that measurement error can impose on model performance.

1. Download : Download high-res image (304KB)
2. Download : Download full-size image

Fig. 1. Graphical representation for $d=100$ of the hub (left) and random (right) structure, respectively. While the hub structure is fixed for a given $d$, the random graph is subject to change due to the generation process.

To inspect model performance, we examine both the estimated precision matrix and the ability to do variable selection using BAGUS on the true data (true), BAGUS on the contaminated data (naïve), and our IRO-BAGUS methodology on the contaminated data (corrected). For each estimated precision matrix ${\hat{\mathbf{\Omega }}}_x$, estimation error is measured by ${\Vert {\hat{\mathbf{\Omega }}}_x-{\mathbf{\Omega }}_x\Vert }_F$, and variable selection, evaluated by different metrics involving the true positives (TP), false positives (FP), true negatives (TN), and false negatives (FN), are reported: specificity (SPE), sensitivity (SEN), precision (PRE), accuracy (ACC), and Matthews correlation coefficient (MCC); these values are defined as $\text{SPE}=\frac{\text{TN}}{\text{TN + FP}},\text{SEN}=\frac{\text{TP}}{\text{TP + FN}},$$\text{PRE}=\frac{\text{TP}}{\text{TP + FP}},\text{ACC}=\frac{\text{TP + TN}}{\text{TP + FP + TN + FN}},$$\text{MCC}=\frac{\text{TP}\times \text{TN}-\text{FP}\times \text{FN}}{\text{(TP + FP)(TP + FN)(TN + FP)(TN + FN)}}.$Additionally, we also report the area under the ROC curve (AUC). Since the AUC marginalizes over all possible thresholds and tuning parameters, it gives us insight into the amount of separation of the non-zero and zero elements for each method. This reflects the overall ability to recover the graph structure. These different metrics give insight into the tradeoffs and gains of each setting.