ℓ₁-Penalized censored Gaussian graphical model

Summary

1. Introduction

An important aim in genomics is to understand interactions among genes, characterized by the regulation and synthesis of proteins under internal and external signals. These relationships can be represented by a genetic network, i.e. a graph where nodes represent genes and edges describe the interactions among them. Genetic networks can be used to gain new insights into the activity of biological pathways and to deduce unknown functions of genes from their dependence on other genes.

Gaussian graphical models (Lauritzen, 1996) have been widely used for reconstructing a genetic network from expression data. The reason of their widespread use relies on the statistical properties of the multivariate Gaussian distribution, which allow the topological structure of a network to be related to the non-zero elements of the concentration matrix, i.e. the inverse of the covariance matrix. Thus, the problem of network inference can be recast as the problem of estimating a concentration matrix. The graphical lasso (Yuan and Lin, 2007) is a popular method for estimating a sparse concentration matrix, based on the idea of adding an ${\ell }_1$-penalty to the likelihood function of the multivariate Gaussian distribution. Nowadays, this estimator is widely used in applied research (e.g. Menéndez and others, 2010; Vinciotti and others, 2016) and widely studied in the computational as well as theoretical (e.g. Bickel and Levina, 2008; Friedman and others, 2008; Witten and others, 2011) literature. The interested reader is referred to Augugliaro and others (2016) for an extensive review.

Despite the widespread literature on the graphical lasso estimator, there is a great number of fields in applied research where modern measurement technologies make the use of this graphical model theoretically unfounded, even when the assumption of a multivariate Gaussian distribution is satisfied. A first example of this is Reverse Transcription quantitative Polymerase Chain Reaction (RT-qPCR), a popular technology for gene expression profiling (Derveaux and others, 2010). This technique is used to measure the expression of a set of target genes in a given sample through repeated cycles of sequence-specific DNA amplification followed by expression measurements. The cycle at which the observed expression first exceeds a fixed threshold is commonly called the cycle-threshold (McCall and others, 2014). If a target is not expressed, the threshold is not reached after the maximum number of cycles (limit of detection) and the corresponding cycle-threshold is undetermined. For this reason, the resulting data is naturally right-censored (McCall and others, 2014; Pipelers and others, 2017). Another example is given by the flow cytometer, which is an essential tool in the diagnosis of diseases such as acute leukemias and malignant lymphomas (Brown and Wittwer, 2000). A flow cytometer measures a limited range of signal strength and records each marker value within a fixed range, such as between 0 and 1023. If a measurement falls outside this range, then the value is replaced by the nearest legitimate value; that is, a value smaller than 0 is censored to 0 and a value larger than 1023 is censored to 1023. A direct application of the graphical lasso for network inference from data such as these is theoretically unfunded since it does not consider the effects of the censoring mechanism on the estimator of the concentration matrix.

In this article, we propose an extension of the graphical lasso estimator that takes into account the censoring mechanism of the data explicitly. Our work can be related to Städler and Bühlmann (2012), who propose an ${\ell }_1$-penalized estimator of the inverse covariance matrix of a multivariate Gaussian model based on the assumption that the data are missing at random. As we shall see in the following of this article, failure to take into account the censoring mechanism causes a poor behavior of this approach when compared with our proposal. Our proposal can also be related to the work of Perkins and others (2013), Hoffman and Johnson (2015), and Pesonen and others (2015), who provide a maximum likelihood estimator of the covariance matrix under left-censoring. However, these works do not address the estimation of the precision matrix under a sparsity assumption and for this reason they are applicable only when the sample size is larger than the number of nodes in the network.

The remaining part of this article is structured as follows. In Section 2, we extend the notion of Gaussian graphical model to censored data and in Section 3, we propose the extension of the graphical lasso estimator to a censored graphical lasso estimator and two Expectation-Maximization (EM) algorithms for inference of parameters in the censored Gaussian graphical model. Section 4 is devoted to the evaluation of the behavior of the proposed algorithms and the comparison of the proposed estimator with existing competitors. Finally, in Section 5, we study a real dataset and in Section 6, we draw some conclusions.

2. The censored Gaussian graphical model

In order to describe the technical details of the proposed method, let $\mathit{X}=(X_1,\dots ,X_p{)}^{\mathrm{\top }}$ be a $p$-dimensional random vector. Graphical models allow to represent the set of conditional independencies among these random variables by a graph $\mathcal{G}=\{\mathcal{V},\mathcal{E}\}$⁠, where $\mathcal{V}$ is the set of nodes associated to $\mathit{X}$ and $\mathcal{E}\subseteq \mathcal{V}\times \mathcal{V}$ is the set of ordered pairs, called edges, representing the conditional dependencies among the $p$ random variables (Lauritzen, 1996).

As shown in Lauritzen (1996), the pattern of non-zero elements in $\mathrm{\Theta }$ defines the corresponding graph, namely the undirected edge $(h,k)$ is an element of the edge set $\mathcal{E}$ of the corresponding conditional independence graph if and only if ${\theta }_{hk}\ne 0$⁠.

A closer look at Definition 1 reveals that the proposed notion of censored Gaussian graphical model is characterized by a high degree of generality since it covers also the special cases of the classical Gaussian graphical model (⁠$l_h=-\mathrm{\infty }$⁠, $u_h=+\mathrm{\infty }$ for any $h$⁠), and the cases in which there is only left-censored data (⁠$u_h=+\mathrm{\infty }$ for any $h$⁠) or right-censored data (⁠$l_h=-\mathrm{\infty }$ for any $h$⁠).

3. ${\ell }_1$-Penalized estimator for censored Gaussian graphical model

3.1. The censored graphical lasso estimator

Consider a sample of $n$ independent observations drawn from the censored Gaussian graphical model $\{\mathit{X},R(\mathit{X};\mathit{l},\mathit{u}),\phi ({\mathit{x}}_o,\mathit{r};\mathit{\mu },\mathrm{\Theta }),\mathcal{G}\}$⁠. For ease of exposition, we shall assume that $\mathit{l}$ and $\mathit{u}$ are fixed across the $n$ observations, but the extension to the cases where the censoring vectors are specific to each observation is straightforward and does not require a specific treatment. If these values are unknown we suggest to use the smallest and largest observed value for each variable as possible estimates for $l_h$ and $u_h$⁠, respectively.

Like in the standard graphical lasso estimator, the nonnegative tuning parameter $\rho$ is used to control the amount of sparsity in the estimated concentration matrix ${\hat{\mathrm{\Theta }}}^{\rho }=({\hat{\theta }}_{hk}^{\rho })$ and, consequently, in the corresponding estimated graph ${\hat{\mathcal{G}}}^{\rho }=\{\mathcal{V},{\hat{\mathcal{E}}}^{\rho }\}$⁠, where ${\hat{\mathcal{E}}}^{\rho }=\{(h,k):{\hat{\theta }}_{hk}^{\rho }\ne 0\}$⁠. When $\rho$ is large enough, some ${\hat{\theta }}_{hk}^{\rho }$ are shrunken to zero resulting in the removal of the corresponding link in ${\hat{\mathcal{G}}}^{\rho }$⁠; on the other hand, when $\rho$ is equal to zero and the sample size is large enough the estimator ${\hat{\mathrm{\Theta }}}^{\rho }$ coincides with the maximum likelihood estimator of the concentration matrix, which implies a fully connected estimated graph.

3.2. Fitting the censored graphical lasso model

1. 1.

   Let $\{{\hat{\mathit{\mu }}}_{ini}^{\rho },{\hat{\mathrm{\Theta }}}_{ini}^{\rho }\}$ be initial estimates;

2. 2.

   E-step: Compute $\overline{\mathit{x}}({\hat{\mathit{\mu }}}_{ini}^{\rho },{\hat{\mathrm{\Theta }}}_{ini}^{\rho })$ and $S({\hat{\mathit{\mu }}}_{ini}^{\rho },{\hat{\mathrm{\Theta }}}_{ini}^{\rho })$⁠;

3. 3.

   M-step: Let ${\hat{\mathit{\mu }}}^{\rho }=\overline{\mathit{x}}({\hat{\mathit{\mu }}}_{ini}^{\rho },{\hat{\mathrm{\Theta }}}_{ini}^{\rho })$ and solve the problem (3.5) using $S({\hat{\mathit{\mu }}}_{ini}^{\rho },{\hat{\mathrm{\Theta }}}_{ini}^{\rho })$⁠;

4. 4.

   If a convergence criterion is met then return $\{{\hat{\mathit{\mu }}}^{\rho },{\hat{\mathrm{\Theta }}}^{\rho }\}$ else let ${\hat{\mathit{\mu }}}_{ini}^{\rho }={\hat{\mathit{\mu }}}^{\rho }$ and ${\hat{\mathrm{\Theta }}}_{ini}^{\rho }={\hat{\mathrm{\Theta }}}^{\rho }$⁠;

5. 5.

   Repeat steps 2–4.

As shown by Guo and others (2015), the approach works well when the real concentration matrix is sparse or the tuning parameter is sufficiently large. The main advantage coming from the approximation (3.6) is that now, in order to evaluate the quantities $\overline{\mathit{x}}({\hat{\mathit{\mu }}}_{ini}^{\rho },{\hat{\mathrm{\Theta }}}_{ini}^{\rho })$ and $S({\hat{\mathit{\mu }}}_{ini}^{\rho },{\hat{\mathrm{\Theta }}}_{ini}^{\rho })$⁠, we only need to compute $E_{c_i\mid o_i}(X_{ih}\mid {\mathit{X}}_{i{c}_i}\in D_{c_i})$ and $E_{c_i\mid o_i}(X_{ih}^2\mid {\mathit{X}}_{i{c}_i}\in D_{c_i})$⁠, for $h=1,\dots ,p$⁠. These quantities can be computed by using exact formulas (Johnson and others, 1994) and require only the evaluation of the cumulative distribution function of the univariate Gaussian distribution. In the following of this article, we denote by $\overline{S}({\hat{\mathit{\mu }}}_{ini}^{\rho },{\hat{\mathrm{\Theta }}}_{ini}^{\rho })$ the approximate estimate of the empirical covariance matrix resulting from this approach. By an extensive simulation study, in Section 4.1, we shall study the behavior of the proposed algorithm based on the usage of the matrix $\overline{S}({\hat{\mathit{\mu }}}_{ini}^{\rho },{\hat{\mathrm{\Theta }}}_{ini}^{\rho })$ in steps 2 and 3.

As with graphical lasso, our proposed method requires a sequence of $\rho$-values, which should be suitably defined so as to reduce the computational cost needed to compute the entire path of the estimated parameters. Theorem 2 gives the exact formula for the derivation of the largest $\rho$-value, denoted by ${\rho }_{max}$⁠, and the corresponding cglasso estimator.

A proof of Theorem 2 is reported in the supplementary material available at Biostatistics online. Theorem 2 shows how the maximum value of $\rho$ can be efficiently computed using $p$ optimization problems which do not involve the tuning parameter, and then calculating $S(\hat{\mathit{\mu }},\hat{\mathrm{\Theta }})$⁠. Once ${\rho }_{max}$ is calculated, the entire path of the cglasso estimates is calculated as follow:

1. 1.

   Compute ${\rho }_{max}$ as specified in Theorem 2 and let ${\rho }_{min}$ be the smallest $\rho$-value;

2. 2.

   Compute a decreasing sequence $\{{\rho }_{(k)}{\}}_{k=1}^K$ of distinct $\rho$-values starting from ${\rho }_{max}$ to ${\rho }_{min}$⁠;

3. 3.

   for $k=2$ to $K$ do

4. 4.

     Let ${\hat{\mathit{\mu }}}_{ini}^{{\rho }_{(k)}}={\hat{\mathit{\mu }}}^{{\rho }_{(k-1)}}$ and ${\hat{\mathrm{\Theta }}}_{ini}^{{\rho }_{(k)}}={\hat{\mathrm{\Theta }}}^{{\rho }_{(k-1)}}$

5. 5.

     Use the EM algorithm to compute $\{{\hat{\mathit{\mu }}}^{{\rho }_{(k)}},{\hat{\mathrm{\Theta }}}^{{\rho }_{(k)}}\}$ with $\{{\hat{\mathit{\mu }}}_{ini}^{{\rho }_{(k)}},{\hat{\mathrm{\Theta }}}_{ini}^{{\rho }_{(k)}}\}$ as starting values;

6. 6.

   end for

The previous description shows that the entire path can be computed using the estimates obtained for a given $\rho$-value as warm starts for fitting the next cglasso model. This strategy is commonly used also in other efficient lasso algorithms and R packages (Friedman and others, 2010, Friedman and others, 2018). In our model, this strategy turns out to be remarkably efficient since using a sufficiently fine sequence of $\rho$-values, the starting values defined in Step 4 will be sufficiently close to the estimates computed in Step 5, thus increasing the speed of convergence of the resulting algorithm.

3.3. Tuning parameter selection

4. Simulation studies

4.1. Evaluating the behavior of the approximate EM algorithm

In a first simulation, we evaluate the effects of the approximation (3.6) on the accuracy of the estimators obtained with the EM algorithm. As in the real dataset studied in Section 5, we consider right-censored data and we set the right-censoring value $u_h$ equal to 40, for any $h=1,\dots ,p$⁠. For this simulation, where we plan to use the full EM algorithm, we set the number of variables $p$ to 10 and the sample size $n$ to 100. To simulate a sample from a sparse right-censored model, we use the following procedure. First, using the method implemented in the R package huge (Zhao and others, 2015), we simulate a sparse concentration matrix with a random structure, where we set the probability that a ${\theta }_{hk}$ is different from zero to 0.1. Then, the elements of the parameter $\mathit{\mu }$ are selected to obtain a fixed probability of right censoring for any given random variable. More specifically, we randomly draw a subset $\mathcal{D}$ from $\mathcal{I}=\{1,\dots ,p\}$ and for each $h\in \mathcal{D}$ the value of the parameter ${\mu }_h$ is such that $\text{Prob}\{R(X_h;-\mathrm{\infty },40)=+1\}=0.25$⁠, while for each $h\notin \mathcal{D}$ the parameter ${\mu }_h$ is such that the probability of right censoring is approximately equal to ${10}^{-11}$⁠. The cardinality of the set $\mathcal{D}$⁠, denoted by $|\mathcal{D}|$⁠, is used in the study as a tool to analyze the effects of the number of censored variables on the behavior of the proposed EM algorithm and its approximated version. Finally, we draw a sample from a multivariate Gaussian distribution with parameters given in the previous steps and we treat each value greater than 40 as a missing (censored) value. We simulated 100 datasets from this model and for each simulation we compute a path of cglasso estimates using the proposed EM algorithm and a path using the approximated EM algorithm. In the following of this section, the estimates belonging to the two paths are denoted by $\{{\hat{\mathit{\mu }}}_e^{\rho };{\hat{\mathrm{\Theta }}}_e^{\rho }\}$ and $\{{\hat{\mathit{\mu }}}_a^{\rho };{\hat{\mathrm{\Theta }}}_a^{\rho }\}$⁠, respectively. For each path, the largest value of the tuning parameter was computed using the results given in Theorem 2 while the smallest value was set equal to $1\times {10}^{-3}$⁠.

The first part of the Table 1 reports the average CPU times for computing the path. As expected, the computational time needed to compute a path is always an increasing function of the number of the censored variables but, when we use the matrix $S({\hat{\mathit{\mu }}}_{ini}^{\rho },{\hat{\mathrm{\Theta }}}_{ini}^{\rho })$ in the E-step, the table reveals that the form of the relationship between the CPU time and $|\mathcal{D}|$ is almost exponential implying that the computation of the cglasso estimator is infeasible also for relatively small datasets. In contrast to this, when we use the matrix $\overline{S}({\hat{\mathit{\mu }}}_{ini}^{\rho },{\hat{\mathrm{\Theta }}}_{ini}^{\rho })$ to approximate the current estimate of the empirical covariance matrix, the table shows an almost linear dependence of the CPU time, which implies a significant reduction of the computational complexity compared with the full EM algorithm. Although the results showed in Table 1 strongly suggest the use of the approximated EM algorithm to compute the cglasso estimator, they do not provide information about the difference between the two estimates. For this reason, we also computed the largest Euclidean distance between ${\hat{\mathit{\mu }}}_e^{\rho }$ and ${\hat{\mathit{\mu }}}_a^{\rho }$⁠, denoted as $\underset{\rho }{max}\Vert \mathrm{\Delta }{\hat{\mathit{\mu }}}^{\rho }{\Vert }^2$⁠, and the largest Frobenius distance between ${\hat{\mathrm{\Theta }}}_e^{\rho }$ and ${\hat{\mathrm{\Theta }}}_a^{\rho }$⁠, denoted by $\underset{\rho }{max}\Vert \mathrm{\Delta }{\hat{\mathrm{\Theta }}}^{\rho }{\Vert }_F^2$⁠. The average and standard deviations of these values are reported in the second part of the Table 1. All the results clearly show that the two estimators are sufficiently close to each other, and point to the use of the approximated EM algorithm for the derivation of the cglasso estimator.

4.2. Comparison of methods on data simulated from a censored Gaussian graphical model

In a second simulation study, we compare our proposed estimator with MissGlasso (Städler and Bühlmann, 2012), which performs ${\ell }_1-$penalized estimation under the assumption that the censored data are missing at random, and with the glasso estimator (Friedman and others, 2008), where the empirical covariance matrix is calculated by imputing the missing values with the limit of detection. These estimators are evaluated in terms of both recovering the structure of the true graph and the mean squared error. We use a similar approach to the previous simulation for generating right censored data. In particular, we set the right censoring value to 40 for any variable and the sample size $n$ to 100. We generate a sparse concentration matrix with random structure and set the probability of observing a link between two nodes to $k/p$⁠, where $p$ is the number of variables and $k$ is used to control the amount of sparsity in $\mathrm{\Theta }$⁠. Finally, we set the mean $\mathit{\mu }$ in such a way that ${\mu }_h=40$ for the $H$ censored variables, i.e. $\text{Prob}\{R(X_h;-\mathrm{\infty },40)=+1\}=0.50$⁠, while for the remaining variables ${\mu }_h$ is sampled from a uniform distribution on the interval $[10;35]$⁠. At this point, we simulate a sample from the latent $p$-variate Gaussian distribution and treat all values greater than 40 as missing. The quantities $k$⁠, $p$⁠, and $H$ are used to specify the different models used to analyze the behavior of the considered estimators. In particular, we consider the following cases:

• Model 1: $k=3$⁠, $p=50$ and $H\in \{25,35\}$⁠. This setting is used to evaluate the effects of the number of censored variables on the behavior of the proposed estimators when $n>p$⁠.

• Model 2: $k\in \{1,5\}$⁠, $p=50$ and $H=30$⁠. This setting is used to evaluate the effects of the sparsity of the matrix $\mathrm{\Theta }$ on the considered estimators when $n>p$⁠.

• Model 3: $k=3$⁠, $p=200$ and $H=100$⁠. This setting is used to evaluate the impact of the high dimensionality on the estimators (⁠$p\gg n$⁠).

Fig. 1.

Open in new tabDownload slide

Results of the comparative simulation study based on Model 1 with $H=25$⁠. (a) The average precision-recall curves are shown; (b) The box-plots of $\underset{\rho }{min}\text{MSE}({\hat{\mathrm{\Theta }}}^{\rho })$ for the considered estimators is shown.

The curves show how cglasso gives a better estimate of the concentration matrix both in terms of precision and recall, for any given value of the tuning parameter. Figure 1b shows the distributions of the quantity $\underset{\rho }{min}\text{MSE}({\hat{\mathrm{\Theta }}}^{\rho })$⁠, which gives the minimum value of the mean squared error attained along the path of solutions. These box-plots emphasize that, not only the cglasso has a mean squared error much smaller than glasso and MissGlasso, but also that it is much more stable than its competitors. The same behavior was also observed in the other models used in our numerical study, and can be found in the Supplementary Material available at Biostatistics online. Table 2 reports the summary statistics from all the simulations. In addition to the quantity described above (second meta-column), the first meta-column shows the mean squared error in the estimation of the mean (denoted by $\underset{\rho }{min}\text{MSE}({\hat{\mathit{\mu }}}^{\rho })$⁠) and the third meta-column reports the Area Under the precision-recall Curve (AUC). Note that the considered measures allow us to study the behavior of the estimators along the entire path. The distribution of the minimum value of the mean squared errors shows that, not only our estimator is able to recover the structure of the graph but also outperforms the competitors in terms of both estimation of $\mathit{\mu }$ and $\mathrm{\Theta }$⁠. We did not report $\underset{\rho }{min}\text{MSE}({\hat{\mathit{\mu }}}^{\rho })$ for glasso since this method does not allow to estimate the parameter $\mathit{\mu }$⁠. The results on the AUC suggest that cglasso can be used as an efficient tool for recovering the structure of the true concentration matrix of a Gaussian graphical model from censored data.

4.3. Testing the robustness of the method on more realistic biological data

In a third simulation study, we test the robustness of the method on more realistic biological data. In particular, we consider expression data from Wille and others (2004) on the extensively studied Arabidopsis thaliana biological system. The study reports data from $n=118$ experiments on $p=39$ genes, whose regulatory network is of interest. Although the data are fully observed, we test our method on a dataset where observations are made artificially right censored, similarly to Städler and Bühlmann (2012). In particular, we produce three datasets at different levels of right censoring, by recording as missing the 10%, 20%, and 30% of the highest values, respectively. Then, we compare our method, cglasso, with MissGlasso (Städler and Bühlmann, 2012) and with three additional methods, which impute missing values first and then infer the network using graphical lasso on the imputed values. As in Städler and Bühlmann (2012), we consider k-nearest neighbor imputation, which we denote with missknn, and imputation using random forests, which we denote with MissForest. These methods are based on a missing at random assumption. We further consider a method which relies on multivariate Gaussian data under a censoring mechanism which we denote with imputeLCMD (Lazar, 2015).

Table 3 shows the Euclidean norm between the true observed data and the imputed one for the five methods considered and across the three different levels of censoring. The results show how the performance of all methods decreases the higher the level of censoring and how cglasso is overall superior to the other methods across all comparisons, followed closely in some cases by imputeLCMD. The areas under the precision-recall curves of the networks inferred by the five methods compared with the graphical lasso network selected by ${\overline{\text{BIC}}}_{0.5}({\hat{\mathcal{E}}}^{\rho })$ from the fully observed data show how the proposed cglasso method leads to a significant gain in network recovery over MissGlasso, missknn, and MissForest even with 30% of censored data.

5. Application to single cell-data: megakaryocyte-erythroid progenitors

Recent advances in single-cell techniques have provided the opportunity to finely dissect cellular heterogeneity within known populations and to uncover rare cell types. In a study about the formation of blood cells, Psaila and others (2016) have recently identified three distinct subpopulations of cells, which are all derived from hematopoietic stem cells through cell differentiation. One of these sub-populations, denoted by MK-MEP, is a previously unknown, rare population of cells that are bipotent but primarily generate megakaryocytic progeny. In this section, we look closely at this sub-population and investigate the molecular mechanisms of regulation within these cells. To this end, we used data available from Psaila and others (2016) on 87 genes and 48 single human MK-MEP cells profiled by multiplex RT-qPCR.

As discussed in the Introduction, RT-qPCR data are typically right-censored. In this particular study, the limit of detection is fixed by the manufacturer to 40. Raw data have been mean normalized using the method proposed in Pipelers and others (2017) and using the two housekeeping genes provided. Figure 2a shows the relationship between the proportion of right-censored data and the mean of the normalized cycle-threshold. Two main conclusions can be drawn from this figure. Firstly, the proportion of censored data is an increasing function of the mean of the normalized cycle-threshold: this is expected and it means that the assumption of missing-at-random is not justified on this dataset. Secondly, there are some genes with a very high proportion of censoring (see the points above the dashed line): these may be genes whose transcription failed to amplify, in which case they should be treated as missing rather than censored. Following this explorative analysis and considering also the possible computational problems caused by the inclusion of these genes, we filtered out the genes with a proportion of censoring above 85% for subsequent analysis. The resulting data set contains the expression of 63 genes measured on 48 cells.

Fig. 2.

Open in new tabDownload slide

Real data analysis. (a) The proportion of right-censored data versus the mean of the normalized cycle-threshold; the black line is obtained by fitting a logistic regression model, while the dashed line identifies the threshold used to filter out the genes from the study. (b) Path of the ${\overline{\text{BIC}}}_{0.5}({\hat{\mathcal{E}}}^{\rho })$ measure; the vertical dashed line identifies the optimal value of the tuning parameter. (c) The estimated graph.

The normalized cyclic thresholds of the remaining genes are used to fit a right-censored GGM by using the proposed cglasso estimator. Figure 2b shows the path of the ${\overline{\text{BIC}}}_{0.5}({\hat{\mathcal{E}}}^{\rho })$ measure and the vertical dashed line is traced in correspondence of the optimal $\rho$-value. The resulting estimated graph (see Figure 2c) contains about $0.6\mathrm{\%}$ of all possible edges and about $19\mathrm{\%}$ of the considered genes. Consistent with the existing knowledge about this subpopulation, the estimated graph shows the central role of two genes, CD42 (glycoprotein 1b) and MYB (Psaila and others, 2016). The first one, CD42, is a megakaryocyte gene whose expression was found to be significantly high in the MK-MEP subpopulation. In particular, it is expressed later during megakaryocyte differentiation and has been associated with unipotent megakaryopoietic activity in mouse models. The link between CD42 and VWF, another highly expressed genes in this subpopulation (Psaila and others, 2016), was discovered by Chan and others (2017) as well as other connections, such as the link between the genes CD42 and CD61 and that between CD42 and TGFB1. The second central gene, MYB, is a transcription factor that is known to enhance erythroid differentiation at the expense of megakaryopoiesis.

6. Conclusion

In this article, we have proposed a computational approach to fit Gaussian Graphical models in the presence of censored data. The approach includes both the cases of right- and left-censored data. Since classical Gaussian graphical models cannot be used in high-dimensional settings, we also introduced ${\ell }_1$-penalization to produce sparsity (model selection) and parameter estimation simultaneously. The resulting estimator is called censored graphical lasso (cglasso). The computational problem of estimating the mean and conditional independence graph of a Gaussian graphical model from censored data is solved via an EM-algorithm. An extensive simulation study showed that the proposed estimator overcomes the existing estimators both in terms of parameter estimation and of network recovery. The analysis of a real RT-qPCR dataset showed how the method is able to infer the regulatory network underlying blood development under high levels of censoring.

In addition, we have proposed an approximated EM-algorithm, which is computationally more efficient and makes the method feasible for high-dimension settings. This approach relies on an efficient approximation of the mixed moments of the latent variables conditional on the observed. The approximation is exact under conditional independence and has thus been found to work well under sparse settings, as shown by Guo and others (2015) and our own simulations. Future work will investigate the theoretical justifications for this as well as refining the scenarios under which a good performance is to be expected.

Acknowledgments

Conflict of Interest: None declared.

Funding

The project was partially supported by the European Cooperation in Science and Technology (COST) [COST Action CA15109 European Cooperation for Statistics of Network Data Science (COSTNET)].

References