Bayesian Regularization of Gaussian Graphical Models with Measurement Error

TL;DR

This paper introduces a Bayesian method to accurately estimate Gaussian graphical models in high-dimensional data contaminated with measurement error, improving the identification of variable relationships.

Contribution

It presents a novel Bayesian approach combining spike-and-slab Lasso and Imputation-Regularization Optimization to correct measurement error effects in precision matrix estimation.

Findings

Outperforms naive methods ignoring measurement error.

Provides more accurate identification of variable relationships.

Successfully applied to gene network inference.

Abstract

We consider a framework for determining and estimating the conditional pairwise relationships of variables when the observed samples are contaminated with measurement error in high dimensional settings. Assuming the true underlying variables follow a multivariate Gaussian distribution, if no measurement error is present, this problem is often solved by estimating the precision matrix under sparsity constraints. However, when measurement error is present, not correcting for it leads to inconsistent estimates of the precision matrix and poor identification of relationships. We propose a new Bayesian methodology to correct for the measurement error from the observed samples. This Bayesian procedure utilizes a recent variant of the spike-and-slab Lasso to obtain a point estimate of the precision matrix, and corrects for the contamination via the recently proposed Imputation-Regularization Optimization procedure designed for missing data. Our method is shown to perform better than the naive method that ignores measurement error in both identification and estimation accuracy. To show the utility of the method, we apply the new method to establish a conditional gene network from a microarray dataset.

Tables

Table 1. Table 1 : Simulation results for the hub graph structure, as specified in Section 5.1 . For each signal-to-noise ratio and d 𝑑 d , the true, naive, and corrected models are shown for metrics defined in Section 5.1 .

$\gamma$	d	Model	SEN	SPE	PRE	ACC	MCC	FROB	AUC
0.1	100	True	1.00	0.65	0.85	0.99	0.73	5.11	0.95
		Naive	1.00	0.50	0.76	0.99	0.61	6.81	0.93
		Corrected	1.00	0.51	0.78	0.99	0.62	6.17	0.97
	200	True	1.00	0.67	0.77	0.99	0.71	7.36	0.94
		Naive	1.00	0.51	0.69	0.99	0.59	9.63	0.92
		Corrected	1.00	0.51	0.71	0.99	0.60	8.68	0.96
0.25	100	True	1.00	0.66	0.84	0.99	0.74	5.09	0.95
		Naive	0.99	0.38	0.60	0.98	0.47	8.54	0.90
		Corrected	1.00	0.36	0.68	0.98	0.49	7.71	0.94
	200	True	1.00	0.67	0.78	1.00	0.72	7.29	0.94
		Naive	1.00	0.40	0.52	0.99	0.45	12.10	0.90
		Corrected	1.00	0.37	0.62	0.99	0.48	10.68	0.95
0.5	100	True	1.00	0.66	0.85	0.99	0.74	5.03	0.95
		Naive	1.00	0.20	0.50	0.98	0.31	9.67	0.84
		Corrected	1.00	0.20	0.70	0.98	0.37	8.74	0.89
	200	True	1.00	0.68	0.77	0.99	0.72	7.36	0.94
		Naive	1.00	0.20	0.37	0.99	0.27	13.70	0.83
		Corrected	1.00	0.17	0.59	0.99	0.31	12.53	0.89

Table 2. Table 2 : Simulation results for the random graph structure, as specified in Section 5.1 . For each signal-to-noise ratio and d 𝑑 d , the true, naive, and corrected models are shown for metrics defined in Section 5.1 .

Amt. ME	d	Model	SEN	SPE	PRE	ACC	MCC	FROB	AUC
0.1	100	True	1.00	0.42	0.84	0.98	0.59	4.61	0.89
		Naive	1.00	0.32	0.80	0.98	0.50	5.34	0.88
		Corrected	1.00	0.32	0.80	0.98	0.50	5.04	0.91
	200	True	1.00	0.36	0.76	0.99	0.52	6.72	0.86
		Naive	1.00	0.30	0.66	0.99	0.44	7.61	0.85
		Corrected	1.00	0.28	0.68	0.99	0.43	7.25	0.90
0.25	100	True	1.00	0.45	0.86	0.98	0.61	4.66	0.90
		Naive	1.00	0.27	0.70	0.97	0.42	6.68	0.85
		Corrected	1.00	0.23	0.75	0.97	0.40	5.72	0.88
	200	True	1.00	0.37	0.75	0.99	0.52	6.77	0.86
		Naive	1.00	0.23	0.52	0.99	0.34	9.04	0.82
		Corrected	1.00	0.18	0.60	0.99	0.32	7.95	0.86
0.5	100	True	1.00	0.43	0.85	0.98	0.59	4.65	0.89
		Naive	1.00	0.14	0.55	0.97	0.26	7.71	0.79
		Corrected	1.00	0.09	0.67	0.97	0.24	6.49	0.79
	200	True	1.00	0.37	0.76	0.99	0.53	6.74	0.86
		Naive	1.00	0.12	0.39	0.98	0.21	10.42	0.77
		Corrected	1.00	0.06	0.56	0.99	0.18	8.92	0.78

Equations

$${\mathbf{\Sigma }}_w={\mathbf{\Sigma }}_{\mathbf{x}}+{\mathbf{\Sigma }}_u$$

$${\mathbf{\Omega }}_w={\mathbf{\Sigma }}_w^{-1}={\left({\mathbf{\Sigma }}_{\mathbf{x}}+{\mathbf{\Sigma }}_u\right)}^{-1}={\mathbf{\Omega }}_x-{\mathbf{\Omega }}_x(\mathbf{I}+{\mathbf{\Sigma }}_u{\mathbf{\Omega }}_x{)}^{-1}{\mathbf{\Sigma }}_u{\mathbf{\Omega }}_x,$$

$$\pi ({\omega }_{ij})=\frac{\eta }{2v_1}\exp \left\{-\frac{|{\omega }_{ij}|}{v_1}\right\}+\frac{1-\eta }{2v_0}\exp \left\{-\frac{|{\omega }_{ij}|}{v_0}\right\}$$

$$\pi ({\mathbf{\Omega }}_x)=\prod _{i<j}\pi ({\omega }_{ij})\prod _i\pi ({\omega }_{ii})I({\mathbf{\Omega }}_x\succ 0)I(||{\mathbf{\Omega }}_x||\le B).$$

$$\pi ({\mathbf{\Omega }}_x|\mathbf{X})\propto \prod _{i=1}^n\pi ({\mathbf{x}}_i|{\mathbf{\Omega }}_x)\pi ({\mathbf{\Omega }}_x).$$

$$L^{\text{UC}}({\mathbf{\Omega }}_x)=\log \pi ({\mathbf{\Omega }}_x|\mathbf{X})=\frac{n}{2}\left(\text{tr}({\mathbf{X}}^T{\mathbf{\Omega }}_x\mathbf{X})-\text{logdet}({\mathbf{\Omega }}_x)\right)+\sum _{i<j}\pi ({\omega }_{ij})+\sum _i\pi ({\omega }_{ii})+K$$

$$\{\begin{array}{l}{\omega }_{ij}|r_{ij}=0\sim \text{Laplace}(0,v_0)\\ {\omega }_{ij}|r_{ij}=1\sim \text{Laplace}(0,v_1)\end{array}.$$

$$p_{ij}=\frac{v_1}{v_0}\frac{1-\eta }{\eta }\exp \left\{|{\omega }_{ij}|\left(\frac{1}{v_1}-\frac{1}{v_0}\right)\right\}.$$

$${\mathbf{w}}_i|{\mathbf{x}}_i,{\mathbf{\Omega }}_x\sim N_d({\mathbf{x}}_i,{\mathbf{\Sigma }}_u),i=1,\dots ,n$$

$${\mathbf{x}}_i|{\mathbf{\Omega }}_x\sim N_d(0,{\mathbf{\Omega }}_x^{-1}),i=1,\dots ,n$$

$${\omega }_{ij}|r_{ij}=0,v_0\sim \text{Laplace}(0,v_0),i\ne j,i,j=1,\dots ,n$$

$${\omega }_{ij}|r_{ij}=1,v_1\sim \text{Laplace}(0,v_1),i\ne j,i,j=1,\dots ,n$$

$${\omega }_{ii}\sim \text{Exp}(\tau ),i=1,\dots ,n$$

$$r_{ij}|\eta \sim \text{Bern}(\eta ),i\ne j,!i,j=1,\dots ,n.$$

$$|{\mathbf{\Omega }}_x{|}^{\frac{1}{2}}\exp \left\{-\frac{1}{2}{\mathbf{X}}^T{\mathbf{\Omega }}_x\mathbf{X}\right\}\prod _{i<j}\pi ({\omega }_{ij}|r_{ij})\pi (r_{ij}|\eta )\prod _i\pi ({\omega }_{ii}|\tau )I({\mathbf{\Omega }}_x\succ 0)I(||{\mathbf{\Omega }}_x||\le B).$$

$$Q({\mathbf{\Omega }}_x|{\mathbf{\Omega }}_x^{(t)})={\mathbb{E}}_{\mathbf{R}|{\mathbf{\Omega }}_x^{(t)}}\log \pi ({\mathbf{\Omega }}_x,\mathbf{X}|\mathbf{W},{\mathbf{\Sigma }}_u),$$

$$\text{BIC}=n(\text{tr}(\mathbf{S}{\hat{\mathbf{\Omega }}}_x)-\text{logdet}({\hat{\mathbf{\Omega }}}_x))+\log (n)\times q,$$

$$\begin{array}{rl}& \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)}}.\end{array}$$

$$\begin{array}{rl}\text{(a) }& \frac{1}{n{v}_1}=a_1\sqrt{\frac{\log (d)}{n}}(1-{ϵ}_0),\\ \text{(b) }& \frac{1}{n{v}_0}>a_5\sqrt{\frac{\log (d)}{n}},\\ \text{(c) }& \frac{v_1^2(1-\eta )}{v_0^2\eta }\le d^{{ϵ}_1},\\ \text{(d) }& \tau \le a_1\frac{n}{2}\sqrt{\frac{\log (d)}{n}}.\end{array}$$

$$\frac{1}{k_1}+2b(a_1+a_4)M_{\Gamma }\sqrt{\frac{\log (d)}{n}}<B<\sqrt{2n{v}_0}.$$

$$\log \pi (\mathbf{X},\mathbf{W}|{\mathbf{\Omega }}_x)=\log \pi (\mathbf{X}|{\mathbf{\Omega }}_x)+\log \pi (\mathbf{W}|\mathbf{X},{\mathbf{\Omega }}_u).$$

$$\underset{n\ge 1}{\sup }\underset{\mathbf{X},{\mathbf{\Omega }}_x}{\sup }|\int m_n(\mathbf{X},\mathbf{W})\mathbb{I}(m_n(\mathbf{X},\mathbf{W})>\xi )\pi (\mathbf{X}|\mathbf{W},{\mathbf{\Omega }}_x)|\to 0.$$

$$\begin{array}{rl}\log \pi (\mathbf{X},\mathbf{W}|{\mathbf{\Omega }}_x)& =\sum _{i=1}^n\left[\log \pi ({\mathbf{w}}_i|{\mathbf{x}}_i,{\mathbf{\Omega }}_u)+\log \pi ({\mathbf{x}}_i|{\mathbf{\Omega }}_x)\right]\\ & =-\frac{1}{2}\sum _{i=1}^n\left[({\mathbf{w}}_i-{\mathbf{x}}_i{)}^T{\mathbf{\Omega }}_u({\mathbf{w}}_i-{\mathbf{x}}_i)+{\mathbf{x}}_i^T{\mathbf{\Omega }}_x{\mathbf{x}}_i\right]+\frac{1}{2}\log \text{det}({\mathbf{\Omega }}_x)+C,\end{array}$$

$$\begin{array}{rl}|\log \pi (\mathbf{X},\mathbf{W}|{\mathbf{\Omega }}_x)|& \le \frac{1}{2}\sum _{i=1}^n\left[({\mathbf{w}}_i-{\mathbf{x}}_i{)}^T{\mathbf{\Omega }}_u({\mathbf{w}}_i-{\mathbf{x}}_i)+K_1{\mathbf{x}}_i^T{\mathbf{x}}_i\right]+K_2\\ & =\sum _{i=1}^{n}m({\mathbf{x}}_i,{\mathbf{w}}_i)=m(\mathbf{X},\mathbf{W}),\end{array}$$

$$\begin{array}{rl}\tilde{m}(\mathbf{W},{\mathbf{\Omega }}_x^{(t)})& =\int m(\mathbf{X},\mathbf{W})\pi (\mathbf{X}|\mathbf{W},{\mathbf{\Omega }}_x^{(t)})d\mathbf{X}\\ & =\int \sum _{i=1}^{n}m({\mathbf{x}}_i,{\mathbf{w}}_i)\left[\prod _{j=1}^n\pi ({\mathbf{x}}_i|{\mathbf{w}}_i,{\mathbf{\Omega }}_x^{(t)})\right]d{\mathbf{x}}_1,\dots ,d{\mathbf{x}}_n\\ & =\sum _{i=1}^n\int m({\mathbf{x}}_i,{\mathbf{w}}_i)\pi ({\mathbf{x}}_i|{\mathbf{w}}_i,{\mathbf{\Omega }}_x^{(t)})d{\mathbf{x}}_i,\end{array}$$

$${\mathbb{E}}_{{\mathbf{x}}_i|{\mathbf{w}}_i,{\mathbf{\Omega }}_x^{(t)}}\left[m({\mathbf{x}}_i,{\mathbf{w}}_i)\right]=\frac{1}{2}{\mathbf{w}}_i^T{\mathbf{\Omega }}_u{\mathbf{w}}_i+\frac{1}{2}\text{tr}(({\mathbf{\Omega }}_u+K_1{\mathbf{I}}_d){\mathbf{\Lambda }}^{(t)})-\underset{\ge 0}{\underbrace{{\mathbf{w}}_i^T{\mathbf{\Omega }}_u{\mathbf{\Lambda }}^{(t)}{\mathbf{\Omega }}_u{\mathbf{w}}_i}},$$

$$\tilde{m}(\mathbf{W},{\mathbf{\Omega }}_x^{(t)}\le \frac{1}{2}\sum _{i=1}^n{\mathbf{w}}_i^T{\mathbf{\Omega }}_u{\mathbf{w}}_i+K_3=m^{\ast }(\mathbf{W}).$$

$$\begin{array}{rl}\log \pi ({\mathbf{x}}_i,{\mathbf{w}}_i|{\mathbf{\Omega }}_x)& =-\frac{1}{2}({\mathbf{w}}_i-{\mathbf{x}}_i{)}^T{\mathbf{\Omega }}_u({\mathbf{w}}_i-{\mathbf{x}}_i)-\frac{1}{2}{\mathbf{x}}_i{\mathbf{\Omega }}_x{\mathbf{x}}_i\\ & =-\frac{1}{2}{\mathbf{x}}_i^T({\mathbf{\Omega }}_x+{\mathbf{\Omega }}_u){\mathbf{x}}_i+{\mathbf{x}}_i^T{\mathbf{\Omega }}_u{\mathbf{w}}_i+C_1,\end{array}$$

$$\begin{split}&\int\log\pi(\bm{x}_{i},\bm{w}_{i}|\bm{\Omega}_{x})\pi(\bm{X}|\bm{w}_{i},\bm{\Omega}_{x}^{(t)})\\ =&\int\big{[}\log\pi(\bm{w}_{i},\bm{x}_{i}|\bm{\Omega}_{x})+\log\pi(\bm{X}_{-i}|\bm{\Omega}_{x})\big{]}\pi(\bm{x}_{i}|\bm{w}_{i},\bm{\Omega}_{x}^{(t)})\pi(\bm{X}_{-i}|\bm{\Omega}_{x}^{(t)})d\bm{x}_{i}d\bm{X}_{-i}\\ =&\underbrace{\int\log\pi(\bm{w}_{i},\bm{x}_{i}|\bm{\Omega}_{x})\pi(\bm{x}_{i}|\bm{w}_{i},\bm{\Omega}_{x}^{(t)})d\bm{x}_{i}}_{=A}\underbrace{\int\pi(\bm{X}_{-i}|\bm{\Omega}_{x}^{(t)})d\bm{X}_{-i}}_{=1}\\ &\qquad+\underbrace{\int\log\pi(\bm{X}_{-i}|\bm{\Omega}_{x}\pi(\bm{X}_{-i}|\bm{\Omega}_{x}^{(t)})d\bm{X}_{-i}}_{=C_{2}}\underbrace{\int\pi(\bm{x}_{i}|\bm{w}_{i},\bm{\Omega}_{x}^{(t)})d\bm{x}_{i}}_{=1}.\end{split}$$

$${\mathbb{E}}_{{\mathbf{x}}_i|{\mathbf{w}}_i,{\mathbf{\Omega }}_x}({\mathbf{x}}_i({\mathbf{\Omega }}_x+{\mathbf{\Omega }}_u){\mathbf{x}}_i)=\text{tr}(({\mathbf{\Omega }}_x+{\mathbf{\Omega }}_u){\Lambda }^{-1,(t)})+{\mathbf{w}}_i{\mathbf{\Omega }}_u{\mathbf{\Lambda }}^{(t)}({\mathbf{\Omega }}_x+{\mathbf{\Omega }}_u){\mathbf{\Lambda }}^{(t)}{\mathbf{\Omega }}_u{\mathbf{w}}_i$$

$${\mathbb{E}}_{{\mathbf{x}}_i|{\mathbf{w}}_i,{\mathbf{\Omega }}_x}({\mathbf{x}}_i^T{\mathbf{\Omega }}_u{\mathbf{w}}_i)={\mathbf{w}}_i^T{\mathbf{\Omega }}_u{\mathbf{\Lambda }}^{(t)}{\mathbf{\Omega }}_u{\mathbf{w}}_i.$$

$$-\frac{1}{2}{\mathbf{x}}_i^T({\mathbf{\Omega }}_x+{\mathbf{\Omega }}_u){\mathbf{x}}_i+{\mathbf{x}}_i^T{\mathbf{\Omega }}_u{\mathbf{w}}_i-\frac{1}{2}{\mathbf{w}}_i^T{\mathbf{\Omega }}_u{\mathbf{\Lambda }}^{(t)}({\mathbf{\Omega }}_x+{\mathbf{\Omega }}_u){\mathbf{\Lambda }}^{(t)}{\mathbf{\Omega }}_u{\mathbf{w}}_i+{\mathbf{w}}_i^T{\mathbf{\Omega }}_u{\mathbf{\Lambda }}^{(t)}{\mathbf{\Omega }}_u{\mathbf{w}}_i+C,$$

$$v_1=\frac{1}{a_1(1-{ϵ}_0)\sqrt{n\log (d)}}\to 0$$

Full text

1 Introduction

A core problem in statistical inference is estimating the conditional relationship among random variables. Naturally, a full description of the underlying connections among the numerous random variables is valuable information across many disciplines, such as in biology where the relationships among hundreds of genes involved in a metabolic process is desired to be uncovered. In fact, under the assumption that the variables follow a multivariate Gaussian distribution, the inverse covariance matrix, known as the precision matrix, characterizes conditional dependence between two dimensions. This is accomplished by noting that if an element of the precision matrix is 0, then the two variables are conditionally independent; see [12] 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 naive 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$ [9]. The common approach is to assume that the precision matrix is sparse [3]; 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 one approach utilizes a sparse regression technique that estimates the precision by iteratively regressing each variable on the remaining variables, for instance [10], 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 [5], CLIME [1], and TIGER [14]. 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 [6].

There are many practical issues associated with Gaussian graphical models, such as hyperparameter tuning [24], missing data [13], and repeated trials [22], which practitioners need to adjust for a successful analysis. We address another practical issue involved with these models, measurement error. Measurement error occurs when the variables of interest are not observed directly; instead, the observations are the desired variables that have been additionally perturbed with noise from some measurement process. This happens when, for instance, an inaccurate device is used to measure some sort of health metric. Measurement error models have been studied extensively for classical settings such as density deconvolution and regression [2], but, to our knowledge, have not yet been well studied in the context of Gaussian graphical models, especially in high dimensional setting.

We propose a Bayesian methodology to correct for measurement error in estimating a sparse precision matrix; our new method extends the optimization procedure of [6]. While directly incorporating the estimate of the uncontaminated variable is possible, we find the incorporation of the imputation-regularization technique of [13] to provide more desirable results. Our procedure imputes the mismeasured random variables, then performs BAGUS on this imputation; these steps are performed for a small number of cycles, requiring more computation but giving better results than the naive estimator. We prove consistency of the estimated precision matrix with the imputed procedure, and illustrate the performance in a simulation study. We conclude with an application to a microarray data.

2 Contaminated Gaussian Graphical Models

Given a $d$-dimensional random vector, $\bm{x}=\{x^{1},\ldots,x^{d}\}$, the conditional dependence of two variables $x^{i}$ and $x^{j}$, for any pair $(i,j)$ with $1\leq i\leq j\leq d$, given all the remaining variables is of interest. This conditional dependence structure is usually represented by an undirected graph $G=(V,E)$, where $V=\{1,\ldots,d\}$ is the set of nodes and $E\subseteq V\times V=\{(1,1),(1,2),\ldots,(d,d)\}$ is the set of edges [12]. 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 $\bm{x}$ follows a multivariate normal distribution with mean $\bm{0}$ and covariance matrix $\bm{\Sigma}_{x}$, $\bm{x}\sim N_{d}(\bm{0},\bm{\Sigma}_{x})$, every edge corresponds to a non-zero entry in the precision matrix $\bm{\Omega}_{x}=\bm{\Sigma}_{x}^{-1}$, see [12]. 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 $\bm{\Omega}_{x}$ are non-zero.

When measurement error is present, denote $\bm{U}=(\bm{u}_{1},\ldots,\bm{u}_{n})^{T}$ as measurement errors that are independent from data $\bm{X}=(\bm{x}_{1},\ldots,\bm{x}_{n})^{T}$. For $i=1,\ldots,n$, the amount of measurement error is drawn from another multivariate normal distribution with mean $\bm{0}$ and covariance matrix $\bm{\Sigma_{u}}$, $\bm{u}_{i}\sim N_{d}(\bm{0},\bm{\Sigma_{u}}).$ We assume $\bm{\Sigma_{u}}$ to be diagonal, and hence the amount of measurement error on each variable is uncorrelated. We make a common assumption that $\bm{\Sigma_{u}}$ is known or estimable from ancillary data, such as replicate measurements. The contaminated variables $\bm{w}=\bm{x}+\bm{u}$ in general have a different conditional dependence structure from that of $\bm{x}$. Indeed, the covariance and precision matrix of $\bm{w}$ is given by

[TABLE]

and

[TABLE]

respectively; here, $\bm{I}$ denotes the $d\times d$ identity matrix. Equation (1) follows from the Kailath variant formula in [17]. Furthermore, equation (1) suggests that $\bm{\Omega}_{w}$ and $\bm{\Omega}_{x}$ are equal if the product $\bm{\Omega}_{x}(\bm{I}+\bm{\Sigma}_{u}\bm{\Omega}_{x})^{-1}\bm{\Sigma}_{u}\bm{\Omega}_{x}$ is equal to a zero matrix. This is generally not the case when the matrix $\bm{\Sigma}_{u}$ is not zero.

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

When no measurement error is present, i.e the $\bm{x}_{i}$ are directly observed, the sample covariance matrix $\bm{S}=n^{-1}\sum_{i=1}^{n}(\bm{x}_{i}-\bar{\bm{x}})(\bm{x}_{i}-\bar{\bm{x}})^{\top}$, with $\bar{\bm{x}}$ being the sample mean, is a consistent estimator for $\bm{\Sigma}_{x}$. However it has the rank of at most $n<d$, so it is not invertible to estimate $\bm{\Omega}_{x}$. When measurement error is present, we assume the covariance matrix of measurement error $\bm{\Sigma}_{u}$ is known or estimable from replicates. A naive approach is first to estimate $\bm{\Sigma}_{x}$ by $\tilde{\bm{\Sigma}}_{x}=\bm{S}_{w}-\bm{\Sigma}_{u}$, where $\bm{S}_{w}$ denotes the sample covariance from contaminated data $\bm{W}$, and then to invert $\tilde{\bm{\Sigma}}_{x}$ to estimate $\bm{\Omega}_{x}$. The main issue with this approach is that $\tilde{\bm{\Sigma}}_{x}$ is generally not positively definite. This implies its inverse is also not positively definite, which is necessary to find a consistent estimate $\bm{\Omega}_{x}$. Hence, a correction procedure to estimate $\bm{\Omega_{x}}$ need not rely upon the sample covariance matrix $\tilde{\bm{\Sigma}}_{x}$ directly. Furthermore, the procedure is also 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, [13] develop a methodology to efficiently handle high dimension 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^{(0)}$. During the $t^{\text{th}}$ iteration, the IRO algorithm generates Y from the distribution given by the current estimate of $\theta$, i.e. $Y\sim\pi(Y|X,\theta^{(t-1)})$. Then, using $X$ and $Y$, maximizes $\theta$, under regulation, using the full likelihood. [13] show that this procedure results in a consistent estimate of $\theta^{(t)}$, 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. $\bm{w}\sim N_{d}(\bm{x},\bm{\Sigma_{u}})$ and $\bm{x}\sim N_{d}(\bm{0},\bm{\Sigma_{x}})$. The IRO algorithm proceeds iteratively between the two following steps:

• •

Imputation step: At iteration $t$, draw $\bm{X}^{(t)}=(\bm{x}_{1}^{(t)},\ldots,\bm{x}_{d}^{(t)})$ from the posterior full conditional of $\bm{X}$, using the current estimate of $\bm{\Omega_{x}}^{(t-1)}$. Specifically, for $i=1,\ldots,n$, draw $\bm{x}_{i}^{(t)}|\bm{w},\bm{\Omega_{x}^{(t-1)}}\sim N_{d}(\bm{\Lambda^{-1}}\bm{\Omega_{u}}\bm{w}_{i},\bm{\Lambda}^{-1})$ where $\bm{\Lambda}=(\bm{\Omega_{x}}^{(t-1)}+\bm{\Omega_{u}})$. Note that the posterior distribution of $\bm{x}_{i}$ depends only on $\bm{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 $\bm{X}^{(t)}$ and obtain a new MAP estimate of $\bm{\Omega_{x}}^{(t)}$.

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. The next subsections 3.1-3.3 outline prior specifications, the full model, and variable selection for BAGUS. After that, section 3.4 discusses consistency of the IRO-BAGUS estimate.

3.1 The Spike-and-Slab Lasso Prior Specification

Denote the elements $\bm{\Omega}_{x}$ to be $\omega_{ij}$. Recently, a non-convex, continuous relaxation penalty for the spike-and-slab prior was created for the standard lasso problem [20]. This prior was extended to the case of graphical models by [6], and is given by

[TABLE]

for the off diagonal elements ($i\neq j$), where $0<v_{0}<v_{1}$ and $0<\eta<1$. This prior can be interpreted as a mixture of the spike-and-slab prior. The first component of the mixture has prior probability $\eta$, and is associated with the slab component, i.e. $\omega_{ij}\neq 0$. Conversely, with prior probability $1-\eta$ the element is from the spike component, suggesting $\omega_{ij}=0$.

Traditionally, the spike-and-slab prior has a point mass component at 0 and some other continuous distribution for the slab component. This is to represent setting unwanted terms exactly to 0. Here, both the spike and the slab components are distributed according to a Laplace distribution; both are centered at 0, but the spike is more tightly centered by a smaller variance term than the slab. This relaxation of the spike-and-slab prior allows for efficient gradient based algorithms, while still being theoretically sound as shown in [19].

Shrinkage is not desired on the diagonal elements, so a weakly informative exponential prior is given instead, $\pi(\omega_{ii})=\tau\exp\left\{-\tau\omega_{ii}\right\}.$ Another consideration for the prior of $\bm{\Omega}_{x}$ is to ensure the whole matrix to be positive definite, denoting as $\bm{\Omega}_{x}\succ 0$. Moreover, in line with [6], we require the spectral norm to be bounded above by some value $B$, $||\bm{\Omega}_{x}||\leq B$. This assumption will be important going forward. The full prior distribution for $\bm{\Omega}_{x}$ is then given by

[TABLE]

3.2 The Full Model

Without measurement error, the posterior distribution is specified as

[TABLE]

The full conditionals can be derived for (4), but, to avoid costly MCMC sampling for this large dimensional problem, [6] opted to instead find the mode of the the posterior distribution, often referred to as the MAP. The MAP can be found by minimizing the uncontaminated (UC) objective

[TABLE]

with respect to $\bm{\Omega}_{x}$, where $K$ is the normalizing constant in (4). To this end, [6] proved the local convexity of (4) when $||\bm{\Omega}_{x}||\leq B<\infty$, which allows an easy optimization procedure that converges asymptotically to the correct precision matrix.

3.3 Variable Selection

Many practictioners use Gaussian graphical models for the purpose of identifying non-zero entries of $\bm{\Omega}_{x}$, which signify conditional dependencies among the two different variables. The spike-and-slab lasso formulation allows for this quite easily by viewing the optimization as an instance of the EM-algorithm and defining the hierarchical prior

[TABLE]

Here, $r_{ij}$ is the random indicator that the element of the precision matrix follows from the spike or the slab component, where $r_{ij}\sim\text{Bern}(\eta)$. A further hierarchical level can be added by treating $\eta$ as random instead of a fixed hyperparameter. Recent work from [4] illustrates this and is line with the spike-and-slab Lasso setting of [20]. Given our purpose is to study the effect of the measurement error, we choose to treat it as a fixed.

The conditional posterior distribution for $r_{ij}$ is also Bernoulli, with probability of success

[TABLE]

We will use the MAP estimate of $\omega_{ij}$ in (7) as the approximate probability of inclusion. With the inclusion probability a hard threshold will be specified to determine the final inclusion for the purpose of simulation and model selection. Denote $\bm{R}$ and $\bm{P}$ to be the matrix of indicators and conditional posterior probability of inclusion for each element of $\bm{\Omega}_{x}$. We note that for final inference it may be better to forego this inclusion threshold, and instead rank-order the $p_{ij}$ for purposes of downstream investigation; however, this will depend on the application at hand.

3.4 Consistency of the IRO-BAGUS algorithm

The entire data generation process for the contaminated sample is summarized below:

[TABLE]

Instead of approximating the posterior distribution of all the parameters, the IRO-BAGUS algorithm iteratively generates realizations of uncontaminated data, $\bm{X}$, then optimizes $\bm{\Omega}_{x}$ with these generated values. Under some technical conditions, the IRO algorithm is shown to produce a consistent estimate after each iteration in the context of missing data when the regularization step results in a consistent estimate [13]. We show that these conditions are also held in the case of contaminated data, so the IRO-BAGUS algorithm results in a consistent estimate. Theorem 1 is the analogue statement of consistency as in the missing data case. The proof is given in the appendix.

Theorem 1**.**

Assuming $||\bm{\Omega}_{x}||\leq B$, then the estimate $\bm{\Omega}_{x}^{(t)}$ is uniformly consistent to $\bm{\Omega}_{x}$ when $\log(t)=\mathcal{O}(n)$.

It can be seen that the nature of the IRO algorithm is similar to that of MCMC. Under a few more conditions, [13] note that the IRO results in a Markov chain with a stationary distribution, and hence the average of the maximization steps are consistent estimates of the underlying parameters. Our final estimates are the averaged regularized optimization steps given by BAGUS from the imputed data at each iteration, removing a small number of the beginning iterations as burn-in. By averaging instead of taking only the final iteration, we make the analysis less variable. In this sense, the relationships that the correction procedure identifies are more likely to be true relationships, cutting down on the number of false positives.

4 Computation for the IRO-BAGUS algorithm

4.1 Finding MAP estimate for $\bm{\Omega}_{x}$

Here we consider some computational aspects of our proposed methodology. First, we focus to the optimization of $\bm{\Omega}_{x}$. In our procedure, once $\bm{X}$ is generated, the objective function to be optimized is $L^{\text{uc}}$, as was shown in Equation (7); we note this is due to the conditional independence of $\bm{W}$ and $\bm{\Omega}_{x}$ in the hierarchical structure of the contamination process. Optimizing $L^{\text{uc}}$ is difficult to do directly; therefore, the latent factors $r_{ij}$ from Section 3.3 are introduced into the process as in [6]. This allows an E-step similar to the spike-and-slab Lasso and an M-step similar to the Graphical Lasso.

The optimization seeks to find the MAP of the posterior proportional to

[TABLE]

The E-step takes the conditional expectation of $r_{ij}$ in the proportional posterior. Each $r_{ij}$ is conditionally Bernoulli with probability as given in Equation (7), which allows for easy calculation of the desired conditional expectations. Then, the desired $Q$ function to maximize in the M-step is given by

[TABLE]

where the expectation is taken element wise for $\bm{R}$ by assumed independence of inclusion. Maximizing $Q$ is done by a block coordinate descent algorithm. The algorithm cycles between column-wise updates of $\Omega_{x}$. We put the details of this procedure in the Appendix.

4.2 Other Computation Considerations

4.2.1 Estimating $\bm{\Sigma}_{u}$

We have assumed the covariance matrix of measurement error $\bm{\Sigma}_{u}$ to be known before applying the IRO-BAGUS algorithm. In practice, the matrix $\bm{\Sigma}_{u}$ is often estimated from ancillary data, such as replicate observations. Assuming measurement error between variables to be independent is reasonable for many problems and often used in the literature [21]. In that case, only the diagonal of $\bm{\Sigma}_{u}$ only needs to be estimated. For the data analysis application we provide in Section 6, we estimated them with the method described in [23], assuming homogeneity of measurement error between observations. After that, we performed the IRO-BAGUS algorithm as previously described.

4.2.2 Starting Values

The starting value plays a significant role in the speed of optimization at each step. To begin, we perform a naive analysis on the raw contaminated data, $\bm{W}$, giving estimate $\bm{\Omega}_{x}^{(0)}$. This value is then used to start the IRO procedure by generating $\bm{X}$. Each optimization has a warm start from the previous iteration’s estimated precision matrix.

4.2.3 Addressing the constraint, $||\bm{\Omega}_{x}||\leq B$

The constaint that $||\bm{\Omega}_{x}||\leq B$ needs to be addressed. [6] suggest using a threshold on the largest absolute value of the elements of the column being updated in the block coordinate descent. We use the same threshold, and find no performance issues when used with the IRO algorithm.

4.2.4 Positive-Definiteness of $\bm{\Omega}_{x}$

Many procedures to estimate a sparse precision can not guarantee postive-definiteness, however [6] show that the output of BAGUS from the EM algorithm is always symmetric and positive definite. It is easy to show that the imputation step, with final results averaged, also results in this nice property.

Theorem 2**.**

The estimated precision matrix $\hat{\bm{\Omega}}_{x}=T^{-1}\sum_{t=1}^{T}\bm{\Omega}_{x}^{(t)}$ is symmetric and positive definite if the initial value of $\bm{\Omega}_{x}$ for BAGUS at each iteration $t$ was also positive definite.

Proof.

By Theorem 5 in [6], if the initial value to optimize BAGUS is positive definite, then $\bm{\Omega}_{x}^{(t)}$ is also positive definite. The set of positive definite matrices form a cone, and hence the average will also be in this cone. ∎

4.2.5 Parameter Tuning

There are four hyperparameters in BAGUS, $\eta,\tau,v_{0},\text{and }v_{1}$. As with [6], we always set $\eta=0.5$ and $\tau=v_{0}$, which leaves two hyperparameters to tune. Again, we follow [6], who suggest a BIC-like criteria to select the best model from a grid of hyperparameters. This criteria is

[TABLE]

where $\hat{\bm{\Omega}}_{x}$ is the estimated precision matrix and $q$ is the number of non-zero elements of the estimated in the upper diagonal of the precision matrix. We use this in similar fashion for the IRO procedure, but instead we use the averaged $\bm{\Omega}_{x}^{(t)}$ in the BIC calculation.

5 Simulation

5.1 Simulation Setup

We investigate the performance of our methodology under several different settings. For each setting we generate $\bm{x}_{i}$ following a $d$-variate Gaussian distribution with mean $\bm{0}$ and precision matrix $\bm{\Omega}_{x}$ according to some graph structure; we refer to this as the true data. Then, the contaminated observation $\bm{w}_{i}$ was generated from $\bm{w}_{i}=\bm{x}_{i}+\bm{u}_{i}$, where $\bm{u}_{i}\sim N_{d}(\bm{0},\bm{\Sigma}_{u}),~{}i=1,\ldots,n$. The measurement error covariance matrix $\bm{\Sigma}_{u}$ is assumed to be a diagonal matrix, with element $\left[\bm{\Sigma}_{u}\right]_{ii}=\gamma\left[\bm{\Sigma}_{x}\right]_{ii}$, where $\left[\bm{\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 [25]. 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 $\bm{\Omega}_{x}$.

Hub: For $d/20$ groups, $\omega_{ij}=\omega_{ji}=1$ if in the same group. $\omega_{ij}=0$ otherwise. 2. 2.

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

We illustrate the stuctures in Figure 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. Hyperparameter tuning was done as described in Section 4.2. 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 into the model performance.

To inspect model performance, we examine both the estimated precision matrix and the ability to do variable selection of BAGUS on the true data (true), BAGUS on the contaminated data (naive), and our IRO-BAGUS methodology on the contaminated data (corrected). For each estimated precision matrix $\hat{\bm{\Omega}}_{x}$, estimation error is measured by $||\hat{\bm{\Omega}}_{x}-\bm{\Omega}_{x}||_{F}$, and variable selection is 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

[TABLE]

Additionally, we also report the area under the ROC curve (AUC), which gives insight into the amount of seperation of the classification. These different metrics give insight into the tradeoffs and gains of each setting.

5.2 Simulation Results

Table 1 and Table 2 present the results for the hub and random structure, respectively. To begin, we note the effect of the increasing measurement error. This can be observed by examining the growing difference in the performance of the true and naive model when holding $d$ fixed and increasing the amount of contamination. Focusing on the hub structure, a decrease in the quality of selection and estimation can be observed for each setting, which grows worse with more contamination. The selection accuracy metrics with respect to the prespecified 0.5 cut-off show drops in performance of around $50\%$. The estimated precision matrix from the naive grows worse with measurement error, and is also about $50\%$ worse when the signal-to-noise is 0.5.

We now turn attention to the performance of the correction step. First, take note of the first five metrics which are based on the confusion matrix for the 0.5 cutoff threshold. Averaging across the IRO iterations was expected to result in an analysis that favored identifying relationships that were more certain, which can be observed by inspection of the precision (PRE). The gains from the precision are most notable as $d$ grows larger, and more pair-wise relationships exist; when $d=200$, we note nearly 10% and 50% performance gains in the precision for signal-to-noise ratios of 0.25 and 0.5, respectively. In both the hub and random structure the naive and corrected models perform similarly in terms of the sensitivity, specificity, accuracy, and MCC.

It seems at first glance that the selection performance, ignoring the precision, of the correction procedure is comparable to the naive, but these discrepencies can be attributed to the prespecified inclusion cut-off on the $\bm{P}$ matrix. In practice it can often be more reasonable to rank order the inclusion probabilities to identify relationships to further investigate in future experiments. With this in mind, we turn to the performance with respect to the AUC where consistent improvements can be seen for the hub and random structure in most all settings. The AUC helps understand the amount of seperation found in the model across all thresholds, which helps justify that the correction step is making improvements in seperating the classes for the true relationships as AUC improvements are seen in all but the random graph with $d=200$ and signal-to-noise ratio of $0.5$.

We note two items in regard to the AUC. First, the AUC of the corrected model sometimes outperforms the true model, too. In particular, this happens in the hub structure when the amount of measurement error is $0.1$. This can be attributed to the measurement error in models that are easily identified. Second, in the random structure with the amount of measurement error being $0.5$, the corrected model does not make substantial improvements in results over the naive model. We note the difficulty of this setting, as the random structure often performs worse than other structures in identification, and now we add more noise via the contamination. With a relatively small sample size, this noise is difficult to overcome.

Finally, we note the quality of the estimated precision matrix, as measured by Frobenious norm of the difference. In every setting for both the hub and random matrices, the corrected model outperforms the naive model’s estimate of the precision matrix. In the hub structure this improvement is often of the order of $15$-$20\%$ better, while in the random structure a $10\%$ improvement is generally observed. If the intent of the analysis is to use the estimated precision matrix in downstream analysis, this can result in more refined results.

6 Data Analysis

A common source of noise in analysis involving gene expression datasets is measurement error [18]. Gaussian graphical models are often used to inspect the relationship of different genes in varying experiments [11]. We illustrate our methodology using an Affymetrix microarray dataset containing 144 subjects of favorable histology Wilms tumors hybridized to the Affymetrix Human Genome U133A Array [8]. The data is publicly available on the GEO website, dataset GSE10320 uploaded 1/30/2009. A feature of Affymetrix data, and many other gene expression measurement platforms, is the use of multiple probes for each gene for each patient, giving replicate measurements for each patient’s gene measurement. The replicates for each patient enable an estimate of the measurement error, where we again assume the amount of contamination is independent across genes.

We follow the preprocessing steps taken in [21] and [15], which used this study in the context of measurement error in variable selection for linear models. The process begins by processing the raw data with the Bayesian Gene Expression (BGX) package [23]. BGX creates a posterior distribution for the log-scale expression level of each gene in each sample. The study recorded measurements for 22283 different genes.

To remove unnecessary computational burden, we reduced the number of genes by applying four different filters in the following order. The first filter removes expression values that do not have a corresponding Entrez gene ID in the NCBI database [16]. The second filter removes expression values with low variability by requiring at least 25% of samples to have intensities above 100 fluorescence units. The third filter removes expression values with low variability by requiring the interquartile range to be at least 0.6 on the log scale. The last filter removes expression values that have have an error to signal to noise ratio greater than 0.5, which we discuss in more depth below. After filtering, there were 273 expression values remaining for the analysis.

Now, we discuss how we estimate the measurement error of each gene. We assume that the measurement error variance is constant across patients for a given gene. We also assume that the measurement error is independent for each gene, and need not be equal for each gene. Let $\bm{\hat{\mu}}=(\hat{\mu}_{1j},\ldots,\hat{\mu}_{nj})^{T}$ denote the estimated vector of the patient’s gene expression levels for gene $j$. Further, let $\bar{\mu}=n^{-1}\sum_{j=1}^{n}\hat{\mu}_{ij}$ and $\hat{\sigma}_{j}^{2}=n^{-1}\sum_{j=1}(\hat{\mu}_{ij}-\bar{\mu}_{j})^{2}$ denote the mean and variance of each gene, respectively. For patient $i$, standardized measurements are given by $\bm{W}_{i}=(W_{i1},\ldots,W_{ip})$, calculated as $W_{ij}=\hat{\sigma}_{j}^{-1}(\hat{\mu}_{ij}-\bar{\mu}_{j})$ for each $j=1,\ldots,273$.

Let $\text{var}(\hat{\mu}_{ij})$ denote the posterior variance of the estimated distribution of patient $i$’s gene $j$. These estimates are then combined as $\hat{\sigma}_{u,j}^{2}=n^{-1}\sum_{i=1}^{n}\text{var}(\hat{\mu}_{ij})$. The measurement error covariance matrix of the standardized data $\bm{W}$ is then estimated by diagonal matrix $\hat{\bm{\Sigma}}_{u}$, where $(\hat{\bm{\Sigma}}_{u})_{j,j}=\hat{\sigma}_{u,j}^{2}/\hat{\sigma}_{j}^{2}$ for $j=1,\ldots,p$ and off-diagonal elements are 0. The fourth filter can be now formalized, where genes are removed if $\hat{\sigma}_{u,j}^{2}\geq 0.5\hat{\sigma}_{j}^{2}$; i.e. only genes with a noise-to-signal ratio less than 1 are kept for the analysis.

The original BAGUS algorithm and the corrected BAGUS algorithm were run for the remaining genes found after filtering. As with the simulations, the corrected BAGUS found fewer conditional pair-wise relationships; for this data set, the original BAGUS and corrected BAGUS found 1045 and 552 conditional pair-wise relationships, respectively. Of the 1045 naive pair-wise relationships, 42% were also found in the corrected pair-wise relationships; similarly, of the 552 corrected conditional pair-wise relationships, 80% were found in the naive model. The large percentage overlap of relationships in the corrected model with relationships in the naive model suggests that most relationships in the corrected model are true relationships. Conversely, the small percentage overlap of relationships in the naive model with those in the correct model suggests that the naive model is finding many false positive relationships. We illustrate the conditional pair-wise dependencies of the genes in Figure 2. The naive analysis is shown on the left and the corrected on the right, where the green edges signify relationships found by both procedures and purple edges signify procedure specific relationships.

7 Conclusion

We proposed a correction methodology for Gaussian graphical models when contaminated with additive measurement error. The core solution to the problem involves using the imputation-regularization algorithm to generate the true values of underlying process with a consistent estimate of the precision matrix. This provides a consistent, positive-definite estimate of the true precision matrix, which, as simulations illustrate, remove many false positive pair-wise relationships. Additionally, we show marked improvements in the AUC of the threshold matrix, indicating better separation of the underlying relationships. From a practitioner point of view, this allows for more reliable downstream analysis and further investigations to be undergone.

To our knowledge, the novel imputation-regularization algorithm has yet to be used for problems pertaining to contaminated data. This provides an avenue of future research for more a practical issue in high-dimensional problems, measurement error, which is starting to gain attention. Moreover, many practical issues still remain in the Gaussian graphical model context, such as the tuning of hyperparameters and the interpretation of the output from the Gibbs sampler-like IRO algorithm. Another potiental avenue of research to pursue is when the amount of measurement error is unknown and not assumed independent. In this case, sparsity would need to be imposed on $\bm{\Omega}_{u}$ in conjunction with $\Omega_{x}$, posing a challenging, but useful, computational procedure.

8 Appendix

8.1 Proofs

The proof for Theorem 1 in Section 3 is established here. Work done for the IRO algorithm laid the foundation for certain conditions to be met to establish consistency, see the appendix of [13]. We follow closely with their development, and prove the necessary conditions to establish consistency in our context of contaminated GGMs. These conditions include two main parts: (1) the consistency of the regularization step, specifically the BAGUS procedure in our context, and (2) some technical conditions regarding the log-likelihood $\pi(\bm{X},\bm{W})$. To that end, Assumptions 1 and 2 below ensures the consistency of the BAGUS procedure, while Assumption 3 ensures the metric entropy of the log likelihood not to grow too fast. Discussion of Assumptions 1 and 2 can be found in [6], while Assumption 3 has been commonly used in the literature of high-dimensional statistics, see the Remark 1 in the appendix of [13].

Assumption 1**.**

$\lambda_{\text{max}}(\bm{\Omega}_{x})\leq 1/k_{1}\leq\infty$, where $\lambda_{\text{max}}(\bm{\Omega}_{x})$ is the largest eigenvalue of $\bm{\Omega}_{x}$ and $k_{1}$ is a constant such that $k_{1}>0$.

For the following assumption we need to define the following values. Let the column sparsity for $\bm{\Omega}_{x}$ be denoted $b=\max_{i=1,\ldots,d}\sum_{j=1}^{d}\bm{1}(\omega_{ij}\neq 0)$. For a $m\times q$ matrix $A$ let $|||A|||_{\infty}=\max_{1\leq j\leq q}\sum_{i=1}^{d}|a_{ij}|$ be the maximum absolute row sum. Define $M_{\Sigma_{x}}=|||\bm{\Sigma}|||_{\infty}$ and $M_{\Gamma}=|||\bm{\Gamma}_{s,s}^{-1}|||_{\infty}$ where $\bm{\Gamma}=\bm{\Sigma}_{x}\otimes\bm{\Sigma}_{x}$ and $\bm{\Gamma}_{s,s}$ denotes the subset of $\bm{\Gamma}$ by indices $s=\{(i,j):\bm{\Omega}_{x}\neq 0\}$. Let $a_{1}>0$ and $a_{2}>0$ be any predefined constant value. Also, let $a_{3}$ and $k_{2}$ be defined such that $\frac{\log(d)}{n}<a_{3}<\frac{1}{4}$ and $\mathbb{E}(e^{tx^{(j)^{2}}})\leq k_{2}$ for all $|t|\leq a_{3}$ and $j=1,\ldots,d$. We define $a_{4}=a_{1}(2+a_{2}+a_{3}^{-1}k_{2}^{2})$, $a_{5}=(a_{4}+2M_{\Sigma_{x}}^{2}(a_{1}+a_{4})M_{\Gamma}+6(a_{1}+a_{4})bM_{\Gamma}^{2}M_{\Sigma}^{3}/M$. Finally, define constants $\epsilon_{0}>0$ and $\epsilon_{1}>0$, where $\epsilon_{1}$ is small.

Assumption 2**.**

For the previously defined constants, the following three statements hold:

The hyperparameters $v_{0},v_{1},\eta,$ and $\tau$ satisfy

[TABLE] 2. 2.

For the bound $||\bm{\Omega}_{x}||<B$, we have that $B$ satisfies

[TABLE] 3. 3.

For $M=\max\{2b(a_{1}+a_{4})M_{\Gamma}\max\{3M_{\Sigma},3M_{\Gamma}M_{\Sigma}^{3},\frac{2}{k_{1}^{2}}\},\frac{2a_{1}\epsilon_{0}}{k_{1}^{2}}\}$, we have $\sqrt{n}\geq M\sqrt{\log(p)}$.

Assumption 3**.**

The parameter space of $\bm{\Omega}_{x}$, or an $L_{1}$-ball containing the space of $\bm{\Omega}_{x}$, grows at a rate of $\mathcal{O}(n^{\alpha})$ for some $0\leq\alpha\leq\frac{1}{2}$.

Under these assumptions, we show that the developed procedure to correct for measurement errors satisfy the general conditions for the consistency of the IRO estimate. We state each condition and prove it to hold with our procedure.

Condition 1**.**

$\log\pi(\bm{X},\bm{W}|\bm{\Omega}_{x})$* is a continuous function of $\bm{\Omega}_{x}$ for each $\bm{x},\bm{w}\in\mathbb{R}^{d}$ and a measurable function of $(\bm{X},\bm{W})$ for each $\bm{\Omega}_{x}$.*

Proof.

We have the expansion

[TABLE]

Hence, the log posterior is continuous for symmetric positive-definite $\bm{\Omega}_{x}$ since $\bm{x}\sim N(\bm{0}_{d},\bm{\Omega}_{x}^{-1})$. The log posterior is also measurable for $(\bm{X},\bm{W})$ due to properties of the Gaussian distribution. ∎

Condition 2**.**

Three conditions for the Glivenko-Cantelli theorem to hold.

There exists a function $m_{n}(\bm{X},\bm{W})$ such that $\sup_{\bm{\Omega}_{x},\bm{X}}|\log\pi(\bm{X},\bm{W}|\bm{\Omega}_{x})|\leq m_{n}(\bm{X},\bm{W})$. 2. 2.

There exists $m_{n}^{*}(\bm{W})$, such that:

1. (a)

$0\leq\int m_{n}(\bm{X},\bm{W})\pi(\bm{X}|\bm{W},\bm{\Omega}_{x}^{(t)})d\bm{X}\leq m_{n}^{*}(\bm{W})$ for all $\bm{\Omega}_{x}^{(t)}$, 2. (b)

$\mathbb{E}[m_{n}^{*}(\bm{W})]<\infty$ , 3. (c)

$\sup_{n\in\mathbb{Z}^{+}}\mathbb{E}[m_{n}^{*}(\bm{W})\mathbb{I}(m_{n}^{*}(\bm{W})\geq\xi)]\rightarrow 0$ as $\xi\rightarrow\infty$.

Also, as $\xi\rightarrow\infty$,

[TABLE] 3. 3.

Define $\mathcal{F}_{n}=\{\int\log\pi(\bm{X},\bm{W}|\bm{\Omega}_{x})\pi(\bm{X}|\bm{W},\bm{\Omega}_{x}^{(t)})d\bm{X}\}$ and $\mathcal{G}_{n,M}=\{q\bm{1}\{m_{n}^{*}(\bm{W})\leq M\}|q\in\mathcal{F}_{n}\}$. Suppose that, for every $\epsilon$ and $M>0$, the metric entropy $\log(N(\epsilon,\mathcal{G}_{n,M},L_{1}(\mathbb{P}_{n})))=\mathcal{O}(n),$ where $\mathbb{P}_{n}$ is the emprical measure of $\bm{W}$ and $N(\epsilon,\mathcal{G}_{n,M},L_{1}(\mathbb{P}_{n}))$ is the covering number with respect to the $L_{1}(\mathbb{P})$-norm.

Proof.

We begin with part (1). Note that

[TABLE]

where $C$ contains constants not related to $(\bm{X},\bm{W},\bm{\Omega}_{x})$. Hence,

[TABLE]

where $K_{1}$ and $K_{2}$ are constants depending on upper bound $B$.

To prove part (2) note

[TABLE]

where the last equality follows from conditional independence of each $\bm{x}_{i}$. Let $\bm{\Lambda}^{(t)}=(\bm{\Omega}_{x}^{(t)}+\bm{\Omega}_{u})^{-1}$, and notice this the sum of expectations of $m(\bm{x}_{i},\bm{w}_{i})$ with respect to Gaussian random variables following $N(\bm{\Lambda}^{(t)}\bm{\Omega}_{u}\bm{w}_{i},\bm{\Lambda}^{(t)})$ for each $i=1,\ldots,n$. Now,

[TABLE]

which, since $||\bm{\Lambda}^{(t)}||\leq K_{3}$, implies

[TABLE]

Marginally $\bm{w}_{i}\sim N(\bm{0}_{d},\bm{\Sigma}_{x},\bm{\Sigma}_{u})$, and hence $m^{*}(\bm{W})$ is the sum of scaled chi-square distributions. Conditions (b) and (c) easily follow from the properties of the chi-square distribution.

To prove part (3), we make use of Remark 1 found in the Appendix of [13]. Since all elements in $\cup_{n\geq 1}\mathcal{F}_{n}$ are uniformly Lipschitz, see [7], the metric entropy can be measured on the basis of the parameter space of $\bm{\Omega}_{x}$. The functions in $\mathcal{G}_{n,M}$ are bounded and the parameter space can be contained by the $L_{1}$ ball due to the continuity of $\log\pi(\bm{X},\bm{W}|\bm{\Omega}_{x})$. By Assumption 3, then $\log(N(\epsilon,\mathcal{G}_{n,M},L_{1}(\mathbb{P}_{n})))=\mathcal{O}(n^{2\alpha}\log(d)).$ ∎

Condition 3**.**

Define $Z_{t,i}=\log\pi(\bm{x}_{i},\bm{w}_{i}|\bm{\Omega}_{x})-\int\log\pi(\bm{x}_{i},\bm{w}_{i}|\bm{\Omega}_{x})\pi(\bm{X}|\bm{w}_{i},\bm{\Omega}_{x}^{(t)})$. $Z_{t,i}$ are subexponential random variables.

Proof.

First, we note that

[TABLE]

where $C_{1}$ is a constant free of $\bm{X}$. Also note $\log\pi(\bm{w}_{i},\bm{X}|\bm{\Omega}_{x})=\log\pi(\bm{w}_{i},\bm{x}_{i}|\bm{\Omega}_{x})+\log\pi(\bm{X}_{-i}|\bm{\Omega}_{x})$. The integral can then be shown to be

[TABLE]

The value of $A$ is the expectation of $\log\pi(\bm{w}_{i},\bm{x}_{i}|\bm{\Omega}_{x})$ with respect to the full conditional of $X$ at iteration $t$, $\bm{x}_{i}|\bm{w}_{i},\bm{\Omega}_{x}^{(t)}\sim N_{d}(\bm{\Lambda}^{-1,(t)}\bm{\Omega_{u}}\bm{w}_{i},\bm{\Lambda}^{-1,(t)})$ where $\bm{\Lambda}^{(t)}=(\bm{\Omega}_{x}^{(t)}+\bm{\Omega_{u}})$. This expectation is composed of two parts,

[TABLE]

and

[TABLE]

Hence, $Z_{t,i}$ is

[TABLE]

where $C=C_{1}+C_{2}$ is free of $\bm{x}_{i}$ and $\bm{w}_{i}$, which is the sum of scaled chi-squared distributions and thus is subexponential. ∎

Condition 4**.**

For $t=1,\ldots,T$, $Q(\bm{\Omega}_{x}|\bm{\Omega}_{x}^{(t)})$ has a unique maximum at $\tilde{\bm{\Omega}}_{x}^{(t)}$; for any $\epsilon>0$, $\sup_{\bm{\Omega}_{x}\backslash B_{t}(\epsilon)}Q(\bm{\Omega}_{x}|\bm{\Omega}_{x}^{(t)})$ exists, where $B_{t}(\epsilon)=\{\bm{\Omega}_{x}:|\bm{\Omega}_{x}-\tilde{\bm{\Omega}}_{x}^{(t)}|<\epsilon\}$.

Proof.

As noted in [13], this is satisfied if $\bm{\Omega}_{x}$ is restricted to a compact set. So, since BAGUS is strictly convex when restricted by the condition that $||\bm{\Omega}_{x}||\leq B$, then the condition is satisfied. ∎

Condition 5**.**

The penalty function is non-negative, ensures the existence of $\bm{\Omega}_{x}^{(t+1)}$ for $t=2,\ldots,T$, and converges to 0 uniformly as $n\rightarrow\infty$.

Proof.

BAGUS is a non-negative penalty that exists for any $\bm{X}$, and, due to the adaptive nature of the penalty, converges to 0 as $n\rightarrow\infty$. To see the penalty converges to 0, note Assumption 2.1a implies

[TABLE]

as $n\rightarrow\infty$, which, with a similar argument for $v_{0}$, results in the penalty being 0 as $n\rightarrow\infty$. ∎

8.2 Computing BAGUS with the EM-Algorithm

Here we review the optimization of the uncontaminated objective distribution. The direct optimization of $L^{UC}$ in (5) is not easy due to the sum inside the logarithm. [6] use the EM-algorithm to get around this issue by introducing the latent factors $r_{ij}$ from section 3.3. This allows an E-step similar to the spike-and-slab Lasso and an M-step similar to the Graphical Lasso. In this section, if not specified, $\bm{\Sigma}$ and $\bm{\Omega}$ refer to $\bm{x}$’s covariance and precision matrix, respectively.

The optimization seeks to find the MAP of the posterior proportional to

[TABLE]

where the latent indicator $r_{ij}$, as defined in Section 3.3, is incorporated into the off-diagonal elements in the prior specification. The E-step takes the conditional expectation of $r_{ij}$ in the proportional posterior. Each $r_{ij}$ is conditionally Bernoulli with probability

[TABLE]

allowing for easy calculation of the conditional expectation. Then, the desired $Q$ function to maximize in the M-step is given by

[TABLE]

where the expectation is taken element wise for $\bm{R}$ by assumed independence of inclusion.

The M-step optimizes each column of $Q$ seperately with coordinate descent. The last column’s update is now explained, with the other columns following in the same pattern. Partition the covariance matrix as

[TABLE]

and let similar partitions be available for $\bm{\Omega}_{x},\bm{P},\bm{R},\text{and }\bm{S}$. Also note that

[TABLE]

where $c=\omega_{22}-\bm{\omega}_{12}^{T}\bm{\Omega}_{11}^{-1}\bm{\omega}_{12}.$ The update for the last column of $\bm{\Sigma}_{x}$ is the solution from setting subgradient of $Q$ with respect to $\left[\bm{\sigma}_{12}\ \bm{\sigma}_{22}\right]^{T}$ to 0. The update for $\sigma_{22}$ can is easily attained from the setting the subgradient of $\omega_{22}$ to 0,

[TABLE]

We note that each column update requires the matrix $\bm{\Omega}_{11}^{-1}$. This can be computed as $\bm{\Sigma}_{11}-\bm{\sigma}_{12}\bm{\sigma}_{12}^{T}/\sigma_{22}$.