# Variable Selection and Joint Estimation of Mean and Covariance Models with an Application to eQTL Data

Authors: JungJun Lee; SungHwan Kim; Jae-Hwan Jhong; Ja-Yong Koo
Journal: Computational and Mathematical Methods in Medicine
 (2018)

---

## Abstract

In genomic data analysis, it is commonplace that underlying regulatory relationship over multiple genes is hardly ascertained due to unknown genetic complexity and epigenetic regulations. In this paper, we consider a joint mean and constant covariance model (JMCCM) that elucidates conditional dependent structures of genes with controlling for potential genotype perturbations. To this end, the modified Cholesky decomposition is utilized to parametrize entries of a precision matrix. The JMCCM maximizes the likelihood function to estimate parameters involved in the model. We also develop a variable selection algorithm that selects explanatory variables and Cholesky factors by exploiting the combination of the GCV and BIC as benchmarks, together with Rao and Wald statistics. Importantly, we notice that sparse estimation of a precision matrix (or equivalently gene network) is effectively achieved via the proposed variable selection scheme and contributes to exploring significant hub genes shown to be concordant toa priori biological evidence. In simulation studies, we confirm that our model selection efficiently identifies the true underlying networks. With an application to miRNA and SNPs data from yeast (a.k.a. eQTL data), we demonstrate that constructed gene networks reproduce validated biological and clinical knowledge with regard to various pathways including the cell cycle pathway.


---

## Body

## 1. Introduction

Generally, joint estimation of mean and covariance has been developed to address problems related to biomedical data. In longitudinal data analysis, for instance, identifying correct correlation structures within each subject is a major focus and many studies come up with variants of joint mean and covariance estimation to enhance statistical efficiency (see Ye and Pan [1] and references therein). With regard to graphical models, particularly, conditional Gaussian graphical models (cGGM) aiming to elucidate conditional dependence structures subject to the mean have increasingly received much attention [2–5] and conceptually can also be viewed as examples of joint mean and covariance estimation in essence. On a methodological side, joint estimation has been found to be practically applicable in the sense that applications of joint estimation methods cover bioinformatics, such as hormone and transcriptome [1, 4, 6–8], and stock price analysis [9]. A majority of existing methods involve variable selection of covariates for the mean vector and sparsity of covariance mainly based on the penalization methods. However, in this paper, we purposely focus on variable selection and joint estimation with a pure implementation of the likelihood-based method, unlike common penalization techniques, to better analyze eQTL data in light of genetic feature selection.Over the decades, genomics research has focused on comprehensive understanding of regulatory networks in the context of system biology. Commonly we are interested in a gene network, which pictures the interplay among genetic factors (e.g., gene regulation and activation). Particularly, it is important to investigate how a given genotype (genetic variants) at a particular quantitative trait locus (QTL) affects measured phenotypes and traits at that locus. For instance, gene expression quantitative loci (eQTL) make use of gene expressions as quantitative traits. eQTL analysis has been widely applied to figure out the effect of genetic perturbations associated with diseases as well as to construct regulatory networks describing how genes regulate expressions of other genes [10, 11]. More precisely, the location of single nucleotide polymorphisms (SNPs) may affect multiple gene expression levels, and this accidentally causes misleading inference for dependency structure among genes [4]. Many popular methods [2–5, 12] have been introduced to identify gene networks, aiming at learning networks subject to perturbation effects by genetic variants on the basis of population gene expression and genotype data.Yin and Li [4] proposed the sparse conditional graphical Gaussian model with l1 and the adaptive lasso penalty function. Liet al. [3] suggested the two-stage estimation framework: (1) estimating a nonsparse conditional covariance matrix of genes based on a conditional variance operator between the reproducing kernel Hilbert spaces of marker genes and then (2) using l1 and adaptive lasso penalty to obtain sparse estimates of a precision matrix under the cGGM. Caiet al. [12] studied the covariate-adjusted precision matrix estimation (CAPME) method using the constrained l1 optimization.While most recent studies encourage sparsity estimation of a precision matrix based on the penalized likelihood, we instead rely on classical variable methods based on the standard likelihood. Strictly speaking, the penalized likelihood, according to its definition, cannot be viewed as the likelihood. It naturally poses a question whether the likelihood-based method performs better than the penalized likelihood approaches. To address this question, we consider a joint mean and constant covariance model (JMCCM) inspired by Pourahmadi [13] and propose methods for variable selection which effectively identify sparse conditional gene networks and covariates relevant to gene regulations. We employ the modified Cholesky decomposition to guarantee positive definiteness of an estimated precision matrix [13]. With this reparametrization of the precision matrix, the log-likelihood function corresponding to our model can be decomposed into an additive form of each response in terms of Cholesky parameters. This facilitates a coordinate descent type implementation we adopt for the precision matrix estimation. The combination of both generalized cross-validation (GCV) and Bayesian information criterion (BIC) performs variable selection as a benchmark. Rao and Wald statistics are also utilized to add and delete genetic markers for each gene expression and the Cholesky factors for the precision matrix. To the best of our knowledge, only a few works use Rao and Wald statistics for variable selection in the joint estimation problem, particularly with an application to eQTL analysis.In simulation studies, extensive simulation scenarios experimentally confirm improved estimation efficiency and precise variable selection of the proposed method. For real data applications, we perform eQTL analysis via the JMCCM with gene expressions and SNPs of yeast data, in pursuit of detecting the effects of variant perturbations and underlying gene regulations. We find that the JMCCM effectively uncovers biological pathways that may potentially account for known biological processes. Taken together, the JMCCM is shown to be effective in identifying conditional dependence structures among variables compared to the existing graphical models with penalization methods such as the sparse cGGM [5] and CAPME [12].The paper is outlined as follows. In Section2, we describe our JMCCM with modified Cholesky decomposition and maximum likelihood estimates. The variable selection algorithm using Rao and Wald statistics for SNPs and Cholesky factors is explained in Section 3. Section 4 deals with simulation studies to demonstrate performance of variable selection and estimation of the proposed model. The yeast cell cycle pathway genes with SNPs data analysis are presented in Section 5. Concluding remarks and discussion in Section 6 are followed by the Appendix, which provides mathematical details of our method.


## 2. Model And Estimators

### 2.1. JMCCM with the Modified Cholesky Decomposition

Contrary to previous methods [3, 4], the JMCCM primarily aims at simultaneous estimation over the mean and precision matrix of the Gaussian graphical model. Suppose that (x,y) is a pair of the p×1 vector x of genetic markers and the m×1 random vector y of expression levels. Let x denote the vector of SNPs and let y denote the gene expression traits following m-variate multivariate normal distribution with the mean μ(x;β) and covariance matrix Σ as follows:(1)y∣x~Nmμx;β,Σ,where the jth entry of μ(x;β) is x⊤βj for j=1,…,m, and βj=(βj1,…,βjp) is the p×1 linear regression coefficient vector indicating effects of SNPs perturbations to gene expressions and β=(β1,…,βm). Importantly, note that Σ does not depend on x. The coefficient β is assumed to be sparse, since each gene is known to have only a few genetic regulators according to Caiet al. [12]. The precision matrix represents a gene network (graph), as in an undirected GGM, by corresponding the nonzero (i,j)th element of the precision matrix with an edge between two vertices i and j [14]. This edge represents conditional dependence of genes i and j given all other gene expression levels. Thus, our goal is to identify nonzero entries of the precision matrix in order to construct a conditional dependence genetic network after the effects of SNPs perturbations are removed. The precision matrix Σ-1 is also expected to be sparse [4].One of our primary interests is to estimate the precision matrixΣ-1, which is symmetrically positive definite, so we need to ensure that the estimate of the precision matrix also satisfies symmetrical positive definiteness. To this end, we apply the modified Cholesky decomposition [13] to the precision matrix, denoted by K, as follows:(2)Kϕ,τ=Cϕ⊤DτCϕ,where C is an upper triangular matrix with diagonal entries 1 and above-diagonal elements consisting of negative of ϕ=(ϕ12,…,ϕ1m,ϕ23,…,ϕ2m,…,ϕm-1,m) and D is a diagonal matrix containing τ=(τ1,…,τm) with τj>0 for all j=1,…,m as diagonal entries. Here, the superscript "⊤" denotes the transpose of a matrix. Positive definiteness of K is shown in Appendix A. Throughout this paper, a vector in the parenthesis is considered as a column vector. Let ϕj=(ϕj,j+1,…,ϕjm) for j=1,…,m-1. Then we can write ϕ=(ϕ1,…,ϕm-1). The parameter space for JMCCM is defined by Θ={θ=(β,ϕ,τ):β∈Rmp,ϕ∈Rm(m-1)/2andτ∈R+m}, where R+m represents the set of m-dimensional vectors of positive real numbers.


### 2.2. Maximum Likelihood Estimation

Suppose that we have theN independent observations (xi,yi), i=1,…,N, sampled from (1), where yi=(y1i,…,ymi) and xi=(x1i,…,xpi) represent the ith observation of x and y, respectively. With Σ-1=K=K(ϕ,τ), μ(x)=μ(x;β), and μi=μ(xi), the log-likelihood function corresponding to (1) is given by(3)lθ=N2-mlog⁡2π+∑j=1mlog⁡τj-trKVβ,where V(β)=1/N∑i=1N(yi-μi)yi-μi⊤ and "tr" denotes the trace of a matrix. A derivation of (3) can be found in Appendix B. The maximum likelihood estimator (MLE) of θ is defined by (4)θ^=argmaxθ∈Θ⁡lθ.Some notations are needed to express MLE ofθ in conjunction with variable selection. For j=1,…,m, let Gj be an index set of active (or significant) coefficients of βj and denote by gj the number of elements in Gj. Let G={G1,…,Gm} and g=g1+⋯+gm. Obviously, g≤mp and if g=mp, then all variables are significant for all y1,…,ym. Define β(G)=β1(G1),…,βm(Gm)∈Rg and βj(Gj)=βjll∈Gj∈Rgj for j=1,…,m. Write yk=(yk1,…,ykN) for k=1,…,m, X=1/N(x1,…,xp), and X(Gj)=[Xl]l∈Gj∈RN×gj, where Xl denotes lth column of X. Let Y=(y1⊤,…,ym⊤)∈RmN, X=⊕j=1mX(Gj)∈RmN×g, and W=K⊗IN∈RmN×mN, where ⊕ and ⊗ represent direct sum and Kronecker product, respectively, and IN is the N×N identity matrix. Then the MLE β^(G) of β(G) can be expressed as(5)β^G=X⊤WX-1X⊤WY.Observe that estimating β(G) involves K.We can estimateK by obtaining MLE of ϕ and τ. For j=1,…,m-1, let θj=(β,ϕj,τj) with θm=(β,τm). Let Vjj, Vj,21, and Vj,22 represent the (1,1), (2,1), and (2,2) components of the lower jth principal submatrix Vj(β) of V(β). Then we can express(6)lθ=∑j=1mljθj,where (7)ljθj=N2-log2π+log⁡τj-τjVjj-2ϕj⊤Vj,21+ϕj⊤Vj,22ϕj. By (6), we can obtain (ϕ^,τ^) from (ϕ^j,τ^j) for each j with ϕ^j=argmaxϕjlj(θj) and τ^j=argmaxτjlj(θj). Thus, instead of finding a solution to the optimization problem with (3), we optimize each lj(θj) with respect to ϕj and τj. Thus, the MLE ϕ^j and τ^j are expressed in terms of V(β) in a way that(8)ϕ^j=Vj,22-1Vj,21,1τ^j=Vjj-Vj,12Vj,22-1Vj,21forj=1…,m-1with β^ from (5) and 1/τ^m=Vmm. We estimate Σ-1 by K^ and it is defined by (9)K^=Kϕ^,τ^. Derivations of (5), (6), and (8) are presented in Appendix C.


## 2.1. JMCCM with the Modified Cholesky Decomposition

Contrary to previous methods [3, 4], the JMCCM primarily aims at simultaneous estimation over the mean and precision matrix of the Gaussian graphical model. Suppose that (x,y) is a pair of the p×1 vector x of genetic markers and the m×1 random vector y of expression levels. Let x denote the vector of SNPs and let y denote the gene expression traits following m-variate multivariate normal distribution with the mean μ(x;β) and covariance matrix Σ as follows:(1)y∣x~Nmμx;β,Σ,where the jth entry of μ(x;β) is x⊤βj for j=1,…,m, and βj=(βj1,…,βjp) is the p×1 linear regression coefficient vector indicating effects of SNPs perturbations to gene expressions and β=(β1,…,βm). Importantly, note that Σ does not depend on x. The coefficient β is assumed to be sparse, since each gene is known to have only a few genetic regulators according to Caiet al. [12]. The precision matrix represents a gene network (graph), as in an undirected GGM, by corresponding the nonzero (i,j)th element of the precision matrix with an edge between two vertices i and j [14]. This edge represents conditional dependence of genes i and j given all other gene expression levels. Thus, our goal is to identify nonzero entries of the precision matrix in order to construct a conditional dependence genetic network after the effects of SNPs perturbations are removed. The precision matrix Σ-1 is also expected to be sparse [4].One of our primary interests is to estimate the precision matrixΣ-1, which is symmetrically positive definite, so we need to ensure that the estimate of the precision matrix also satisfies symmetrical positive definiteness. To this end, we apply the modified Cholesky decomposition [13] to the precision matrix, denoted by K, as follows:(2)Kϕ,τ=Cϕ⊤DτCϕ,where C is an upper triangular matrix with diagonal entries 1 and above-diagonal elements consisting of negative of ϕ=(ϕ12,…,ϕ1m,ϕ23,…,ϕ2m,…,ϕm-1,m) and D is a diagonal matrix containing τ=(τ1,…,τm) with τj>0 for all j=1,…,m as diagonal entries. Here, the superscript "⊤" denotes the transpose of a matrix. Positive definiteness of K is shown in Appendix A. Throughout this paper, a vector in the parenthesis is considered as a column vector. Let ϕj=(ϕj,j+1,…,ϕjm) for j=1,…,m-1. Then we can write ϕ=(ϕ1,…,ϕm-1). The parameter space for JMCCM is defined by Θ={θ=(β,ϕ,τ):β∈Rmp,ϕ∈Rm(m-1)/2andτ∈R+m}, where R+m represents the set of m-dimensional vectors of positive real numbers.


## 2.2. Maximum Likelihood Estimation

Suppose that we have theN independent observations (xi,yi), i=1,…,N, sampled from (1), where yi=(y1i,…,ymi) and xi=(x1i,…,xpi) represent the ith observation of x and y, respectively. With Σ-1=K=K(ϕ,τ), μ(x)=μ(x;β), and μi=μ(xi), the log-likelihood function corresponding to (1) is given by(3)lθ=N2-mlog⁡2π+∑j=1mlog⁡τj-trKVβ,where V(β)=1/N∑i=1N(yi-μi)yi-μi⊤ and "tr" denotes the trace of a matrix. A derivation of (3) can be found in Appendix B. The maximum likelihood estimator (MLE) of θ is defined by (4)θ^=argmaxθ∈Θ⁡lθ.Some notations are needed to express MLE ofθ in conjunction with variable selection. For j=1,…,m, let Gj be an index set of active (or significant) coefficients of βj and denote by gj the number of elements in Gj. Let G={G1,…,Gm} and g=g1+⋯+gm. Obviously, g≤mp and if g=mp, then all variables are significant for all y1,…,ym. Define β(G)=β1(G1),…,βm(Gm)∈Rg and βj(Gj)=βjll∈Gj∈Rgj for j=1,…,m. Write yk=(yk1,…,ykN) for k=1,…,m, X=1/N(x1,…,xp), and X(Gj)=[Xl]l∈Gj∈RN×gj, where Xl denotes lth column of X. Let Y=(y1⊤,…,ym⊤)∈RmN, X=⊕j=1mX(Gj)∈RmN×g, and W=K⊗IN∈RmN×mN, where ⊕ and ⊗ represent direct sum and Kronecker product, respectively, and IN is the N×N identity matrix. Then the MLE β^(G) of β(G) can be expressed as(5)β^G=X⊤WX-1X⊤WY.Observe that estimating β(G) involves K.We can estimateK by obtaining MLE of ϕ and τ. For j=1,…,m-1, let θj=(β,ϕj,τj) with θm=(β,τm). Let Vjj, Vj,21, and Vj,22 represent the (1,1), (2,1), and (2,2) components of the lower jth principal submatrix Vj(β) of V(β). Then we can express(6)lθ=∑j=1mljθj,where (7)ljθj=N2-log2π+log⁡τj-τjVjj-2ϕj⊤Vj,21+ϕj⊤Vj,22ϕj. By (6), we can obtain (ϕ^,τ^) from (ϕ^j,τ^j) for each j with ϕ^j=argmaxϕjlj(θj) and τ^j=argmaxτjlj(θj). Thus, instead of finding a solution to the optimization problem with (3), we optimize each lj(θj) with respect to ϕj and τj. Thus, the MLE ϕ^j and τ^j are expressed in terms of V(β) in a way that(8)ϕ^j=Vj,22-1Vj,21,1τ^j=Vjj-Vj,12Vj,22-1Vj,21forj=1…,m-1with β^ from (5) and 1/τ^m=Vmm. We estimate Σ-1 by K^ and it is defined by (9)K^=Kϕ^,τ^. Derivations of (5), (6), and (8) are presented in Appendix C.


## 3. Consecutive Variable Selection Algorithm

As mentioned above,β and K are commonly believed to be sparse in genomic data analysis. To address sparsity, the lasso-type penalty imposed on both regression coefficients and precision matrix has been popularly applied to diverse graphical models [4, 5]. Stepping aside the lasso-type approach, we develop a variable selection technique that mainly relies on the combination of classical variable selection methods. Generally, the numbers of SNPs and genes tend to be considerably huge so that computational costs normally become prohibitive. In order to address this problem, the proposed variable selection algorithm proceeds with largely two stages: (1) preliminary variable selection for mean and precision matrix and (2) secondary variable selection in the middle of the joint model estimation. It is important to note that the first stage leads possible variables to be limited in scope (i.e., working parameters in the model) in order to circumvent high computational complexity.


### 3.1. Preliminary Variable Selection

Preliminary variable selection is largely twofold: variable selection for the mean part and covariance part. The idea behind that is to add variables (or equivalently parameters) to the joint model with the maximum Rao statistic and to delete ones with the minimum Wald statistic. You may refer to Koo [15] for the basis selection method or Kooperberget al. [16] that explain variable selection schemes based on Rao and Wald statistics.


#### 3.1.1. Mean Part

When it comes to the mean part, we carry out selecting predictor variables (i.e., SNPs) for each response variable (i.e., gene expression), dealing with a univariate multiple regression problem. In the addition stage, we start off with a model including only an intercept term. The MLE is used as estimator forβ and maximum Rao statistic (Rao [17]) is the criterion for adding a predictor together with GCV (Friedman [18]) as a stopping rule. In the deletion stage, Wald statistic is calculated to exclude predictor variables such that the updated model minimizes Wald statistics. The final model for the mean regression is chosen by the minimum GCV. Details are summarized in Algorithm 1. Once Algorithm 1 is done, the number of predictor variables included in the joint model no longer increases, while variable reduction can happen in Algorithm 3 (Section 3.2).Algorithm 1:Variable selection in the mean vector estimation.
( 1 )  for j in 1:mdo ( 2 ) Set the minimal linear regression model: Yj=βj0. ( 3 ) Compute β^j0 and GCV. ( 4 )  while GCV decreases do ( 5 ) Among predictors not in the current model, add a predictor Xb having the maximum Rao statistic. ( 6 ) Update β^j(Gj) using (5) and compute GCV. ( 7 )  end while ( 8 )  while GCV decreases do ( 9 ) Among predictors in the current model, delete a predictor Xb′ having the minimum Wald statistic. ( 10 ) Update β^j(Gj) using (5) and compute GCV. ( 11 )  end while ( 12 ) The optimal model is chosen by the minimum GCV. ( 13 ) end for


#### 3.1.2. Covariance Part

The rationale behind variable selection in the precision matrix estimation is that eachYj, j=1,…,m-1, is regressed on Yj+1,…,Ym, and the regression coefficients are negative of off-diagonal entries of C(ϕ), the by-product of the modified Cholesky decomposition (2) [13, 19, 20]. Clearly, this is one of the major benefits of reparametrization via Cholesky decomposition [13] in pursuit of improved interpretation.Subsequent to selection for predictor variables, we computeV(β^(G)). Given initial ϕ^j(0)=0, we start with computing τ^j for j=1,…,m as in (8) using V(β^(G)) obtained from Algorithm 1. Then compute Rao statistic to add one variable out of Yj+1…,Ym to the current model that builds on the response variable Yj. We repeatedly update for Yj, where j=1,…,m-1. Next, we choose one variable, say ϕlb, for l=1,…,m-1 and b=l+1,…,m, such that Rao statistic is maximized and thereby update ϕ^l and τ^l, respectively. When calculating τ^ and ϕ^, the BIC is computed, and the addition process stops if the BIC no longer decreases. At the completion of addition, we build a model with all selected variables as a full model and subsequently begin deletion from a full model. Deletion process is similar to the addition process except that the minimum Wald statistic is used for deletion in place of the maximum Rao statistic. Successive deletion continues until the BIC stops decreasing. In the last stage, the final model is also selected by the BIC. Details are presented in Algorithm 2. Once Algorithm 2 is finished, the number of ϕjk included in the joint model no longer increases, while variable reduction could occur in the joint estimation (Algorithm 3, Section 3.2). Sparsity of K^ is achieved through this variable selection scheme along with Algorithm 3.Algorithm 2:Variable selection in the precision matrix estimation.
( 1 ) Compute V(β^(G)) using β^(G) obtained from Algorithm 1. ( 2 ) Compute τ^1,…,τ^m by (8). ( 3 ) Combining initial ϕj(0)=0 for all j=1,…,m-1 with τ^, compute K^ and BIC. ( 4 )  while BIC decreasesdo ( 5 )  for j in 1:(m-1)do ( 6 )  if no elements of ϕj=(ϕj,j+1,…,ϕjm) are included in the joint model then ( 7 ) compute Rao statistics for ϕjk by Vkj2/Vkk, k=j+1,…,m. ( 8 )  else ( 9 ) Set a linear regression model with response Yj and predictors Yb’s whose corresponding coefficients ϕjb are already included in the joint model. Here b∈ { j + 1 , … , m }. Compute Rao statistic for adding one predictor among Yj+1,…,Ym whose corresponding ϕjk’s are not in this linear model. ( 10 )  end if ( 11 )  end for ( 12 ) Among all ϕjk’s not in the joint model, add one with the maximum Rao value. Denote this by ϕlb for l=1,…,m-1 and b=l+1,…,m. ( 13 ) Update ϕ^l, τ^l using (8) as well as K^ and compute BIC. ( 14 ) end while ( 15 ) while BIC decreases do ( 16 ) Compute Wald statistics for all ϕ in the current model. ( 17 ) Delete one ϕl′b′ with the minimum Wald. ( 18 ) Update ϕ^l′, τ^l′ using (8) as well as K^ and compute BIC. ( 19 ) end while ( 20 ) The optimal model is chosen by the minimum BIC.Algorithm 3:Variable selection in the joint estimation.
( 1 )  repeat ( 2 ) Update β^(G) by (5) and compute V(β^(G)) with K^ from Algorithm 2. ( 3 ) Update τ^, ϕ^ and K^ using V(β^(G)) obtained from the previous step. ( 4 )  until the log-likelihood function converges. ( 5 )  while BIC decreases do ( 6 ) Compute Wald statistic of deleting each βjp and ϕjb in the current model. ( 7 ) Delete one variable (i.e., parameter) with the minimum Wald statistic. ( 8 )  repeat ( 9 ) Update β^(G) by (5) and compute V(β^(G)) with K^ from the previous step. ( 10 ) Update τ^, ϕ^ and K^ using V(β^(G)) obtained from the previous step. ( 11 )  until the log-likelihood function converges. ( 12 ) end while ( 13 ) The optimal model is chosen by BIC.


### 3.2. Secondary Variable Selection in the Middle of the Joint Model Estimation

With variables (i.e., parameters) selected by the preliminary stage above, we implement simultaneous parameter estimation for both mean and precision in the joint model and additional variable selection. OnceK^ is computed via Algorithm 2 with fixed β^(G), we start joint estimation by updating β^(G), which is formed with weight least square estimates as in (5) and weights coming from K^. Then ϕ^ and τ^, ultimately K^, are newly computed using updated V(β^(G)), and again this updated K^ serves as a weight for updating β^(G). Over the updates, the interplay between V(β^(G)) and K^ continues until the log-likelihood function converges. Afterwards, deletion to current parameters begins by Wald statistic, excluding one parameter with the minimum Wald statistic. Finally, the BIC is used as a stopping rule for deletion and selection to finalize model estimation. Algorithm 3 contains details about this procedure. While the preliminary variable selection works for each mean and covariance part under the assumption of β^(G) or K^ are fixed, this joint estimation procedure is designed to improve estimation and selection accuracy by reflecting the changes of β^(G) and K^.


## 3.1. Preliminary Variable Selection

Preliminary variable selection is largely twofold: variable selection for the mean part and covariance part. The idea behind that is to add variables (or equivalently parameters) to the joint model with the maximum Rao statistic and to delete ones with the minimum Wald statistic. You may refer to Koo [15] for the basis selection method or Kooperberget al. [16] that explain variable selection schemes based on Rao and Wald statistics.


### 3.1.1. Mean Part

When it comes to the mean part, we carry out selecting predictor variables (i.e., SNPs) for each response variable (i.e., gene expression), dealing with a univariate multiple regression problem. In the addition stage, we start off with a model including only an intercept term. The MLE is used as estimator forβ and maximum Rao statistic (Rao [17]) is the criterion for adding a predictor together with GCV (Friedman [18]) as a stopping rule. In the deletion stage, Wald statistic is calculated to exclude predictor variables such that the updated model minimizes Wald statistics. The final model for the mean regression is chosen by the minimum GCV. Details are summarized in Algorithm 1. Once Algorithm 1 is done, the number of predictor variables included in the joint model no longer increases, while variable reduction can happen in Algorithm 3 (Section 3.2).Algorithm 1:Variable selection in the mean vector estimation.
( 1 )  for j in 1:mdo ( 2 ) Set the minimal linear regression model: Yj=βj0. ( 3 ) Compute β^j0 and GCV. ( 4 )  while GCV decreases do ( 5 ) Among predictors not in the current model, add a predictor Xb having the maximum Rao statistic. ( 6 ) Update β^j(Gj) using (5) and compute GCV. ( 7 )  end while ( 8 )  while GCV decreases do ( 9 ) Among predictors in the current model, delete a predictor Xb′ having the minimum Wald statistic. ( 10 ) Update β^j(Gj) using (5) and compute GCV. ( 11 )  end while ( 12 ) The optimal model is chosen by the minimum GCV. ( 13 ) end for


### 3.1.2. Covariance Part

The rationale behind variable selection in the precision matrix estimation is that eachYj, j=1,…,m-1, is regressed on Yj+1,…,Ym, and the regression coefficients are negative of off-diagonal entries of C(ϕ), the by-product of the modified Cholesky decomposition (2) [13, 19, 20]. Clearly, this is one of the major benefits of reparametrization via Cholesky decomposition [13] in pursuit of improved interpretation.Subsequent to selection for predictor variables, we computeV(β^(G)). Given initial ϕ^j(0)=0, we start with computing τ^j for j=1,…,m as in (8) using V(β^(G)) obtained from Algorithm 1. Then compute Rao statistic to add one variable out of Yj+1…,Ym to the current model that builds on the response variable Yj. We repeatedly update for Yj, where j=1,…,m-1. Next, we choose one variable, say ϕlb, for l=1,…,m-1 and b=l+1,…,m, such that Rao statistic is maximized and thereby update ϕ^l and τ^l, respectively. When calculating τ^ and ϕ^, the BIC is computed, and the addition process stops if the BIC no longer decreases. At the completion of addition, we build a model with all selected variables as a full model and subsequently begin deletion from a full model. Deletion process is similar to the addition process except that the minimum Wald statistic is used for deletion in place of the maximum Rao statistic. Successive deletion continues until the BIC stops decreasing. In the last stage, the final model is also selected by the BIC. Details are presented in Algorithm 2. Once Algorithm 2 is finished, the number of ϕjk included in the joint model no longer increases, while variable reduction could occur in the joint estimation (Algorithm 3, Section 3.2). Sparsity of K^ is achieved through this variable selection scheme along with Algorithm 3.Algorithm 2:Variable selection in the precision matrix estimation.
( 1 ) Compute V(β^(G)) using β^(G) obtained from Algorithm 1. ( 2 ) Compute τ^1,…,τ^m by (8). ( 3 ) Combining initial ϕj(0)=0 for all j=1,…,m-1 with τ^, compute K^ and BIC. ( 4 )  while BIC decreasesdo ( 5 )  for j in 1:(m-1)do ( 6 )  if no elements of ϕj=(ϕj,j+1,…,ϕjm) are included in the joint model then ( 7 ) compute Rao statistics for ϕjk by Vkj2/Vkk, k=j+1,…,m. ( 8 )  else ( 9 ) Set a linear regression model with response Yj and predictors Yb’s whose corresponding coefficients ϕjb are already included in the joint model. Here b∈ { j + 1 , … , m }. Compute Rao statistic for adding one predictor among Yj+1,…,Ym whose corresponding ϕjk’s are not in this linear model. ( 10 )  end if ( 11 )  end for ( 12 ) Among all ϕjk’s not in the joint model, add one with the maximum Rao value. Denote this by ϕlb for l=1,…,m-1 and b=l+1,…,m. ( 13 ) Update ϕ^l, τ^l using (8) as well as K^ and compute BIC. ( 14 ) end while ( 15 ) while BIC decreases do ( 16 ) Compute Wald statistics for all ϕ in the current model. ( 17 ) Delete one ϕl′b′ with the minimum Wald. ( 18 ) Update ϕ^l′, τ^l′ using (8) as well as K^ and compute BIC. ( 19 ) end while ( 20 ) The optimal model is chosen by the minimum BIC.Algorithm 3:Variable selection in the joint estimation.
( 1 )  repeat ( 2 ) Update β^(G) by (5) and compute V(β^(G)) with K^ from Algorithm 2. ( 3 ) Update τ^, ϕ^ and K^ using V(β^(G)) obtained from the previous step. ( 4 )  until the log-likelihood function converges. ( 5 )  while BIC decreases do ( 6 ) Compute Wald statistic of deleting each βjp and ϕjb in the current model. ( 7 ) Delete one variable (i.e., parameter) with the minimum Wald statistic. ( 8 )  repeat ( 9 ) Update β^(G) by (5) and compute V(β^(G)) with K^ from the previous step. ( 10 ) Update τ^, ϕ^ and K^ using V(β^(G)) obtained from the previous step. ( 11 )  until the log-likelihood function converges. ( 12 ) end while ( 13 ) The optimal model is chosen by BIC.


## 3.1.1. Mean Part

When it comes to the mean part, we carry out selecting predictor variables (i.e., SNPs) for each response variable (i.e., gene expression), dealing with a univariate multiple regression problem. In the addition stage, we start off with a model including only an intercept term. The MLE is used as estimator forβ and maximum Rao statistic (Rao [17]) is the criterion for adding a predictor together with GCV (Friedman [18]) as a stopping rule. In the deletion stage, Wald statistic is calculated to exclude predictor variables such that the updated model minimizes Wald statistics. The final model for the mean regression is chosen by the minimum GCV. Details are summarized in Algorithm 1. Once Algorithm 1 is done, the number of predictor variables included in the joint model no longer increases, while variable reduction can happen in Algorithm 3 (Section 3.2).Algorithm 1:Variable selection in the mean vector estimation.
( 1 )  for j in 1:mdo ( 2 ) Set the minimal linear regression model: Yj=βj0. ( 3 ) Compute β^j0 and GCV. ( 4 )  while GCV decreases do ( 5 ) Among predictors not in the current model, add a predictor Xb having the maximum Rao statistic. ( 6 ) Update β^j(Gj) using (5) and compute GCV. ( 7 )  end while ( 8 )  while GCV decreases do ( 9 ) Among predictors in the current model, delete a predictor Xb′ having the minimum Wald statistic. ( 10 ) Update β^j(Gj) using (5) and compute GCV. ( 11 )  end while ( 12 ) The optimal model is chosen by the minimum GCV. ( 13 ) end for


## 3.1.2. Covariance Part

The rationale behind variable selection in the precision matrix estimation is that eachYj, j=1,…,m-1, is regressed on Yj+1,…,Ym, and the regression coefficients are negative of off-diagonal entries of C(ϕ), the by-product of the modified Cholesky decomposition (2) [13, 19, 20]. Clearly, this is one of the major benefits of reparametrization via Cholesky decomposition [13] in pursuit of improved interpretation.Subsequent to selection for predictor variables, we computeV(β^(G)). Given initial ϕ^j(0)=0, we start with computing τ^j for j=1,…,m as in (8) using V(β^(G)) obtained from Algorithm 1. Then compute Rao statistic to add one variable out of Yj+1…,Ym to the current model that builds on the response variable Yj. We repeatedly update for Yj, where j=1,…,m-1. Next, we choose one variable, say ϕlb, for l=1,…,m-1 and b=l+1,…,m, such that Rao statistic is maximized and thereby update ϕ^l and τ^l, respectively. When calculating τ^ and ϕ^, the BIC is computed, and the addition process stops if the BIC no longer decreases. At the completion of addition, we build a model with all selected variables as a full model and subsequently begin deletion from a full model. Deletion process is similar to the addition process except that the minimum Wald statistic is used for deletion in place of the maximum Rao statistic. Successive deletion continues until the BIC stops decreasing. In the last stage, the final model is also selected by the BIC. Details are presented in Algorithm 2. Once Algorithm 2 is finished, the number of ϕjk included in the joint model no longer increases, while variable reduction could occur in the joint estimation (Algorithm 3, Section 3.2). Sparsity of K^ is achieved through this variable selection scheme along with Algorithm 3.Algorithm 2:Variable selection in the precision matrix estimation.
( 1 ) Compute V(β^(G)) using β^(G) obtained from Algorithm 1. ( 2 ) Compute τ^1,…,τ^m by (8). ( 3 ) Combining initial ϕj(0)=0 for all j=1,…,m-1 with τ^, compute K^ and BIC. ( 4 )  while BIC decreasesdo ( 5 )  for j in 1:(m-1)do ( 6 )  if no elements of ϕj=(ϕj,j+1,…,ϕjm) are included in the joint model then ( 7 ) compute Rao statistics for ϕjk by Vkj2/Vkk, k=j+1,…,m. ( 8 )  else ( 9 ) Set a linear regression model with response Yj and predictors Yb’s whose corresponding coefficients ϕjb are already included in the joint model. Here b∈ { j + 1 , … , m }. Compute Rao statistic for adding one predictor among Yj+1,…,Ym whose corresponding ϕjk’s are not in this linear model. ( 10 )  end if ( 11 )  end for ( 12 ) Among all ϕjk’s not in the joint model, add one with the maximum Rao value. Denote this by ϕlb for l=1,…,m-1 and b=l+1,…,m. ( 13 ) Update ϕ^l, τ^l using (8) as well as K^ and compute BIC. ( 14 ) end while ( 15 ) while BIC decreases do ( 16 ) Compute Wald statistics for all ϕ in the current model. ( 17 ) Delete one ϕl′b′ with the minimum Wald. ( 18 ) Update ϕ^l′, τ^l′ using (8) as well as K^ and compute BIC. ( 19 ) end while ( 20 ) The optimal model is chosen by the minimum BIC.Algorithm 3:Variable selection in the joint estimation.
( 1 )  repeat ( 2 ) Update β^(G) by (5) and compute V(β^(G)) with K^ from Algorithm 2. ( 3 ) Update τ^, ϕ^ and K^ using V(β^(G)) obtained from the previous step. ( 4 )  until the log-likelihood function converges. ( 5 )  while BIC decreases do ( 6 ) Compute Wald statistic of deleting each βjp and ϕjb in the current model. ( 7 ) Delete one variable (i.e., parameter) with the minimum Wald statistic. ( 8 )  repeat ( 9 ) Update β^(G) by (5) and compute V(β^(G)) with K^ from the previous step. ( 10 ) Update τ^, ϕ^ and K^ using V(β^(G)) obtained from the previous step. ( 11 )  until the log-likelihood function converges. ( 12 ) end while ( 13 ) The optimal model is chosen by BIC.


## 3.2. Secondary Variable Selection in the Middle of the Joint Model Estimation

With variables (i.e., parameters) selected by the preliminary stage above, we implement simultaneous parameter estimation for both mean and precision in the joint model and additional variable selection. OnceK^ is computed via Algorithm 2 with fixed β^(G), we start joint estimation by updating β^(G), which is formed with weight least square estimates as in (5) and weights coming from K^. Then ϕ^ and τ^, ultimately K^, are newly computed using updated V(β^(G)), and again this updated K^ serves as a weight for updating β^(G). Over the updates, the interplay between V(β^(G)) and K^ continues until the log-likelihood function converges. Afterwards, deletion to current parameters begins by Wald statistic, excluding one parameter with the minimum Wald statistic. Finally, the BIC is used as a stopping rule for deletion and selection to finalize model estimation. Algorithm 3 contains details about this procedure. While the preliminary variable selection works for each mean and covariance part under the assumption of β^(G) or K^ are fixed, this joint estimation procedure is designed to improve estimation and selection accuracy by reflecting the changes of β^(G) and K^.


## 4. Experimental Studies

In order to assess the performance of our proposed method, we carry out experimental simulations and compare the sparse conditional Gaussian graphical model (SCGGM), Zhang and Kim [5], covariate-adjusted precision matrix estimation (CAPME), Caiet al. [12] and joint model with lasso penalty (JML), and Jhonget al. [21], all of which are based on penalized likelihood approaches. The competing methods are run with their default setting regarding tuning parameters. To evaluate similarity, the estimated precision matrix and true matrix are benchmarked by the Steins loss function: (10)δSteinK,K^=trKK^-1-logKK^-1-m,where K^ is an estimate of the true precision matrix K and · denotes the determinant of a matrix. The Frobenius norm of difference between K and K^, denoted by ΔF, where Δ=K-K^, is also considered. In addition to the Steins loss function, to measure how efficiently our model recovers the true conditional dependent relationship among genes, specificity (SPE) and sensitivity (SEN) are used, as defined by (11)SPE=TNTN+FP,SEN=TPTP+FN, where TN, TP, FN, and FP are the numbers of true negatives, true positives, false negatives, and false positives with regard to off-diagonal elements of a precision matrix. Here, we treat a nonzero entry of a precision matrix as “positive.” To combine sensitivity and specificity, Youden’s index (=SPE+SEN-1) is used. The smaller values of δStein and ΔF are better, whereas the larger values of SPE, SEN, and Youden are better.Inspired by Yin and Li [4], we generate simulation data sets in the form of eQTL data sets such that nonzero entries of a precision matrix, commonly called a link (or edge), are randomly assigned with probability c1/m, where m is the number of genes and c1 is some positive constant. For a link generated at the (i,j)th entry of the precision matrix, denoted by kij, the corresponding element is sampled from the uniform distribution over [-1,0.5]∪[0.5,1]. For each row, off-diagonal elements are divided by the sum of their absolute values multiplied by 1.5. And we obtain the true precision matrix K by symmetrizing and setting diagonals as 1. To create the p×m regression coefficients matrix β, we first generate a p×m indicator matrix that has 1 as entry with probability c2/p for some positive constant c2. If the (l,j)th element of this indicator matrix is 1, βjl is randomly generated from Unif([dm,1]∪[-1,-dm]), where dm is the smallest absolute value of K generated.Producingβ and K, we generate SNPs, X=(X1,…,Xp) with Xl~Bernoulli(1,0.5) for l=1,…,p. Finally, we simulated gene expressions by generating y from the multivariate normal distribution given x, Y∣X~N(Xβ,K-1). We generate a data set of N i.i.d. random vectors (X,Y), and simulations are repeated 50 times. In Table 1, we outline the six simulation scenarios of small-scale setup. Table 2 contains the three simulation scenarios of large-scale setup.Table 1
Six scenarios for small-scale experimental study.
Model m p P ( k i j ≠ 0 ) P ( β j l ≠ 0 ) 1 10 10 2 / m 3.5 / p 2 20 10 2 / m 3.5 / p 3 40 10 2 / m 3.5 / p 4 20 20 2 / m 4 / p 5 30 30 2 / m 4 / p 6 40 40 2 / m 4 / pTable 2
Three scenarios for large-scale experimental study.
Model m p P ( k i j ≠ 0 ) P ( β j l ≠ 0 ) 7 100 100 2 / m 3 / p 8 200 200 2.5 / m 15 / p 9 400 200 1.5 / m 20 / pThe small-scale simulation results in Table3 suggest that the JMCCM produces better estimates than all other methods across all six models. Due to computational issues, we drop some results of JML from Table 3. By comparing the results of model 3 with 6 and 2 with 4, when N=1000 and the number of genes is fixed, we can see that JMCCM and CAPME show less changes in δStein and ΔF than they appear in the SCGGM as the number of SNPs increases. This indicates that JMCCM and CAPME are less subject to the increment of the number of SNPs than SCGGM, when modest size of genes is involved. Estimation performance of CAPME seems to be affected by the number of genes more easily than JMCCM because Stein loss and Frobenius norm of CAPME for models 1, 2, and 3 with fixed p increase more rapidly than those of JMCCM do. For identifying structures of the precision matrix, JMCCM surpasses SCGGM in discovering nonzero elements (higher sensitivity) as complexity of the model rises. This is possibly due to the fact that SCGGM tends to produce sparse estimates more than the true precision matrix. Higher sensitivity implies that our proposed model is less likely to miss the influential conditional dependency among genes. While CAPME and JML score near 1 in sensitivity, they mark poor number in specificity because 5-fold cross-validation for CAPME and validation approach for JML with default setting for tuning parameters selection tend to choose small ones, leading to dense estimates. JMCCM produces higher Youden’s index across all simulation scenarios compared to all other methods, and the performance gap between JMCCM and others is increasingly widened as the models increase in sample size and complexity. The results for the large-scale simulations (models 7–9) are summarized in Table 4 and Table S3. Due to long computing time, the results of CAPME for models 8 and 9 are not reported. Overall, the results are consistent with the small-scale simulations. The gap between JMCCM and SCGGM in estimation performance widened as the complexity of the model increases (Table 4).Table 3
Comparisons of the performance of JMCCM with SCGGM, CAPME, and JML for models 1–6. Standard errors are presented in parenthesis.
Model Method δ S t e i n Δ F SPE SEN Youden N = 500 1 JMCCM 0.139 (0.004) 0.519 (0.008) 0.864 (0.010) 0.723 (0.014) 0.587 (0.012) SCGGM 0.171 (0.007) 0.556 (0.010) 0.865 (0.030) 0.624 (0.027) 0.506 (0.014) CAPME 0.113 (0.003) 0.512 (0.008) 0.002 (0.001) 1.000 (0.000) 0.002 (0.001) JML 0.099 (0.003) 0.450 (0.005) 0.160 (0.026) 0.980 (0.005) 0.140 (0.024)  N = 500 2 JMCCM 0.241 (0.005) 0.677 (0.009) 0.915 (0.004) 0.801 (0.007) 0.716 (0.006) SCGGM 0.349 (0.008) 0.821 (0.012) 0.661 (0.029) 0.882 (0.016) 0.542 (0.017) CAPME 0.430 (0.007) 1.028 (0.011) 0.005 (0.001) 1.000 (0.000) 0.005 (0.001) JML 0.232 (0.004) 0.670 (0.008) 0.210 (0.009) 0.980 (0.003) 0.200 (0.009)  N = 500 3 JMCCM 0.705 (0.011) 1.145 (0.009) 0.947 (0.002) 0.653 (0.006) 0.600 (0.005) SCGGM 1.158 (0.034) 1.506 (0.032) 0.716 (0.036) 0.706 (0.029) 0.422 (0.010) CAPME 1.703 (0.012) 2.217 (0.015) 0.005 (0.000) 1.000 (0.000) 0.005 (0.000) JML - - - - -  N = 500 4 JMCCM 0.239 (0.006) 0.678 (0.011) 0.910 (0.004) 0.802 (0.007) 0.712 (0.007) SCGGM 0.402 (0.010) 0.869 (0.011) 0.825 (0.031) 0.772 (0.017) 0.597 (0.019) CAPME 0.447 (0.007) 1.081 (0.010) 0.004 (0.001) 1.000 (0.000) 0.004 (0.001) JML 0.284 (0.005) 0.793 (0.008) 0.110 (0.006) 1.000 (0.001) 0.100 (0.006)  N = 500 5 JMCCM 0.463 (0.007) 0.953 (0.007) 0.941 (0.003) 0.731 (0.007) 0.672 (0.006) SCGGM 0.743 (0.011) 1.188 (0.010) 0.978 (0.002) 0.567 (0.007) 0.545 (0.007) CAPME 0.983 (0.008) 1.683 (0.012) 0.005 (0.001) 1.000 (0.000) 0.005 (0.001) JML 0.524 (0.524) 1.074 (0.009) 0.140 (0.006) 0.980 (0.002) 0.120 (0.006)  N = 500 6 JMCCM 0.711 (0.010) 1.153 (0.010) 0.951 (0.002) 0.639 (0.006) 0.590 (0.005) SCGGM 1.126 (0.011) 1.438 (0.007) 0.988 (0.001) 0.440 (0.004) 0.428 (0.004) CAPME 1.781 (0.013) 2.406 (0.013) 0.006 (0.000) 0.999 (0.000) 0.006 (0.000) JML - - - - -  N = 1000 1 JMCCM 0.068 (0.002) 0.364 (0.006) 0.814 (0.007) 0.886 (0.009) 0.700 (0.009) SCGGM 0.139 (0.010) 0.483 (0.018) 0.888 (0.024) 0.668 (0.031) 0.556 (0.016) CAPME 0.055 (0.002) 0.346 (0.006) 0.003 (0.001) 1.000 (0.000) 0.003 (0.001) JML 0.046 (0.001) 0.301 (0.005) 0.150 (0.015) 1.000 (0.001) 0.150 (0.015)  N = 1000 2 JMCCM 0.117 (0.003) 0.466 (0.006) 0.894 (0.003) 0.898 (0.005) 0.793 (0.004) SCGGM 0.168 (0.007) 0.550 (0.010) 0.652 (0.018) 0.960 (0.006) 0.612 (0.016) CAPME 0.215 (0.003) 0.696 (0.007) 0.003 (0.001) 1.000 (0.000) 0.003 (0.001) JML 0.119 (0.002) 0.459 (0.006) 0.320 (0.120) 0.990 (0.003) 0.310 (0.012)  N = 1000 3 JMCCM 0.321 (0.006) 0.766 (0.007) 0.916 (0.002) 0.808 (0.004) 0.724 (0.003) SCGGM 0.645 (0.015) 1.107 (0.011) 0.674 (0.021) 0.847 (0.024) 0.521 (0.008) CAPME 0.824 (0.006) 1.405 (0.007) 0.005 (0.000) 1.000 (0.000) 0.005 (0.000) JML - - - - -  N = 1000 4 JMCCM 0.115 (0.003) 0.467 (0.006) 0.893 (0.003) 0.904 (0.006) 0.797 (0.005) SCGGM 0.286 (0.009) 0.745 (0.012) 0.786 (0.028) 0.822 (0.019) 0.609 (0.015) CAPME 0.215 (0.003) 0.703 (0.006) 0.003 (0.001) 1.000 (0.000) 0.003 (0.001) JML 0.123 (0.002) 0.486 (0.005) 0.160 (0.006) 1.000 (0.002) 0.150 (0.005)  N = 1000 5 JMCCM 0.217 (0.004) 0.643 (0.006) 0.924 (0.002) 0.867 (0.005) 0.792 (0.004) SCGGM 0.671 (0.011) 1.141 (0.001) 0.995 (0.002) 0.542 (0.009) 0.536 (0.008) CAPME 0.473 (0.004) 1.065 (0.007) 0.004 (0.000) 1.000 (0.000) 0.004 (0.000) JML - - - - -  N = 1000 6 JMCCM 0.331 (0.005) 0.781 (0.006) 0.926 (0.002) 0.789 (0.005) 0.715 (0.004) SCGGM 1.024 (0.010) 1.380 (0.006) 0.998 (0.000) 0.414 (0.003) 0.412 (0.003) CAPME 0.844 (0.006) 1.466 (0.008) 0.005 (0.000) 0.999 (0.000) 0.004 (0.000) JML - - - - -Table 4
Comparisons of the performance of JMCCM with SCGGM and CAPME for models 7–9. Standard errors are presented in parenthesis.
Model Method δ S t e i n Δ F SPE SEN Youden N = 500 7 JMCCM 2.234 (0.021) 2.071 (0.011) 0.969 (0.001) 0.638 (0.004) 0.608 (0.003) SCGGM 2.670 (0.017) 2.160 (0.007) 0.994 (0.000) 0.512 (0.002) 0.506 (0.002) CAPME 12.914 (0.038) 9.257 (0.028) 0.002 (0.000) 1.000 (0.000) 0.002 (0.000)  N = 500 8 JMCCM 6.560 (0.073) 3.669 (0.028) 0.982 (0.001) 0.475 (0.004) 0.457 (0.003) SCGGM 10.687 (0.035) 4.318 (0.006) 0.925 (0.000) 0.430 (0.002) 0.355 (0.002) CAPME - - - - -  N = 500 9 JMCCM 16.703 (0.254) 5.998 (0.062) 0.986 (0.001) 0.466 (0.003) 0.452 (0.002) SCGGM 34.025 (0.722) 7.277 (0.048) 0.874 (0.005) 0.453 (0.003) 0.328 (0.003) CAPME - - - - -We also simulate SNPs to mimic the linkage disequilibrium (LD) which is known to be a common phenomenon in DNA sequence and to assess the performance of our approach. We randomly generate two groups of SNPs which lie on LD, each of which contains LD block including 10 correlated SNPs such that correlation is greater than 0.9. Together with these SNPs, 10 SNPs are generated independently from Bernoulli trials with probability of 0.5 for a total of 20 SNPs and 20 genes:m=20 and p=20. Table 5 presents the simulation results using these data of sample size N=500 and N=1000 with 50 repetitions. Compared to the results of model 4 in Table 3, which is based on 20 independent SNPs with N=500, LD is found to have little impact on all methods, except for slight decreases in Stein loss and Frobenius norm of SCGGM.Table 5
Comparisons of the performance of JMCCM with SCGGM and CAPME for the experiment with highly correlated SNPs. 20 genes (m) and 20 SNPs (p) are involved. Standard errors are presented in parenthesis.
Method δ S t e i n Δ F SPE SEN Youden N = 500 JMCCM 0.239 (0.006) 0.678 (0.011) 0.910 (0.004) 0.802 (0.007) 0.712 (0.007) SCGGM 0.402 (0.010) 0.869 (0.011) 0.825 (0.031) 0.772 (0.017) 0.597 (0.019) CAPME 0.447 (0.007) 1.081 (0.010) 0.004 (0.001) 1.000 (0.000) 0.0004 (0.001)  N = 1000 JMCCM 0.115 (0.003) 0.467 (0.006) 0.893 (0.003) 0.904 (0.006) 0.797 (0.005) SCGGM 0.286 (0.009) 0.745 (0.012) 0.786 (0.028) 0.822 (0.019) 0.609 (0.015) CAPME 0.215 (0.003) 0.703 (0.006) 0.003 (0.001) 1.000 (0.000) 0.003 (0.001)


## 5. eQTL Analysis of Yeast Data

In this section, we apply the proposed algorithm to genomic eQTL (i.e., expression quantitative trait loci) data in order to examine whether the proposed method effectively recovers true dependency of gene expressions, which builds on known molecular mechanisms. To this end, we collect a set of yeast data [22], including polymorphic genotypes and mRNA expressions. The yeast data have been widely applied to elucidate the biological interactions between nucleotide polymorphisms and their responding genes (e.g., perturbation effects) [23–26]. Thus, our primary goal is to identify conditional dependency among genes with an adjustment of SNPs perturbations to each gene expression level.The data sets are collected for two yeast parent strains, BY4716 (BY) and RM11-1a (RM), and their 112 segregants. We obtain SNPs for 1,260 loci to the exclusion of the redundant SNPs observed in neighboring genetic regions and leave 3,684 expression genes after screening out genes of missing more than 5%. In order to validate whether or not the estimated gene network is consistent with unknown biological knowledge, we take the true signaling pathway for comparison. Out of 3,684 yeast genes, we purposely focus on a set of 64 genes that are ascertained in the cell cycle yeast pathway available in the KEGG database [24]. Together with 1,260 SNPs as predictors, we construct a gene network model of 64 genes. JMCCM selects 222 nonzero elements of the precision matrix and leaves nonzero regression coefficients, for a total of 111 links among genes. A total of 489 regression coefficients of SNPs over gene expression levels are included in the final model. Figure 2 displays the conditional gene network estimated by the proposed joint model. A pair of two genes is linked if an off-diagonal element of the estimated precision matrix is nonzero, and so 51 genes are shown to be linked and assemble in a module.Given this estimated gene network, we notice that the gene network is found to be concordant with the true cell cycle pathway. For instance, according to the true cell cycle pathway in Figure1, MCM3 is linked to ORC3, ORC5, MCM7, and MCM4. MCM5 is connected to TAH11 and ORC1. PHO4 is closely linked to CLN2 and SIC1, whose molecular function is related to MAP kinase orthologs (i.e., MAPK pathway genes). It is interesting to note that CLN3 and SWI4 are mutually connected together, linking to their downstream gene MBP1. Hence, this undirected graph contributes to recovering lots of links among the 64 genes of the pathway.Figure 1
The yeast cell cycle pathway from the KEGG database. Source:http://www.kegg.jp/kegg-bin/highlight_pathway?scale=1.0&map=map04111&keyword=cell%20cycle.Figure 2
The gene network of the yeast data related to the cell cycle pathway via the JMCCM.With regard to genetic variant effects, JMCCM is shown to effectively identify some of the well-known direct genetic perturbations. Gene expressions are regulated by some genetic variants, which, unless otherwise taken into account, may falsely capture the interplay of genes in a network. Our founding includes Clb-specific Cdk inhibitor (SIC1) as influencing the molecular interface between cyclin-Cdk complexes (e.g., binding to and blocking Cdk1/Clb activity, ultimately to maneuver the timing of DNA replication (see TableS1 in Supplementary Materials; [27]). In addition, our gene networks demonstrate that SIC1 is strongly perturbed by CLB3, while SCGGM did not detect any perturbation effects. More interestingly, many previous experiments validated the association between SIC1 and CLB3 to uncover the underlying mechanism of the cell cycle [28–30]. Clearly, this implies that JMCCM does a better job in accounting for SNPs perturbation as compared to SCGGM.Moreover, the gene module detected from JMCCM is shown in Table6. Pertaining to these gene modules, we hypothesize if gene modules, each containing hub genes and their neighboring genes, are enriched with common biological processes or not. To test this hypothesis, we conduct the gene ontology (GO) enrichment analysis [31] over the detected gene network modules (see Table S2 in Supplementary Materials) from both JMCCM and SCGGM using Fisher’s exact tests. Table 6 demonstrates that the proposed method outstandingly performs detecting modules biologically associated with many molecular processes, while none is detected by SCGGM. More importantly, the pathway of the organic acid metabolic process enriched in gene module 1 is also reported and validated on the basis of the network modules constructed with large-scale integration of yeast data in Zhuet al. [32].Table 6
Gene ontology (GO) enrichment analysis over the genes in detected module from the JMCCM.
Module Module size GO category GO enrichmentp value Module 1 94 Carboxylic acid metabolic process 0.00101 Carboxypeptidase activity 0.00904 Catabolic process 0.00101 Cellular catabolic process 0.00101 Exopeptidase activity 0.00904 Metalloexopeptidase activity 0.00904 Metallopeptidase activity 0.00904 Organic acid metabolic process 0.00101 Peptidase activity 0.00904Putting all things together, we conclude that the proposed method facilitates recovering the true SNPs perturbations in the midst of gene regulations and elucidating the underlying interplay of gene interactions. These fortes of the proposed algorithm are favored for reinforcinga priori biological knowledge and address a novel hypothesis related to clinical and translational potential.


## 6. Conclusion and Discussion

In this paper, we propose JMCCM to efficiently identify conditional dependent structures of gene expressions with adjustments to perturbation effects of SNPs. Contrary to the existing conditional graphical models, the precision matrix commonly used to reveal the true relationship among genes is parameterized via the modified Cholesky decomposition. The maximum likelihood estimates of the precision matrix were computed, while variable selection of SNPs and Cholesky factors are carried out separately and jointly by the GCV and BIC criterion. From experimental studies, it is clearly shown that JMCCM performs better than the existing penalization methods. Besides, JMCCM in the application to yeast cell cycle data successfully recovers many parts of the cell cycle pathway with adjustments of SNPs to each gene expressions level. Notably, the model entails the estimation of precision matrix, of which components are assumed to be constant. So, in the future, we may relax this somewhat strong assumption in the way that the model can parametrize overτ and ϕ in pursuit of more accurate estimation [13]. We leave this for next study.



---

*Source: 4626307-2018-06-25.xml*