Variable Selection for High-dimensional Generalized Linear Models using an Iterated Conditional Modes/Medians Algorithm
Vitara Pungpapong Min Zhang Dabao Zhang
stat.ME






ChulaAddress]Vitara Pungpapong mycorrespondingauthor
vitara@cbs.chula.ac.th
[mycorrespondingauthor]Corresponding author.
PurdueAddress]Min Zhang
PurdueAddress]Dabao Zhang

[ChulaAddress]Department of Statistics, Faculty of Commerce and Accountancy, Chulalongkorn University, Bangkok, Thailand
[PurdueAddress]Department of Statistics, Purdue University, West Lafayette, IN, USA


High-dimensional linear and nonlinear models have been extensively used to identify associations between response and explanatory variables. The variable selection problem is commonly of interest in the presence of massive and complex data. An empirical Bayes model for high-dimensional generalized linear models (GLMs) is considered in this paper. The extension of the Iterated Conditional Modes/Medians (ICM/M) algorithm is proposed to build up a GLM. With the construction of pseudodata and pseudovariances based on iteratively reweighted least squares (IRLS), conditional modes are employed to obtain data-drive optimal values for hyperparameters and conditional medians are used to estimate regression coefficients. With a spike-and-slab prior for each coefficient, a conditional median can enforce variable estimation and selection at the same time. The ICM/M algorithm can also incorporate more complicated prior by taking the network structural information into account through the Ising model prior. Here we focus on two extensively used models for genomic data: binary logistic and Cox's proportional hazards models. The performance of the proposed method is demonstrated through both simulation studies and real data examples. The implementation of the ICM/M algorithm for both linear and nonlinear models can be found in the  R package which is freely available on CRAN.



Empirical Bayes variable selection; High-dimensional data; Generalized linear models; Binary logistic model; Cox's model.








§ INTRODUCTION

Due to the advanced development of computing technology, high-dimensional data are widely found in many fields ranging from genomics to marketing and finance. Such data are known as large p small n data where the number of variables p is relatively larger than the sample size n. Generalized linear models (GLMs) have been extensively used to identify associations between response and explanatory variables. Similar to a normal linear regression model, the classical GLMs approach also suffers in large p small n scenario due to collinearity among covariates. Furthermore, regression coefficients are usually assumed to be sparse. 

In biomedical science, variable selection plays an important role in analysis of high-throughput genotype. A typical microarray data have several thousands predictors. For genome-wide association studies (GWAS), the data may contain thousand to million Single nucleotide polymorphisms (SNPs). However, the sample size is usually only in hundreds or thousands. Collinearity problem occurs especially for those genes or SNPs resided in the same biological pathway. It is also believed that only few genes or SNPs are associated with the phenotype of interest. In addition to continuouse phenotype, discrete or censored outcomes frequently appear in the field of biomedical science. Binary logistic and Cox's proportional hazards models with high-dimensional covariates have received considerable attention over the past decade.

One popular approach in fitting high-dimensional GLMs is through a penalized likelihood method. The most basic and popular penalized likelihood estimator is lasso (<cit.>) where ℓ_1-penalty is employed to produce sparse coefficients. Many works has been proposed to efficiently compute the regularization paths for GLMs with lasso penalty. <cit.> introduced the predictor-corrector method to determine the entire path of lasso for GLMs. <cit.> proposed a fast algorithm exploiting cyclical coordinate descent to compute a regularization path for GLMs which can shorten computing time considerably. An important extension of the lasso method is adaptive lasso (<cit.>) which allows weights of predictors in the penalty term to be varied. <cit.> presented that the adaptive lasso enjoys the oracle properties in linear model and GLMs under mild condition. This means that it has variable selection consistency and the ability to correctly select the nonzero coefficients with probability converging to one. Much efforts have been devoted in investigating theoretical properties of adaptive lasso. <cit.> presented that the adaptive lasso has oracle property even when the number of covariates is much larger than the sample size in linear model. Other recent work including <cit.> and <cit.> also demonstrated the oracle property of the adaptive lasso in high-dimensional setting for GLMs. Another advantage of the adaptive lasso method is that it can be solved by the same efficient algorithm for solving lasso.

Alternatively, Spike-and-slab approaches to Bayesian variable selection has gained increasing attention recently. The spike part concentrates its mass at the values near zero and the slab component has a density spread over parameter spaces. There has been much work on hierarchical linear models with absolutely continuous spikes (i.e., mixture normal priors). Examples include <cit.>, <cit.>, <cit.>, <cit.>. MCMC stochastic search has been extensively used to fit the model with such Bayesian formulation in large-scale data analysis. However, it is computational expensive when the number of predictors p is large. Instead of using MCMC, <cit.> proposed an EMVS algorithm to fit a linear model with mixture normal prior. With conjugate spike-and-slab prior, the EMVS has closed form solutions for both E- and M- steps which can save computational time over MCMC. Although an absolutely continuous spike component has ability to determine the amount of shrinkage adaptively, regression coefficients are not sparse in the exact sense. To enforce some regression coefficients to be exactly zero,  a dirac spike using a point mass at zero should be considered instead. Nonetheless, EMVS algorithm is not in the closed form anymore. Moreover, extensions of the EMVS algorithm to GLMs are not straightforward.




<cit.> presented a contribution of empirical Bayes
thresholding (<cit.>) to select variables in a linear regression framework which can identify and estimate sparse predictors efficiently. A dirac spike-and-slab prior, a mixture prior of an atom of probability at zero and a heavy-tailed density, is put on each regression coefficient. An iterated conditional modes/medians (ICM/M) was also proposed for fast and easy-to-implement algorithm for empirical Bayes variable selection. Similar to an iterated conditional modes algorithm proposed by <cit.>, conditional modes are used to obtain hyperparameters and parameters other than regression coefficients. With a dirac spike-and-slab posterior for each regression coefficient, a conditional median is employed to enforce the variable selection. As demonstrated in <cit.>, empirical Bayes variable selection can also handle the case when the information about structural relationship among predictors is available through the Ising prior. The ICM/M algorithm can be easily implemented for such complicated prior.

The aim of this paper is to generalize empirical Bayes variable selection
 to more class beyond linear model. The general framework for extension of empirical Bayes variable selection to GLMs will be described. The details on the implementation of ICM/M algorithm for binary logistic model and Cox's proportional hazards model will also be discussed here. 

The rest of of this paper is organized as follows. The next section presents a Bayesian model formulation for high-dimensional GLMs. Section 3 introduces the ICM/M algorithm for GLMs. In Section 4, we introduce how to quantify variable importance and select variables based on false discovery rate. Section 5 and Section 6 focus on the implementation of ICM/M algorithm for binary logistic model and Cox's proportional hazards model respectively. Numerical results based on simulation studies are also be shown in Section 5 and Section 6. Application to real data sets are presented in Section 7 and 8. We then conclude this paper with a discussion in Section 9.



§ A BAYESIAN MODEL FORMULATION FOR HIGH-DIMENSIONAL GLMS

Consider a generalized linear model (GLM) where each response variable Y_i is assummed to be independent and has distribution belonging to exponential family taken the form

    f_Y_i(y_i|θ_i, ϕ) = exp{y_iθ_i-b(θ_i)/a(ϕ)+c(y_i, ϕ)}

where a(.),b(.), and c(.) are functions which vary according to distributions. θ=(θ_i, ..., θ_n) is known as the canonical parameter and ϕ is a dispersion parameter assummed to be known here. A link function g(.) is used to connect the linear predictors η_i = β_0+𝐗_iβ to the mean μ_i=E[Y_i]. That is, g(μ_i)=η_i=β_0+𝐗_iβ.

With an n × 1 matrix of responses 𝐘 and n × p matrix containing values of p predictors, the maximum likelihood estimators (MLEs) are obtained to estimate β in low-dimensional setting (n > p). An iteratively reweighted least squares (IRLS) algorithm is typically used to estimate such coefficients. Based on current parameter estimates (β̂_0, β̂), a quadratic approximation to the likelihood is employed to construct pseudodata 𝐙 = (Z_1,..,Z_n)^t and pseudovariances Σ=diag{σ_i^2} as follows:

    𝐙   = η̂+ (𝐘-μ̂)(dη/dμ|η=η̂)
    Σ   = diag{σ_i^2} = diag{ Var[Y_i] (dη/dμ|η=η̂)^2}

where η̂ = β̂_0+𝐗β̂. <cit.> showed that the update of β in IRLS is the result of a weighted least squares regression of 𝐙 on 𝐗 where the weighted matrix is Σ^-1. Indeed, the underlying distribution of pseudodata 𝐙 is approximated by N(β̂_0+𝐗β̂, Σ).

When n≪ p, a unique solution for the classical MLEs does not exist due to the fact that the design matrix 𝐗 is not full rank. Thus, the regression coefficients β cannot be updated using a weighted least squares in IRLS procedure. Furthermore, typical high-dimensional data analysis usually requires that only a subset of important variables are used in the modeling process. Variable selection to get sparse regression coefficients is a crucial task when analyze such massive data.

To introduce sparsity in the modeling process, an independent mixture of an atom of probability at zero and a distribution describing non-zero effect can be put on each of the regression coefficient β_j as prior in the Bayesian framework. With the approximated distribution of pseudodata 𝐙 and assuming that all other parameters except β_j are known, β_j has a sufficient statistic (X_j^tΣ^-1X_j)^-1(X_j^tΣ^-1Z̃_j) where Z̃_j = 𝐙 - η̂ + X_jβ_j and the approximated distribution of the sufficient statistics is N(β_j, (X_j^tΣ^-1X_j)^-1). Following <cit.>, the prior of β_j is in the form

    β_j | ω∼ (1-ω)δ_0(β_j) + ωγ(β_j).

Under the prior distribution in (<ref>), each β_j is zero with probability (1-ω) and β_j is drawn from the nonzero part of prior γ(β_j) with probability ω. Laplace distribution is employed for γ(β_j), that is,

    γ(β_j) = α(X_j^tΣ^-1X_j)^-1/2/2exp{ -α(X_j^tΣ^-1X_j)^1/2 |β_j|},

where α > 0 is a scale parameter. 

When the information of structural relationships among predictors is available, an indicator variable τ = (τ_1,...,τ_p)^t where τ_j = I_{β_j0} is introduced and the underlying relationships represented by an undirected graph is put on τ. Specifically, given τ_j, the prior of β_j is

    β_j | τ_j ∼ (1-τ_j) δ_0(β_j)+τ_j γ(β_j)

Again, the Laplace distribution is employed for γ(β_j). To model the relationship among τ under an undirected graph G=(V,E) comprising a set V of vertices and a set E of edges, the following Ising model is considered:

    P(τ) = 1/Z(a,b)exp{a∑_j τ_j + b ∑_<j,l>∈ Eτ_j τ_l}

where a and b are parameters and Z(a,b) is a normalizing constant, that is, 

    Z(a,b) = ∑_τ∈{0,1}^pexp{a∑_j τ_j + b ∑_<j,l>∈ Eτ_j τ_l}

When b>0, the interaction between nearest neighboring nodes is called ferromagnetic, i.e., neighboring τ_j and τ_l tend to have the same value. When b<0, the interaction is called anitiferromagnetic, i.e., neighboring τ_j and τ_l tend to have different values. When b=0, the prior gets back to independent and identical Bernoulli distribution. The value of a+b indicates the preferred value of each τ_j. In fact, τ_j tends to be one when a+b>0 and τ_j tends to be zero when a+b<0.



§ THE ITERATED CONDITIONAL MODES/MEDIANS ALGORITHM FOR GLMS

<cit.> presented an iteratated conditional modes/medians (ICM/M) algorithm for fast computation of empirical Bayes variable selection in a linear model. Data-driven optimal values for hyperparameters and auxiliary parameters are obtained as the modes of their full conditional distribution functions. Each regression coefficient is obtained as the median of its full conditional distribution function. With a spike and slab prior on β as in (<ref>) and  (<ref>), the regression coefficients obtained from conditional medians are sparse which therefore enforce variable selection and estimation simultaneously. The iterative procedure for updating regression coefficients and other parameters is carried out until convergence.

Adopting the idea of IRLS, we can extend the ICM/M algorithm to GLMs. Based on the approximated distribution of the current pseudo data, the ICM/M algorithm is applied to update all parameters. The procedure consists of the outer and inner loop. The outer loop is taken place to update pseudodata and pseudovariances based on current parameter estimates. The inner loop is where the ICM/M is employed to cycle through all parameters update.  



  
Proposition 1  With current values of β_0 and β, pseudodata and pseudovariances {𝐙, Σ} can be calculated and a sufficient statistic for β_j is (X_j^t Σ^-1 X_j)^-1(X_j^t Σ^-1Z̃_̃j̃) w.r.t. the approximated distribution of 𝐙 which is N(β_0+𝐗β, Σ). Then, the iterative conditional median of β_j in the ICM/M algorithm can be constructed as the posterior median of β_j in the following Bayesian analysis, 

    {[             (X_j^t Σ^-1 X_j)^-1/2(X_j^t Σ^-1Z̃_j) | β_j ∼ N( (X_j^t Σ^-1 X_j)^1/2, 1); β_j | ω∼ (1-ω)δ_0(β_j) +ωα (X_j^t Σ^-1 X_j)^1/2/2exp{-α (X_j^t Σ^-1 X_j)^1/2 |β_j| }. ].



With the independent prior as in (<ref>), the details of the ICM/M algorithm is demonstrated in Algorithm 1. 


  
Algorithm 1 


  * Initialize (β_0^(0), β^(0)) and set k=0.

  * Set k = k+1.

  * Update pseudodata and pseudovariances (𝐙^(k), Σ^(k)).

  * Update weight ω as the mode of its full conditional distribution function.

    ω^(k) = mode P(ω | 𝐘, 𝐗, 𝐙^(k), Σ^(k), β_0^(k), β^(k))


  * For j=1,...,p, update β_j as the posterior median. 

    β_j^(k) = median P(β_j | 𝐘, 𝐗, 𝐙^(k), Σ^(k), ω^(k), β_0^(k-1), β_1:(j-1)^(k), β_(j+1):p^(k-1))


  * Update β_0 as follows:

    β_0^(k) = Z̅_w - X̅_w β^(k)

where Z̅_w and X̅_w are the weighted means of 𝐙^(k) and 𝐗 respectively with weight for each observation being σ_i^-2.

  * Iterate between 2 - 6 until convergence.


To incorporate the information of structural relationships among predictors when it is available, the Ising prior (<ref>) is employed. Since the normalizing constant of the Ising model Z(a,b) is often computationally expensive, the values of hyperparameters a and b can be computed by maximizing the pseudo-likelihood instead of the original prior likelihood. Specifically, the pseudo-likelihood is as follows:

    L_p(a,b) = ∏_j=1^p P ( τ_j | {τ_l: <j,l> ∈ E}, a, b ) = ∏_j=1^p exp{τ_j (a+b∑_<j,l>∈ Eτ_l)}/1+exp{a+b∑_<j,l> ∈ Eτ_l }

Indeed, a and b are the logistic regression coefficients when the binary variable τ_j is regressed on ∑_<j,l> ∈ Eτ_l for j=1,...,p. Thus, the conditional median of β_j can be constructed on the basis of Proposition 2 and the details of the algorithm can be found in Algorithm 2.



  
Proposition 2 With current values of β_0 and β, pseudodata and pseudovariances {𝐙, Σ} can be calculated and a sufficient statistic for β_j is (X_j^t Σ^-1 X_j)^-1(X_j^t Σ^-1Z̃_̃j̃) w.r.t. the approximated distribution of 𝐙 which is N(β_0+𝐗β, Σ). Then, the iterative conditional median of β_j in the ICM/M algorithm can be constructed as the posterior median of β_j in the following Bayesian analysis, 

    {[                     (X_j^t Σ^-1 X_j)^-1/2(X_j^t Σ^-1Z̃_j) | β_j ∼ N( (X_j^t Σ^-1 X_j)^1/2, 1); β_j | ϖ_j ∼ (1-ϖ_j)δ_0(β_j) +ϖ_j α (X_j^t Σ^-1 X_j)^1/2/2exp{-α (X_j^t Σ^-1 X_j)^1/2 |β_j| }, ].

where the probability ϖ_j is specified as follows,

    ϖ_j^-1 = 1+exp{-a-b ∑_l:<j,l> ∈ Eτ_l}




  
Algorithm 2 


  * Initialize (β_0^(0), β^(0)) and set k=0.

  * Set k = k+1.

  * Update pseudodata and pseudovariances (𝐙^(k), Σ^(k)).

  * Update τ^(k) = (τ_1^(k), ..., τ_p^(k)) where τ_j^(k) = I_{β_j^(k-1) 0}.

  * Update hyperparameters (a,b) as the mode of its pseudo-likelihood function.

    (a^(k), b^(k)) = mode{∏_j=1^p P (τ_j^(k) | τ_l^(k) : <j,l> ∈ E ) }


  * For j=1,...,p, update β_j as the posterior median. 

    β_j^(k)   = median P(β_j | 𝐘, 𝐗, 𝐙^(k), Σ^(k), τ^(k), a^(k), b^(k)β_0^(k-1), β_1:(j-1)^(k), β_(j+1):p^(k-1)) 
       = median P(β_j | 𝐘, 𝐗, 𝐙^(k), Σ^(k), ϖ_j^(k), β_0^(k-1), β_1:(j-1)^(k), β_(j+1):p^(k-1))

where the probability ϖ_j^(k) =  [ 1+exp{-a^(k)-b^(k)∑_l:<j,l>∈ Eτ_l^(k)} ]^-1.

  * Update β_0 as follows:

    β_0^(k) = Z̅_w - X̅_w β^(k)

where Z̅_w and X̅_w are the weighted means of 𝐙^(k) and 𝐗 respectively with weight for each observation being σ_i^-2.

  * Iterate between 2 - 7 until convergence.




§ VARIABLE IMPORTANCE AND FALSE DISCOVERY RATE

Although the ICM/M algorithm with a spike-and-slab prior can enforce variable estimation and selection simultaneously, it is often of interest to quantify the importance of variables. <cit.> proposed a local posterior probability to evaluate the importance of variables. Specifically, given all other parameter estimates except for β_j, a local posterior probability of j-th predictor is defined as

    ζ_j    = P(β_j  0 | 𝐘, 𝐗, β̂_0, β̂_1:(j-1), β̂_(j+1):p, .) 
       = ∫ P(β_j | 𝐘, 𝐗, β̂_0, β̂_1:(j-1), β̂_(j+1):p, .).

However, with the Bayesian formulation for GLMs, the local posterior probability in (<ref>) is not in the closed form and it is natural to estimate it by using the pseudodata in the last iteration (𝐙, Σ) and its approximated distribution. That is, 

    ζ_j ≈∫ P(β_j | 𝐙, Σ, 𝐗, β̂_0, β̂_1:(j-1), β̂_(j+1):p, .).

Such probability in (<ref>) has a closed form and can be easily computed. 

With the local posterior probability ζ and true β, a true false discovery rate (FDR) given the data can be computed as

    FDR = ∑_j=1^p 1_{β_j=0} 1_{ζ_j>κ}/∑_j=1^p
    1_{ζ_j>κ},      0 ≤κ < 1

Following <cit.>, the expected FDR given the data in Bayesian scheme is defined as

    FDR = ∑_j=1^p (1-ζ_j) 1_{ζ_j>κ}/∑_j=1^p
    1_{ζ_j>κ}

We then can select a set of important predictors based on the local posterior probability and FDR. By controlling FDR at prespecified level, κ can be chosen and a set of important predictors lists all predictors having the local posterior probability greater than κ. It has been shown that the local posterior probability is a good indicator to quantify the importance of variables.



§ THE ICM/M ALGORITHM FOR BINARY LOGISTIC REGRESSION



 §.§ Impplementation details

When the response variable is binary taking values {0,1}, the logistic regression is commonly used. The logistic regression assumes Y_i is independently Bernoulli distributed with mean E[Y_i] = P(Y_i=1) = π(𝐗_i) and variance Var[Y_i] = π(𝐗_i)(1-π(𝐗_i)) where

    π(𝐗_i) = e^β_0+𝐗_iβ/1+e^β0+𝐗_iβ.

Equivalently, the logistic regression model can be written as

    log( π(𝐗_i)/1-π(𝐗_i)) = β_0+𝐗_iβ.


The pseudodata and pseudovariance for binary logistic regression can be constructed based on the current parameter estimates (β̂_0, β̂) as

    Z_i    = β̂_̂0̂+𝐗_iβ̂ + Y_i -π̂(𝐗_i)/π̂(𝐗_i)(1-π̂(𝐗_i))
    σ_i^2    = 1/π̂(𝐗_i)(1-π̂(𝐗_i))


where π̂(𝐗_i) = e^β̂_0+𝐗_iβ̂/1+e^β̂_0+𝐗_iβ̂ for i=1,...,n. 

Some cautions need to be considered to avoid any divergence issues. The details in the implementation of ICM/M algorithm are discussed here to prevent divergence problem:


  * When a probability is within the range of ϵ =10^-5 of 0 or 1, we set it to 0 and 1 respectively.

  * When a probability is close to 1. Numeric difficulty arises in calculating pseudodata and pseudovariances due to the term π̂(𝐗_i)(1-π̂(𝐗_i)) in the denominator. Hazard or tail functions are employed to improve numerical stability in IRLS (<cit.>). We define

    F(η_i) = e^η_i/1+e^η_i = 1/1+e^-η_i,


    h_+(η_i) = F(η_i)       and      h_-(η_i) = 1 - F(η_i).

The pseudodata and pseudovariance can be computed as

    Z_i    = {[ η̂_i + 1/h_+(η̂_i)      , if  Y_i = 1; η̂_i - 1/h_-(η̂_i)      , if  Y_i = 0 ] . 
    σ^2_i    = 1/h_-(η̂_i)h_+(η̂_i).

In addition, to avoid the overflows, when η̂_i < -30 and η̂_i>30, η̂_i is set to -30 and 30 respectively. Care should also be taken to select the appropriate form of the function in each tail. Specifically, the form involving e^-η should be used for positive value of η and the form involving e^η should be used for negative value of η.

  * The idea of active-set convergence is adopted to check convergence in each iteration. The algorithm stops when a complete cycle to update the coefficients in the inner loop does not change the active set of predictors - those with nonzero coefficients. Active set convergence is also mentioned in <cit.>, <cit.>, and <cit.>.




 §.§ Simulation studies

Simulation studies were conducted to evaluate the performance of our proposed logistic empirical Bayes variable selection via the ICM/M algorithm. Large p small n data sets were simulated from the model (<ref>). We compare the performance of our approach with other two popular methods: the regularized logistic regression with lasso (Lasso) and adaptive lasso (ALasso) penalties. A tuning parameter for both methods was chosen based on 10-fold cross-validation. Using lasso fits as initial values, the ICM/M algorithm was carried out. 

Three cases of prior on structured predictors were considered here. Case 1 assumed that all predictors are mutually independent. For Case 2 and Case 3, the information of structural relationships among predictors was assumed to be known and the ICM/M with the Ising prior was applied in these two cases. Average misclassification rates (MR) calculated on the test data sets, number of false positives (FP), number of false negatives (FN), model sizes (MS), ∥β-β̂∥_1 = ∑_j-1^p |β_j-β̂_j|, and ∥β-β̂∥_2^2=∑_j=1^p |β_j-β̂_j|^2 were calculated among 100 simulated data sets. 



  §.§.§ Case 1: Independent prior

The data for case 1 was generated with n=250 and p=1,000. The intercept term β_0 was set to zero and there are 10 non-zero regression coefficients including  β_1=...=β_5=10 and β_11=...=β_15=-5. The covariates were partitioned into 10 blocks, where each block containing 100 covariates were serially correlated at the same level ρ. The values of ρ were {0, 0.3, 0.5, 0.7, 0.9}.



From Table <ref>, misclassification rates among the three methods were comparable and they tended to get smaller when the correlation among predictors was high.

In order to assess the ability to select the correct variables, average number of false positives and false negatives are also reported in Table 1. Lasso produced large number of false positives across all values of ρ. Adaptive lasso had relatively lower false positive rates than lasso. ICM/M had the smallest number of false positive rates among the three methods. However, it had relatively higher number of false negatives. Furthermore, we also observed that the number of false positives got smaller and the number of false negatives got bigger for higher value of ρ. While the true model contained only 10 non-zero regression coefficients, the average model size from lasso was five to nine times larger than those in the true model. Adaptive lasso selected moderate numbers of non-zero regression coefficients. The average model sizes of ICM/M was closest to the true model. However, several true positives are missed especially when the correlation among predictors was high (ρ = 0.9). 

For estimation errors, ICM/M outperformed the other two methods in terms of both ∥β-β̂∥_1 and ∥β-β̂∥_2^2 except when ρ=0.9 that adaptive lasso was better than ICM/M.




In order to evaluate the importance of variables through the local posterior probability obtained from ICM/M, the true and estimated FDR were plot against κ in Figure <ref> . We first noticed that the estimated FDR was inflated from the true FDR. Furthermore, the FDR was well approximated when the correlation among covariates got higher. By controlling the FDR at level 0.05, the threshold of κ was chosen and all variables having local posterior probabilities suppassed the threshold were selected. Table <ref> shows the results on the value of κ, number of false positives, number of false negatives, and number of important variables for Case 1 simulation. We observed that the value of κ was decreasing for higher value of ρ causing lower number of false positives and higher number of false negatives. Overall, the results in Table <ref> are better than those in Table <ref>.  Thus, the local posterior probability and estimated FDR was a suitable tool to select important variables.






  §.§.§ Case 2: Linear chain prior

For case 2, we again fixed n=250 and p=1000. The intercept term β_0 was set to zero and the values of regression coefficients β depend on an indicator variable τ which follows a Markov linear chain model with transition matrix:

    T =    τ_j+1=0    τ_j+1=1 τ_j=0     0.99     0.01 τ_j=1     0.5     0.5

And τ_1 ∼π, where π = (0.5      0.5). The effect size for non-zero coefficients were drawn from Uniform[3,10]. The covariates 𝐗 were generated from AR(1) with different value of ρ in {0, 0.3, 0.5, 0.7, 0.9}.



As shown in Table <ref>, ICM/M had smallest misclassification rates when predictors were independent or had mild correlation among them while adaptive lasso had smallest misclassification rates for mild/high correlation among predictors. Again, we saw the same pattern as in Case 1 simulation that misclassification rates were lower when correlation among predictors was higher for all three methods. 

ICM/M had much smaller number of false positives comparing to lasso and adaptive lasso across values of ρ. Lasso and adaptive lasso generally selected large number of non-zero predictors as the true model in Case 2 simulation contains 21 non-zero coefficients. ICM/M produced relatively larger number of false negatives and smaller model sizes especially for high correlation among predictors.

When comparing estimation errors, ICM/M yielded the smallest ∥β-β̂∥_1 while adaptive lasso yielded the smallest ∥β-β̂∥_2^2 across all values of ρ. Thus, it depended on the circumstances to choose between ICM/M and adaptive lasso based on estimation errors. If we would like to avoid large errors, adaptive lasso is preferred since it penalized large deviation more heavily. But if we prefer a model that the estimates are generally close to the actual values but miss badly in some of the coefficients then ICM/M is more appropriate than adaptive lasso.



With the Ising prior, the estimated FDR curves obtained from ICM/M were underestimated making the chosen value of κ to be smaller than it supposed to be. When FDR was controlled at level 0.05, the results in Table <ref> had smaller number of false positives and slightly higher number of false negatives. To reduce the number of false negatives, we can increase the FDR level.





  §.§.§ Case 3: Pathway structure prior

 To assess the performance of our proposed method in pathway-based analysis, the data was simulated based on the genotype from publicly available Parkinson's disease (PD) dataset (dbGaP study accession number: phs000089.v3.p2). Parkinson metabolic pathway and other six pathways related to PD were obtained from KEGG database. For each genetic region in these pathways, SNPs resided in the genetic region were mapped to the gene based on their physical locations. With this procedure, most genes had at least one affiliated SNPs and there were few genes with more than thousand associated SNPs. To make the simulated data manageable, we selected all SNPs for genes having 5 or less than 5 associated SNPs. For those genes having more than 5 associated SNPs, only 5 SNPs representing each gene were randomly chosen. This process yielded p=1,152 SNPs representing 341 genes left for analysis. The phenotype 𝐘 was generated using logistic regression with these 1,152 SNPs as covariates. The regression coefficients for all 46 SNPs affiliated to 15 genes resided in Parkinson metabolic pathway were chosen to have non-zero effect size. The location of non-zero coefficients were fixed but the effect sizes were random and drawn from Uniform[1,10]. There were 1,741 individuals in the PD dataset which we randomly selected 871 individuals to be served as training data and 870 individuals as test data. The simulation was run for 100 times. 



The results from Case 3 simulation were consistent with simulation studies in Case 1 and Case 2. ICM/M outperformed the othere two methods in terms of number of false positives although it has higher number of false negatives. ICM/M also had lowest misclassificationr rate in this case. For estimation errors, ICM/M and adaptive lasso also had similar ∥β-β̂∥_1 but adaptive lasso yielded smaller ∥β-β̂∥_2^2.



As shown in Figure <ref>, FDR curve was well approximated in this case. In addition, the results in Table <ref> were similar to the results in Table <ref> when FDR was controlled at level 0.05.



To sum up the simulation section for binary logistic model, it is worth noticing that ICM/M with lasso fits as initial values improved significantly over lasso in terms of number of false positives. When correlation among predictors was very high, ICM/M might miss some important variables. However, it was interesting to see that ICM/M was able to reduce model size from lasso dramatically while increasing predictive ability at the same time. Furthermore, instead of relying on only regression coefficients from ICM/M, the local posterior probability along with the plot of estimated FDR can be used as a tool to select important variable.



§ THE ICM/M ALGORIHM FOR COX'S PROPORTIONAL HAZARDS MODEL



 §.§ Implementation details

the classical survival analysis framework, the observed data consists of (Y_i, 𝐗_i, δ_i) where Y_i is the ovserved time for subject i, 𝐗_i is a p-dimensional vector of covariates, and δ_i is 1 if Y_i is the actual survival time and 0 for right-censoring. The Cox's regression (<cit.>) models survival times through hazard function

    h(Y_i|𝐗_i) = h_0(Y_i)e^𝐗_iβ

where h(Y_i|𝐗_i) is the hazard for subject i and h_0(Y_i) is the baseline hazard function (i.e., the hazard function when all covariates 𝐗_i are 0). The Cox's proportional hazard model assumes that covariates are time-independent and the baseline hazard h_0(Y_i) is an unspecified function.


Let t_1<t_2<...<t_m be the increasing list of distinct failure time and define cumulative baseline hazard at time t_j for subject i to be

    H_0(Y_i) = ∑_j:t_j≤ Y_iΔ H_0(t_j),

where Δ H_0(t_j) is the increment of the cumulative baseline hazard at time t_j.


Following <cit.> and <cit.>, the regression coefficients can be found by maximizing the extended likelihood taking the form

    L(β) = ∏_i=1^n  [ {Δ H_0(Y_i)e^𝐗_iβ}^δ_iexp{-H_0(Y_i)e^𝐗_iβ} ].

Given value of β, the increment of the cumulative baseline hazard that maximizes the likelihood is

    ΔĤ_0(t_j) = ∑_i:Y_i=t_jδ_i/∑_i: Y_i≥ t_j e^𝐗_iβ.

Hence,

    Ĥ_0(Y_i) = ∑_j:t_j≤ Y_iΔĤ_0(t_j).

The explicit form of the pseudodata and pseudovariace can be computed based on current coefficients β̂ as 

    Z_i     = 𝐗_i β̂ + 1/Ĥ_0(Y_i) e^𝐗_i β ( δ_i - Ĥ_0(Y_i) e^𝐗_i β̂ )
    σ_i^2     = 1/Ĥ_0(Y_i) e^𝐗_iβ̂.

The approximated distribution of pseudodata 𝐙 is N(𝐗β̂, Σ) where Σ = diag{σ_i^2}. The ICM/M algorithm can now be used to cycle through parameters update and the active-set convergence is employed as stopping criterion.



 §.§ Simulation studies

In this section, three cases in simulation were conducted to show the performance of the ICM/M algorithm for Cox's proportional hazards model. Here we compare the results with lasso and adaptive lasso for Cox's model in large p small n settings. In all three scenarios, we fixed n=250 and p=1,000. Survival times were simulated from a Cox model with the baseline hazard function of a Weibull distribution with a shape parameter ν=10 and a scale parameter λ=1. The censoring times were generated randomly to achieve censoring rate of 50%.

Case 4 simulation assumed independence among regression coefficients and ICM/M with independent prior in (<ref>) was carried out. Case 5 and Case 6 assumed that regression coefficients followed some prior distribution representing the relationship among them. The ICM/M with the Ising prior was applied for Case 5 and Case 6 simulations. Here we report average number of false positives (FP), number of false negatives (FN), model sizes (MS), ∥β-β̂∥_1 = ∑_j-1^p |β_j-β̂_j|, and ∥β-β̂∥_2^2=∑_j=1^p |β_j-β̂_j|^2 calculated among 100 simulated data sets for each case.



  §.§.§ Case 4: Independent prior

The covariates were simulated the same way as in Case 1 with n=250 and p=1,000. Among 1,000 predictors, the failure times were determined by a linear combination of 20 non-zero coefficients: β_1 = ...=β_10=5 and β_101=...=β_110=2.  



Table <ref> shows that ICM/M performed the best regarding the ability to select true important variables across all levels of correlation among covariates. The average numbers of both false positives and negatives were close to 0. Moreover, the average model size was close to the size of true model which was 20. Although lasso also had small numbers of false negatives, it turned out that lasso tended to select much larger number of non-zero coefficients which was more than 100 variables. This caused large numbers of false positives. Adaptive lasso selected moderate size of predictors but its numbers of false positives and negatives were relatively large comparing to those produced by ICM/M.

In a comparison of ∥β-β̂∥_1, ICM/M surpassed the other two methods for ρ=0, 0.3, and 0.5. For ρ=0.7 and 0.9, adaptive lasso performed slightly better than ICM/M. When comparing estimation errors in terms of ∥β-β̂∥_2^2, adaptive lasso performed well when correlation among covariates was high while the performance of lasso and ICM/M were better than adaptive lasso for small value of ρ. Furthermore, we also observed that lasso and ICM/M yielded similar values of ∥β-β̂∥_2^2.



Figure <ref> shows that the estimated and true FDR curves were very close to each other. In addition, the results in Table <ref> and Table <ref> were consistent. The average number of false positives and negatives were all close to zero making the ICM/M was a competitive method for high-dimensional variable selection in Cox's model.





  §.§.§ Case 5: Linear chain prior

Covariates with n=250 and p=1,000 were generated the same way as in Case 2. The location of non-zero coefficients follows a Markov linear chain with the same transition matrix as in Case 2. The effect sizes of those non-zero coefficients were drawn from Uniform[0.5,5]. 



Ability to select variables correctly were assessed in 3 criteria: number of false positives, false negatives, and model size. As you can see in Table 5, ICM/M was the best method to select correct variables here due to small numbers of false positives and false negatives across all value of ρ.  Average model sizes were also close to the true model size containing 19 non-zero coefficients. We also noticed that ICM/M gave slightly higher number of false positives when ρ increased. In contrast, lasso and adaptive lasso yielded lower number of false positives when ρ increased. Although lasso and adaptive lasso performed reasonable well with regard to number of false negatives, both methods always selected large number of unimportant variables resulting in high false positive rates. Besides the variable selection ability, ICM/M overall outperformed the other two methods regarding to the two criteria of estimation errors including ∥β-β̂∥_1 and ∥β-β̂∥_2^2. Numbers of important variables in Table <ref> when controlling FDR at level 0.05 were slightly higher than the model size in Table <ref> but still very close to the true number of important variables.  








  §.§.§ Case 6: Pathway structure prior

For Case 6, gene expression data within an assumed network were simulated. The network consisted of 10 disjoint pathways. Each of which contained 100 genes resulting in p=1,000 in total. Ten regulated genes were assumed in each pathway. The gene expression values were generated from a standard normal distribution. For those regulated genes in the same pathway, the expression values were generated from normal distribution with a correlation of ρ=0.7 among those 10 regulated genes. Eighteen regulated genes from three pathways were chosen to have non-zero coefficients that were drawn from Uniform[0.5, 5].



As you can see in Table <ref>, ICM/M outperformed the other two methods in terms of the ability to select true relevant genes. Both average numbers of false positives and negatives were close to 0 with the average model size being 17.86 which is close to 18. Lasso and adaptive lasso selected much more variables into the model resulting in large number of false positives. We also observed that adaptive lasso had smallest avarage values for ∥β-β̂∥_1 and ∥β-β̂∥_2^2. 



Similar to other simulation studies for Cox's model, the true FDR curve was well approximated by the estimated FDR curve. Moreover, the results based on the regression coefficients from ICM/M method in Table <ref> were similar to variable selection procedure at FDR level = 0.05 in Table <ref>.



To conclude the simulation results for Cox's proportional hazards model, ICM/M gave the best variable selection performance. Lasso and adaptive lasso tended to select much larger number of variables into the model yielding large value of false positives though they might have lower number of false negatives than ICM/M. This suggested the substantial improvement of ICM/M with lasso coefficients as initial values over lasso. In addition, we also noticed that the results based on ICM/M coefficients were comparable to variable selection based on the local posterior probability when FDR is being controlled at level 0.05.



§ ILLUSTRATION: PARKINSON'S DISEASE GENOME-WIDE ASSOCIATION STUDY

As an application in genomics, a GWAS was carried out to identify a set of susceptible genetic markers to Parkinson's disease (PD). This PD dataset is publicly available (dbGaP study accession number: phs000089.v3.p2). The dataset consisted of 940 PD patients and 801 controls resulting in 1,741 samples in total. All individuals were Americans with European ancestry. The genotyping assays were derived from Illumina Infinium HumanHap300 and HumanHap500 SNP chips. There were 310,860 SNPs in common between these two arrays which we focused in this analysis. The following GWAS preprocessing criteria were applied to the data: (i) missingness per individual < 10%, (ii) missingness per marker < 10%, (iii) minor allele frequency (MAF) ≥ 5%, and (iv) Hardy-Weinberg test at significant level = 0.001. In addition, K-nearest neighbor (KNN) method was employed to impute missing genotype.
Parkinson metabolic pathway and other six pathways related to PD were acquired from the KEGG database. There were 26,101 SNPs representing 341 genes in the seven pathways. To make the computation more manageable, the univariate tests were applied and only 2,664 SNPs having p-values surpassed the significant threshold at 0.1.

Assume known pathway information among SNPs from the KEGG database, ICM/M with Ising prior was applied to fit a high-dimensional binary logistic regression. Again, lasso fits were used as initial values of coefficients for ICM/M algorithm. ICM/M was then compared with lasso and adaptive lasso. 

From Table <ref>, Lasso performed relatively poor due to large misclassification rate and model size. Adaptive lasso had the smallest misclassification rate. Similar to results in simulation studies, ICM/M tended to select smaller number of variables into the model comparing to lasso and adaptive lasso. There was also an improvement of ICM/M in terms of misclassification rate over lasso.



Figure <ref> demonstrates the number of SNPs identified by the three methods. There were 24 SNPs identified by all three methods. The details of these overlapping SNPs are shown in Table <ref>. We noticed that the signs of β̂ for each SNPs were the same for three methods. Overall, ICM/M tended to produce larger effect sizes for regression coefficients than the other two methods. This might due to the thresholding of ICM/M to screen out the small regression coefficients.





For ICM/M, FDR obtained from the local posterior probability was estimated. When controlling FDR at level 0.05, κ was chosen to be 0.86 yielding 114 important variables. If FDR level was increased to 0.10, the corresponding κ = 0.64 yielding 157 important variables. 

Due to large number of identified SNPs, we further investigated our analysis. For comparison purposes, we only focused on the results of ICM/M based on the fitted model for the rest of this section although one might consider using FDR and the local posterior probability to select important variables instead.

Since the correlation among SNPs resided in the same genetic region tended to be high, it was possible that different variable selection methods select different sets of SNPs representing the same gene. We therefore obtained genetic markers by grouping SNPs affiliated to the same gene. The genetic marker for gene g was defined as 𝐗̃_g = 𝐗_g β̂_g where 𝐗_g was a design matrix consisting of all SNPs resided in gene g and β̂_g was a vector of corresponding estimated coefficients obtained from variable selection method when considered all SNPs in original data as predictors. With the construction of genetic markers, it was easy to check which genes revealed genetic risk to PD. For each method, the logistic regression model with disease status as a response and genetic markers as covariates was fitted to obtain the p-values. The numbers of identified SNPs within the genetic regions are reported here. It was possible that some SNPs were affiliated to more than one gene. Table <ref> provides the results for genes identified by ICM/M (Full table can be found in Supplementary Table A).



There were 28, 35, and 24 genes identified by lasso, adaptive lasso, and ICM/M respectively (See Supplementary Table A). Among 13 genes detected for all three methods, 4 of them including MAPT, PARK2, PIK3CD, and SNCA had been well studied and known as PD susceptible genes in a number of literatures (see e.g., <cit.>, <cit.>, <cit.>, <cit.>, and <cit.>). There were another 6 genes that had been reported in Harmonizome database (see <cit.>) to have association to PD based on GWAS and other genetic association datasets from the GWASdb SNP-Disease Associations dataset. These 6 genes included BIRC6, COX7B2, PRKAR2B, TNPO3, and TRIP12. Among these 6 genes, lasso, adaptive lasso, and ICM/M can identify 5, 4, and 4 genes respectively. In addition, UCHL1 had also been reported to be involved with PD risk in recent literatures (<cit.> and <cit.>). Here, only ICM/M and adaptive lasso are able to detect UCHL1 gene.

We observed that there was a relationship between p-values and number of identified SNPs resided in the genetic region. Indeed, the genes having larger number of identified SNPs tended to have smaller p-values. This phenomenon was not unusual due to the fact that genes with larger number of identified SNPs were accounted for more variations to the response. Interestingly, we noticed that the p-value for ICM/M method tended to be smaller than the other two methods despite the number of identified SNPs resided in a gene. 
		


§ ILLUSTRATION: LUNG ADENOCARCINOMA MICROARRAY ANALYSIS

We illustrate our proposed methodology for Cox's model to the lung cancer data from <cit.>. The microarray data consisting of p=7,129 genes were from n=86 patients with primary lung adenocarcinomas. The survival time of patients in early-stage lung adenocarinomas were recorded and the censoring rate in this data set is 27.9%. 

Three pathways related to lung adenocarcinoma were obtained from KEGG database. Genes in the data set were then mapped to genes in these three pathways to generate an undirectedg graph. This resulted in a gene network consisting of 55 nodes and 109 edges. The rest of the genes were either isolated or not in one of these three pathways. Using lasso coefficients as initial values, the ICM/M algorithm with pathway structure prior was applied to the data set. For comparison purposes, the results from lasso and adaptive lasso were also reported.

Lasso , adaptive lasso, and ICM/M selected 24, 23, and 12 genes respectively. Figure <ref> shows the number of overlapping genes that were identified by the three methods. The numbers of selected genes by lasso and adaptive lasso were similar and there were 21 overlapping genes between these two methods. On the other hands, ICM/M selected less genes than the other two methods. There were 9 genes identified by all three methods. Among these 9 genes, 4 of them including BAG1, FXYD3, HPIP, and SLC2A1 were among the top 100 genes related to survival of patients reported in Beer et. al. (2002) (see Table <ref>). Furthermore, we found that another 3 genes were reported to have some relation with lung cancer in literatures. This included INSL4 (<cit.>), PRKACB (<cit.>), and SPRR1B(<cit.>).





In addition to 9 genes in Table <ref>, ICM/M also identified 3 more genes including NME2, STMN1, and ZFP36 (see Table <ref>). In fact, NME2, identified by both Lasso and ICM/M, was among the top 100 genes related to survival of patients reported in Beer et. al. (2002). For STMN1 and ZFP36, only ICM/M was able to identify these 2 genes. Interestingly, both genes were related to lung cancer in recent literatures: STMN1 (<cit.>) and ZFP36 (<cit.>). 



For ICM/M method, estimated FDR were calculated and κ was chosen as 0.81 which corresponded to FDR level 0.05. This yielded 6 genes having local posterior probabilities surpass κ=0.81. If one would like to increase the number of selected genes, FDR level can be increased to any desired level. For example, by controlling FDR level at 0.10, the κ was chosen to be 0.66 yielding 8 important genes. To select all 12 genes, the corresponding FDR level was 0.15 in this case. 



§ DISCUSSION

<cit.> proposed an empirical Bayes methods in selecting massive variables for linear model. A apike-and-slab prior is employed to introduce sparsity in regression coefficient and the Ising prior is used to capture structural relationship among predictors. The ICM/M algorithm is also introduced to cycle through each coordinate of the parameters in the model. This paper intends to extend Empirical Bayes variable selection for GLMs. We here focus on two important models arised in analyzing genomic data including the binary logistic and Cox's proportional hazards models. A main challenge of such an extension has emerged as how to generalize the framework from linear to nonlinear models. Here we borrow the idea of IRLS like in fitting classical GLMs to the ICM/M algorithm. This simple-to-implement algorithm achieves fast computation even in high-dimensional data analysis. The ICM/M algorithm for both linear and nonlinear models is implemented in the  R package.

The lasso is the most widely used method in fitting high-dimensional GLMs. Its popularity stems from its computational feasibility and availability in various packages. Our results show that ICM/M with lasso fits as initial values improves significantly over lasso in terms of number of false positives for both logistic and Cox's proportional hazards model. It can reduce the model size from lasso substantially by thresholding property of the ICM/M algorithm. For binary logistic model, we observed that the predictive ability of our proposed method was also better than lasso. For Cox's model, ICM/M outperformed lasso and adaptive lasso in terms of the ability to select correct variables into the model. ICM/M is also compuational feasible even for large p.

In biomedical field, most researchers would like to integrate expert knowledge such as gene regulartory networks in data analysis. While lasso and adaptive lasso do not utilize any prior information on strutural information among predictors, an important feature of ICM/M is that it can also incorporate more complicated prior easily. Here we demonstrated how to integrate a Markov random field representing relationships among predictors through an Ising model for binary logistic and Cox's models. Simulation studies and real data show a considerably improvement of ICM/M from lasso.

Another advantage of our proposed methodology over lasso and adaptive lasso is that our method does not only provide regression coefficients but also a local posterior probability. A local posterior probability for each predictor is proposed to quantify the importance of variable in linear model (<cit.>). For Bayesian variable selction for GLMs, the local posterior probability is estimated using the pseudodata in the last iteration of IRLS. FDR in Bayesian scheme can also be computed and employed to select the final set of important variables. Therefore, our method provides flexibility to select varibles while controlling FDR at a desired level rather than relying on only regression coefficients.

Finally, we would like to note that our method seems to work better with survival data than binary data in terms of false negative rate. This might be due to the fact that binary data is less informative than continuous data even though the censoring rate is quite high (50%).



§ ACKNOWLEDGEMENT








§ REFERENCES






§ SUPPLEMENTARY TABLE


 l | c c c | c c c 
		*Table A: Results from PD data analysis: details on identified genes (Bolded genes were identified by all methods).  

		2*Gene (Location)     3|c|Number of identified SNPs     3|cp-value 

		     Lasso     ALasso     ICM/M     Lasso     ALasso     ICM/M 

		*Table A: Results from PD data analysis: details on identified genes (Bolded genes were identified by all methods) (continued). 

		2*Gene (Location)     3|c|Number of identified SNPs     3|cp-value 

		     Lasso     ALasso     ICM/M     Lasso     ALasso     ICM/M 

		7rContinue on next page
		ATP5F1 (1p13.2)     1     0     1     0.0088     1      1.01e-6 

		NDUFS2 (1q23)     0     0     1     1     1     0.0010 

		PIK3CD (1p36.2)     1     1     1     5.93e-8     0.0037     1.14e-8 

		CDC20 (1p34.1)     1     0     0     0.0372     1     1 

		AKT3 (1q44)     0     1     0     1     0.0113     1 

		LCLAT1 (2p23.1)     1     0     0     0.0150     1     1 

		BIRC6 (2p22.3)     2     0     1     5.23e-8     1     2.13e-9 

		TRIP12 (2q36.3)     1     1     1     0.0352     0.0638     4.51e-5 

		NDUFB4 (3q13.33)     159     141     62     1.16e-70     3.2e-58     8.69e-66 
	HGD (3q13.33)     159     141     62     1.16e-70     3.2e-58     8.69e-66 

		PRKAR2A (3p21.3-p21.2)     0     0     1     1     1     0.9480 

		IL1RAP (3q28)     2     0    1     0.0021     1     3.17e-5 

		COX7B2 (4p12)     1     1     0     0.153     0.0091     1 

		UCHL1 (4p14)     0     1     1     1     0.0317     0.0146 

		SNCA (4q21)     1     1     1     0.0043     0.0029     0.0020 

		ATP6V1G2 (6p21.3)     48     38     17     2.41e-28     5.99e-24     1.26e-38 

		TNF (6p21.3)     48     38     17     2.41e-28     5.99e-24     1.26e-38 

		PARK2 (6q25.2-q27)     5     5     3     3.29e-8     2.42e-5     1.34e-17 

		DDC (7p12.2)     0     2     0     1     0.8540     1 

		UBE2D4 (7p13)     1     1     0     0.0786     0.0049     1 

		PRKAR2B (7q22)     1     0     1     6.00e-4     1     7.01e-07 

		TNPO3 (7q32.1)     0     1     0     1     0.0950     1 

		ATP6V0A4 (7q34)     2     1     1     0.0098     0.0934     0.0002 

		TCEB1 (8q21.11)     0     1     0     1     0.0149     1 

		TNFRSF10B (8p22-p21)     1     0     1     0.0205     1     0.0183 

		TYRP1 (9p23)     1     1     0     0.1880     0.0062     1 

		CUL2 (10p11.21)     1     1     1     0.0339     0.0034     0.0005 

		ATP5C1 (10p15.1)     0     1     0     1     0.1770     1 

		FAS (10q24.1)     0     0     1     1     1     0.0044 

		LHPP (10q26.13)     0     1     0     1     0.4230     1 

		PRMT3 (11p15.1)     1     2     1     3.97e-5     0.0001     2.58e-7 

		PSMD13 (11p15.5)     0     1     0     1     0.0068     1 

		TYR (11q14-q21)     0     1     0     1     0.1320     1 

		NDUFA12 (12q22)     1     0     0     0.3050     1     1 

		UBE2N (12q22)     0     1     0     1     0.0128     1 

		UBE3A (15q11.2)     0     1     0     1     0.0214     1 

		UBE2Q2 (15q24.2)     0     1     0     1     0.0205     1 

		DET1 (15q25.3)     1     2     0     0.1450     0.0279      1 

		NDUFB10 (16p13.3)     0     1     0     1     0.0270      1 

		MAPT (17q21.1)     107     51     72     7.11e-52     2.78e-28     6.17e-66 

		ANAPC11 (17q25.3)     107     51     72     7.11e-52     2.78e-28     6.17e-66 

		NEDD4L (18q21)     0     1     0     1     0.9180     1 

		ECH1 (19q13.1)     0     1     0     1     0.0252     1 

		AKT2 (19q13.1-q13.2)     0     1     0     1     0.6110     1 

		CBLC (19q13.2)     0     1     0     1     0.0813     1 

		NDUFA11 (19p13.3)     1     0     1     0.0118     1     8.15e-5 

		ATP4A (19q13.1)     1     1     0     0.0682     0.0973     1 

		UBE2G2 (21q22.3)     2     0     0     0.3360     1     1 

		ATP5O (21q22.11)     0     0     1     1     1     7.98e-9 

	



