Updating of the Gaussian graphical model through targeted penalized estimation

Abstract

Updating of the Gaussian graphical model via shrinkage estimation is studied. This shrinkage is towards a nonzero parameter value representing prior quantitative information. Once new data become available, the previously estimated parameter needs updating. Shrinkage provides the means to this end, using the latter as a shrinkage target to acquire an updated estimate. The process of iteratively updating the Gaussian graphical model in shrinkage fashion is shown to yield an improved fit and an asymptotically unbiased and consistent estimator. The workings of updating via shrinkage are elucidated by linking it to Bayesian updating and through the inheritance by the update of eigen-properties of the previous estimate. The effect of shrinkage on the moments and loss of the estimator are pointed out. Practical issues possibly hampering updating are identified and solutions outlined. The presented updating procedure is illustrated through the reconstruction of a gene–gene interaction network using transcriptomic data.

Introduction

The reconstruction of networks from data is a key challenge of our time, as can be witnessed from the number of citations that recently published monographs on graphical models by [38] and [20] have accumulated. Many methods that serve this end have been proposed since, cf., the recently published review by [10]. Virtually all these methods reconstruct such networks de novo. That is, they do not make use of information on these networks already available. For molecular networks such information is condensed not only in repositories like KEGG and STRING [27], [35], but also in the form of omics data from the model system under comparable conditions available through portals like GEO and TCGA [11], [39]. Information from aforementioned sources may enhance the molecular network reconstruction from the data at hand. A thus reconstructed network qualifies itself as novel information on the network. In turn it needs updating when new (experimental) information on the network becomes available in the future. Such iterative reconstruction will progress the knowledge on the network. In this work we explore how this may be achieved.

A network is defined as an undirected graph $\mathcal{G}$ comprising the pair $(\mathcal{V},\mathcal{E})$. $\mathcal{V}$ is the index set of the nodes – representing the (molecular) entities – of the network. $\mathcal{E}=\mathcal{V}\times \mathcal{V}$ is the edge set, consisting of node pairs that are connected. The graphs under study are undirected, which implies that if node pair $(v_1,v_2)\in \mathcal{E}$ automatically $(v_2,v_1)\in \mathcal{E}$ is included too. An edge is operationalized as a pairwise conditional dependency among the random variables represented by the nodes connected by the edge, given all other random variables represented by the remaining nodes in the network. The absence of an edge is vice versa defined as a conditional independency. The conditional independence graph is to be learned from (say) transcriptomic data on the genes that constitute the network. To this end a Gaussian graphical model (GGM) is assumed: $\mathcal{N}({\mathbf{0}}_p,{\mathit{\Omega }}^{-1})$ with $\mathit{\Omega }$ the inverse covariance matrix $\mathit{\Omega }={\mathit{\Sigma }}^{-1}$. Zeroes and nonzeroes in the inverse covariance matrix $\mathit{\Omega }$ correspond to conditional in- and dependencies, respectively, among the corresponding variate pairs. It thus suffices to estimate the inverse covariance matrix from the data by maximization of the log-likelihood, $\propto log\left(\left|\mathit{\Omega }\right|\right)-\text{tr}\left(\mathbf{S}\mathit{\Omega }\right)$ with sample covariance matrix $\mathbf{S}$ obtained from a sample of size $n$, yielding the maximum likelihood (ML) estimator ${\widehat{\mathit{\Omega }}}_{ML}={\mathbf{S}}^{-1}$. Rests to infer its support.

The estimation of a Gaussian graphical model is hampered by the collinearity of omics data. Then, the maximum likelihood estimator of the inverse covariance matrix is not (or ill-) defined as $\mathbf{S}$ is (near-)singular. This is overcome by the augmentation of the log-likelihood loss with a penalty. Here we concentrate on the nonzero centered ridge penalty, $\lambda {\Vert \mathit{\Omega }-{\mathbf{T}}_r\Vert }_F^2$, the sum of the square of the elements of $\mathit{\Omega }-{\mathbf{T}}_r$. With the subscript $r$ referring to ridge, this ${\mathbf{T}}_r$ is a user-specified target matrix towards which the precision estimate is shrunken for large values of $\lambda$. Maximization of this ridge penalized log-likelihood yields the ridge precision estimator (see [42]):

$${\widehat{\mathit{\Omega }}}_r\left(\lambda \right)={\{\frac{1}{2}(\mathbf{S}-\lambda {\mathbf{T}}_r)+{[\lambda {\mathbf{I}}_{pp}+\frac{1}{4}{(\mathbf{S}-\lambda {\mathbf{T}}_r)}^2]}^{1∕2}\}}^{-1}.$$

[42] shows that this estimator satisfies the following properties: (i) ${\widehat{\mathit{\Omega }}}_r\left(\lambda \right)\succ 0$, (ii) ${lim}_{\lambda \downarrow 0}{\widehat{\mathit{\Omega }}}_r\left(\lambda \right)={\mathbf{S}}^{-1}$ (should it exist), (iii) ${lim}_{\lambda \to \infty }{\widehat{\mathit{\Omega }}}_r\left(\lambda \right)={\mathbf{T}}_r$, and (iv) consistency. Furthermore, for a suitable choice of $\lambda$ it – using a diagonal, equivariant target ${\mathbf{T}}_r=\alpha {\mathbf{I}}_{pp}$ with $\alpha \ge 0$ – also outperforms the ML precision estimator in terms of the mean squared error (MSE) [40]. In that MSE sense, even the inverse of the ridge precision estimator, ${\widehat{\mathit{\Sigma }}}_r\left(\lambda \right)={\left[{\widehat{\mathit{\Omega }}}_r\left(\lambda \right)\right]}^{-1}$, outperforms the ML covariance estimator. Although the focus of the present work is the ridge precision estimator, an alternative to the ridge precision estimator would be the inverse of the Ledoit–Wolf shrunken covariance matrix [21], ${\widehat{\mathit{\Omega }}}_{\ell w}\left(\nu \right)={\left[{\widehat{\mathit{\Sigma }}}_{\ell w}\left(\nu \right)\right]}^{-1}$ with ${\widehat{\mathit{\Sigma }}}_{\ell w}\left(\nu \right)=(1-\nu )\mathbf{S}+\nu {\mathbf{T}}_{\ell w}$ for $\nu \in (0,1)$. Note that ${\mathbf{T}}_{\ell w}$ is a covariance target, while ${\mathbf{T}}_r$ is a precision target. Clearly, this Ledoit–Wolf shrinkage covariance estimator shrinks to $\mathbf{S}$ and ${\mathbf{T}}_{\ell w}$ as $\nu \downarrow 0$ and $\nu \uparrow 1$, respectively. Moreover, ${\widehat{\mathit{\Sigma }}}_{\ell w}\left(\nu \right)$ too may – for a suitably chosen shrinkage parameter and a diagonal, equivariant target (although it has also been proven for other target choices, [21], [31]) – outperform the ML covariance estimator in the mean squared error sense.

The target provides means to include a quantitative suggestion of the precision. The number of choices for ${\mathbf{T}}_r$ and ${\mathbf{T}}_{\ell w}$ is in principle unlimited. [16], [31] suggest several low-dimensionally parameterized targets. Some archetypical examples (modified for the ridge precision estimator) are:

• The equivariant and rotational invariant target ${\mathbf{T}}_r=\alpha {\mathbf{I}}_{pp}$ with $\alpha >0$.

• A diagonal ${\mathbf{T}}_r$ with ${\left({\mathbf{T}}_r\right)}_{jj}=1∕{\mathbf{S}}_{jj}$, such that only off-diagonal elements are shrunken.

• A simple model of the covariance and variance by ${\mathbf{T}}_r={[{\sigma }^2(1-\rho ){\mathbf{I}}_{pp}+{\sigma }^2\rho {\mathbf{1}}_{pp}]}^{-1}$, with ${\sigma }^2>0$ and $\rho \in [-1,1]$.

[16], [31] learn the free parameters of ${\mathbf{T}}_r$ from the data at hand through, e.g. the minimization of the mean squared error (under assumptions). As an alternative to these low-dimensionally parameterized and data-derived targets, one may also consider a fully unstructured target obtained through other means than the data at hand, e.g.; ${\mathbf{T}}_r={\mathit{\Omega }}_{\text{known}}$ taken from a related and well-studied model system. Of course, the inverse of these ${\mathbf{T}}_r$ may serve as covariance targets for the Ledoit–Wolf shrinkage covariance estimator ${\widehat{\mathit{\Sigma }}}_{\ell w}\left(\nu \right)$.

Here we describe the problem of updating a Gaussian graphical model after the arrival of novel information. The above introduced shrinkage estimators are employed to accommodate the prior information on the model (in the form of a target that represents the current knowledge/estimate on/of the model at that point in time and formed from the current and preceding data sets) to arrive at an updated estimate from the novel data. Properties of these updating estimators are investigated as well as intuition into their workings is provided. Resolutions for practical complications when updating, batch effects and the simultaneous arrival of multiple sources of novel information, are presented. The paper concludes with an example in which the precision matrix and the therefrom derived gene–gene interaction network of ten pathways are reconstructed from transcriptomic data of ER+ breast cancer samples of a series of five subsequently published data sets. Throughout, proofs are deferred to the appendix, while the supplementary material contains more elaborate versions of the proofs as well as additional results.

Estimation of precision matrices from multiple data sets has been considered in the literature. [6], [8], [29] deal with the simultaneous estimation of multiple precision matrices using so-called ‘fused’ ridge/lasso penalties, while e.g. [25], [28] adopt a Bayesian approach. The methods of [15], [43] have a similar thrust but these works are based on parameter communalities between the different precision matrices. On a different note but closely related is the work of [25], [26], [44], where the focus is on the reconstruction of multiple directed acyclic graphs in order to facilitate causal inference. Implicitly, these methods assume that the data used for the estimation of each precision matrix are on par, in the sense that no time order is present. This absence of an order motivates the use of a fused penalty that encourages to borrow information across data sets (should the data give rise to do so). In particular, this shrinks the different precision estimates towards each other, but not towards a target matrix that is formed from one of the other data sets. Moreover, this ‘across data set’-calibration is done once, instead of repetitively in the envisioned updating context considered here.