Penalized regression and risk prediction in genome-wide association studies

2.1. Data

We use the Crohn's disease and bipolar disorder case and control data provided by the WTCCC. The WTCCC has collected genotype data of about 500 000 SNPs for approximately 2000 samples for each of seven diseases, such as type 1 diabetes, hypertension, bipolar disorder and Crohn's disease, and 3000 controls [16]. For simulations, we use the genotype data of 28 501 SNPs on chromosome 10 for Crohn's disease cases and controls. For quality control purposes, per WTCCC recommendations, we remove some samples and retain 1748 Crohn's disease samples and 2938 control samples; we also exclude some SNPs as recommended. Next, we eliminate the SNPs with a minor allele frequency (MAF) less than 5%. Furthermore, to mimic practical situations while maintaining a reasonable size for repeated simulations, we test each SNP separately by a chi-squared test for its association with Crohn's disease, and remove those with p-values larger than 0.1. At the end, we have about 2300 SNPs left and use them throughout our simulations.

2.2. Model

Let Y _i = 0 or 1 be a binary disease indicator for subject i = 1,…,n, and X_ij subject i's minor allele number (0,1, or 2) for SNP j = 1,…,m. Our aim is to build a model to successfully estimate subject i's risk of disease, P(Y _i = 1|x_i), based on his or her SNP data x_i = (X_i1,…,X_im)^T. As in standard practice for binary outcomes, we use a logistic regression model:

(1)

where β₀ and β_k are unknown regression coefficients to be estimated; p ≤ m indicates any user-specified subset of the SNPs.

In unpenalized logistic regression with MLE, β₀ and β = (β₁,…,β_p)^T are estimated by maximizing the log-likelihood:

(2)

The MLE is asymptotically unbiased with fixed p as n → ∞, but it may not be for a large p. One possible remedy is to introduce regularization or penalization on regression coefficients. The use of certain penalties, such as LASSO [1] and SCAD [18], shrinks many regression coefficient estimates to be 0, effectively selecting a subset of SNPs to be used for prediction. Penalized logistic regression provides coefficient estimates for β₀ and β by maximizing a penalized log-likelihood [21]:

(3)

where λ ≥ 0 is a tuning parameter controlling the extent of penalization imposed by penalty P(β). LASSO regression uses

(4)

which is convex and computationally convenient. However, LASSO estimates are biased and may not be consistent. To avoid these issues, Fan and Li [18] proposed using the SCAD penalty P(β,λ) replacing λP(β):

(5)

for a = 3.7. While maintaining the capability of variable selection, the SCAD penalty does not introduce biased estimates for some larger coefficients. The TLP adaptively determines which larger coefficients will not be penalized by introducing a separate thresholding parameter τ > 0 [19]:

(6)

If a coefficient β_k > τ, it will not be further penalized. The TLP approaches the L₀-loss as τ → 0⁺. A penalized method without the capability of variable selection is ridge regression [2] with penalty

(7)

Certain true models might be best estimated using a hybrid penalty that simultaneously performs variable selection and continuous shrinkage [20]. In these settings elastic net penalized regression may be more suitable. Elastic net penalized regression has been shown to produce a sparse model with good prediction accuracy, possibly superior to LASSO, while simultaneously promoting the grouping of strongly correlated predictors [20]. Its penalty structure is a weighted combination of the LASSO and ridge penalties controlled by a user-specified mixing parameter α, which is restricted to [0, 1]. The naive elastic net penalty [20] is

(8)

where α is selected to match the desired balance of variable selection and coefficient shrinkage. Zou and Hastie [20] suggested that further gains may be possible from using a rescaled version of the elastic net penalty. However, Friedman et al. [21] used the naive version of the penalty in the R package glmnet they developed to perform elastic net penalized regression. Results presented here follow this convention and are not rescaled. For sparse true models (i.e. with few nonzero β_k's) with a large number of candidate predictors (i.e. a large p), variable selection is often beneficial. However, for non-sparse models with many small nonzero |β_k|'s, variable selection will be difficult and may not result in good performance. On the other hand, as the ridge penalty has the grouping function [20], ridge regression performs like model averaging. It is known that neither model selection nor model averaging can dominate the other, and each performs better under different situations [22, 23]. In the current context, especially with non-sparse true models, it is not clear how LASSO, SCAD, and TLP compare to the ridge penalty for risk prediction, or if the elastic net penalty is superior, which is one of our aims.

3.1. Simulation Setups

We use the real SNP data of the WTCCC control cohort to generate disease probabilities, π_i = P(Y _i = 1) . First, we randomly select p₁ causal SNPs (i.e. with corresponding β_k≠0). The true correlations for any two SNPs range from −0.8371 to 1 and approximately fit a symmetric unimodal distribution centered at 0. Table 1 provides summary statistics for all pairwise correlations for example sets of size p = 5,10,50,100,500,1000 randomly selected SNPs. Table 1 demonstrates how the true models with various numbers of p SNPs contain a diverse range of minor, moderate or strong correlations among the SNPs.

We use p₁ = 10 for two sparse models, one with strong effects (i.e. large |β_k|'s) and the other with only moderate effects (i.e. smaller |β_k|'s); we also use p₁ = 300 and p₁ = 900 for two non-sparse models. Second, we set β₀ = log(0.05/0.95) to emulate diseases with low prevalence, and follow Wray et al. [24] to create odds ratios (ORs, OR_k = exp(β_k)) of having disease for the p₁ causal SNPs. Specifically, we set OR_k = 1 + ϵ(OR₀ − 1) with ϵ randomly generated from a standard exponential distribution Exp(1) and OR₀ being the mean OR, which is 2.75 and 1.415 for the two sparse models and 1.17 and 1.125 for the two non-sparse models respectively. We also randomly choose the sign of each β_k to be positive or negative to reflect both risk and protective causal SNPs. Third, the disease probability π_i for each subject i = 1,…,2938 in the WTCCC control cohort, is generated according to logistic regression model (1) with only chosen causal SNPs.

Finally, we use each π_i sequentially to generate disease status Y _i ∼ Bin(π_i) ; this step is repeated until we have n = 2000 cases and n = 2000 controls (while the other cases or controls are ignored) for each simulated dataset. One hundred datasets were generated under each of the four true models.

For each simulated dataset a randomly selected half of both the cases and controls is used as training set for building regression models, while the remaining half is the test set used for unbiased assessment of performance. The performance of each method is evaluated in two distinct settings. In the first setting we rank all SNPs by the p-values of their univariate association with disease. Starting with a few of the most significant SNPs, we fit and refit the logistic model for each method, sequentially adding more and more top ranked SNPs into the model (1) to be fit. The structure of this scenario informs when the inclusion of increasingly less significant SNPs improves or deteriorates the performance. Gail [4] measured the impact of only seven SNPs on classification of one disease, breast cancer, finding a very minor effect. Although they were not directly studying prediction, Yang et al. [25] identified one trait, height, whose heritability could be explained better with models that considered many nonsignificant SNPs. Our first modeling scenario generalizes this previous work to measure the impact of including more and more SNPs (by design including less significant SNPs) on a spectrum of models with less and less true sparsity. Thus, the results can inform about underlying genetic architectures for which penalized regression can use additional SNPs to improve risk classification. The results presented in the following section for the unpenalized regression are from the usual MLE, while those for LASSO, SCAD and ridge use the tuning parameter λ selected via tenfold cross-validation (CV) to have the smallest prediction error for any given number of candidate SNPs.

As exhibited in Eqs. (6) and (8), the elastic net penalty depends on an additional parameter, α, and the TLP penalty requires specification of τ. Elastic net estimates are generated for each of a sequence of penalties defined by a uniformly spaced sequence of values for the mixing parameter, α. The elastic net regression models are fit starting with α = 0, corresponding to ridge regression, and then with α increased by units of 0.10 until α = 1, corresponding to LASSO regression. For the TLP, we apply a range of τ values chosen to yield a series of models with minor to major coefficient shrinkage. To save computing time for tuning parameter selection for the simulated datasets, we use an independent tuning dataset of an equal size generated exactly like the training and test datasets. The idea is similar to the CV except that we only need to fit a model once with the training data, then use the tuning data to calculate the prediction error and thus select λ and τ.

The second setting is designed to compare the performance of the methods with a large number of the candidate SNPs. In penalized regression the regularization parameter λ is systematically varied to generate a solution path of the regression coefficients, from which we identify a global maximum of some performance measurement to represent the best ever performance of the corresponding method.

For each method, the estimated β₀ and β_k from a training set are applied to the corresponding test set to obtain risk estimates, \documentclass{article}\usepackage{amsmath}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{amsfonts}\pagestyle{empty}\begin{document}$\hat{\pi}_{i}$ \end{document}. The correlation of the \documentclass{article}\usepackage{amsmath}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{amsfonts}\pagestyle{empty}\begin{document}$\hat{\pi}_{i}$ \end{document} and the true π_i for the test samples is computed and used to compare the predictive performance of the methods.

This metric has been used in risk and outcome prediction for GWAS data [24, 26]. In addition, we also utilize the AUC for test samples to assess the discriminatory capabilities of the regression methods. The AUC is the gold standard metric that has been most consistently used in the GWAS literature. The use of AUC also permits direct comparison to previous related work. R package glmnet was used to fit the LASSO, ridge and elastic net penalized regression models. SCAD models were fit using the R package ncvreg. TLP models for the simulated datasets were fit using Feature Grouping and Selection Over and Undirected Graph (FGSG) software implemented in Matlab [25] while those for the real data were fit using our own implemented R function. Computational time necessitated using the FGSG software, which was much faster in fitting penalized linear regression models. It is known that linear regression models perform well for binary traits with GWAS data [27]. We also compared the results from penalized logistic regression models fitted by the R function with those from linear models by the FGSG software for the first ten simulated datasets; their differences were within 0.031 in the correlation metric and within 0.01 in the AUC metric.

3.2 Main Results

We first investigate the effect of using an increasing number of top SNPs for risk prediction. Figure 1(a) presents the correlation between true risk and predicted risk, \documentclass{article}\usepackage{amsmath}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{amsfonts}\pagestyle{empty}\begin{document}$Corr(\pi_i, \hat{\pi}_i)$ \end{document}. For each of the four true models, \documentclass{article}\usepackage{amsmath}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{amsfonts}\pagestyle{empty}\begin{document}$Corr(\pi_i, \hat{\pi}_i)$ \end{document} for each method is plotted as a gray curve against the number of the top SNPs used in the candidate model (before penalization) for each of the 100 simulated datasets, and the mean correlation curve over all 100 simulations is plotted as a dark red curve. The elastic net and TLP results are for the data-tuned values of α and τ respectively. In addition, vertical lines mark the number of the SNPs that would meet a Bonferroni adjusted genome-wide significance level at 0.05 when evaluated individually using a chi-squared test. Examination of the curves beyond the vertical lines reveals situations in which better estimates of the disease risk can be obtained by considering more SNPs, including those failing to meet the genome-wide significance level. The horizontal lines mark the correlations obtained from the MLE of the true model (with exactly all the causal SNPs).

In the sparse model with strong effect sizes, all penalized methods predict risk nearly as well or better than the unpenalized method, as shown in Fig. 1(b) where only the maximum correlation across various numbers of top SNPs from each simulated dataset is plotted. For the sparse model with weaker effects and both non-sparse models, all penalized methods by far surpass the MLEs, even ones based on the true models. Among the penalized methods, LASSO, SCAD, TLP, and elastic net outperform ridge regression for sparse models, but the trend reverses for non-sparse models for all but the elastic net. As the number of causal SNPs increases or strength of effect decreases, the relative performance of the elastic net and TLP penalties improves. In fact, the elastic net outperforms all methods for the p₁ = 300 case, which is to be expected as it is a model balanced between extreme sparse and non-sparse models. The best performing elastic net models are at least as good in the non-sparse p₁ = 900 case as those of ridge regression, the best overall performer of the non-mixture penalties. Table 2 provides the mean values of the maximum performance metrics of each regression method for the datasets. The table allows quick comparisons of the various methods in all modeling scenarios. These results reinforce the importance of using a suitable penalty for a given problem, depending on whether the model sparsity assumption holds.

To quantify the impact of including more SNPs, we first examine the performance for the sparse models. The LASSO and SCAD, methods with a variable selection feature, are able to maintain near optimal performance even when the number of candidate SNPs far exceeds that of the true model. Further, the elastic net appears to improve on the LASSO. In contrast, both unpenalized and ridge regressions have their prediction accuracy worsened markedly with the inclusion of more SNPs. For non-sparse models containing many SNPs failing to meet the genome-wide significance level, LASSO, SCAD, and TLP are again able to deal with a large number of SNPs for better risk estimation than the MLE. TLP uses the additional SNPs noticeably better than LASSO and SCAD when the true number of causal SNPs grows. Ridge regression is able to surpass these three penalization methods. In all four models the elastic net performs comparably to the best of the other regressions. This is likely because of its being a hybrid of the sparse and non-sparse regression methods, and our method examined a range of α's corresponding to a range of models from those strongly favoring LASSO to those strongly favoring ridge regression. However, it is noteworthy that the elastic net was not bounded by the performances of LASSO and ridge regressions.

Next, the discriminatory abilities of the methods are assessed because correct classification of disease status is key to personalized medicine. The literature for the clinical application of disease assessments universally reported AUCs as the standard for comparing disease classification methods. Therefore, the current study will assess classification using this metric to enable comparisons to previous work. Figure 2 shows the classification performance of the methods in terms of their AUCs. The main conclusions remain the same: SCAD, closely followed by LASSO , elastic net, and TLP, is the winner for the two sparse models, while elastic net and ridge regression beat other methods for the non-sparse model with p₁ = 900. However, for the non-sparse model with p₁ = 300, ridge regression performs worse than all other penalization methods. Elastic net performs best, followed by LASSO, TLP, and then SCAD. Overall, elastic net is either the top performer or close to the top for all true models, and every type of penalized regression always beats MLE.

The results presented in Fig. 2 demonstrate the value of penalized regression in disease risk estimation and classification, especially in utilizing the information in less significant SNPs that may often go unused. A natural question is whether we can eliminate the need to rank SNPs marginally and examine all SNPs simultaneously. The below simulation results address this question. All the penalized methods start with a full model containing all available SNPs; by varying the tuning parameter λ monotonically, various models are fitted and their performance is assessed. Figure 3(a) provides curves for the correlations between true and predicted risk at any given value of λ for four of the penalized methods: LASSO, SCAD, ridge, and elastic net. Elastic net results for only models with α = 0.5 are shown. Because one value of τ that provides a single intuitive interpretation across all four true models does not exist, TLP results as a regularization path in terms of λ would have limited comparability to the results from the other models in this setting and are not presented here. As before, the result for each simulation is represented by a gray curve, and the mean curve across all simulations is plotted as a dark red curve. For comparison, the horizontal lines mark the correlations obtained from maximum-likelihood estimation using exactly the true causal SNPs. To facilitate plotting, for each penalized method, the value of λ is scaled by its maximum so that it falls inside the interval [0,1].

As before, SCAD, LASSO, and elastic net with α = 0.5 outperform ridge regression for sparse models, while for both non-sparse models ridge regression is the best when judged by their optimal performance shown in Fig. 3(b). Interestingly, LASSO outperforms SCAD in all situations, suggesting the robustness of LASSO to a large number of input variables. The performance of the elastic net with α = 0.5 is between that of LASSO and ridge in all cases as expected. This elastic net's results are closer to the better of ridge and LASSO in all four models; however, the degree to which the best method outperforms the balanced elastic net (α = 0.5) varies by true model. This provides strong evidence that matching the sparsity of the penalty to the model sparsity improves classification. Comparing with earlier results, we can conclude that simultaneous use of too many SNPs will deteriorate the performance of any penalized method, suggesting possible gain in performance by a preliminary screening of a large number of variables. Similar conclusions hold if AUC is used to measure the classification performance of the methods (Fig. 4); however, LASSO, followed closely by the elastic net (α = 0.5), is the overall winner, in particular it beats ridge regression even for the non-sparse model with p₁ = 300, indicating the necessity of variable selection for large p.

3.3. Other Results

Two of the penalties, SCAD and TLP, are non-convex. Thus, there may be multiple local maxima with respect to their corresponding penalized log-likelihood functions, leading to possibly different estimates with different starting values. To examine this issue the authors refit some SCAD and TLP models for the first 20 datasets. The refit models considered the top ranked 1000 SNPs at a few fixed λ values (and a fixed τ = 0.1 for TLP). Eight different sets of initial regression coefficient values were used as the starting values for SCAD and TLP: the estimated coefficient values with the true model, a vector of all zeros, and the coefficients estimated by LASSO at each of the six λ values: 0.01, 0.1, 1, 2.5, 5, and 10. This was performed for the true models with 10 (strong) and 300 causal SNPs to represent one sparse true model and one non-sparse true model. The R package SIS was used to fit the SCAD models as it allowed user-specified initial value sets. FGSG software was used to fit TLP models as before.

Figure 5 presents the findings. Each curve represents the average AUC at a given λ over the 20 datasets for each set of the starting values (with the solid one for the first set). The primary finding is that for the λ generating the best AUC given a set of initial coefficients, all eight sets yield comparable AUCs. Results for many of the SCAD models could not be obtained when λ exceeded 0.1 owing to numerical problems in the R package; the partial curves are still provided. Many AUC values were the same or within 0.01 for the SCAD scenarios, thus, given the scale the curves appear to overlap in the plots. Not surprisingly the AUC is impacted by the starting values used to find regression coefficient estimates. Importantly, the impact appears to be small near tuning parameter values yielding the top performing SCAD or TLP models.

Below is a short summary on computing time needed to fit each type of penalized regression models. We calculated the average CPU time for one value or one set of tuning parameters for each penalized regression method with 1000 candidate SNPs for the true models with 10 (strong) and 300 causal SNPs. For the ten strong SNP scenario, SCAD used approximately 30 s to fit a model, TLP used 20 s, and the model fitting using glmnet ranged from 1.5 s for ridge regression to 7 s for LASSO. For the 300 SNP scenario, SCAD used approximately 44 s per model fitting, TLP used 20 s, and glmnet ranged from 1.2 s for ridge regression to 5 s with LASSO.

4.1. Crohn's Disease

Current research has identified about 80 SNPs associated with Crohn's disease. Figure 6 shows that approximately only the top 50 SNPs are needed to obtain the best risk prediction for all the methods; however, this includes more than just those SNPs meeting the significance level of 0.05/373191. Interestingly, although TLP was the overall winner and all five penalized methods are better than the unpenalized one, the performance difference among the methods is small.

Figure 7 presents the results of the four penalized methods starting with all 5000 SNPs included in a candidate model. With such a large number of candidate SNPs, while the number of the truly predictive SNPs may be small, the ridge penalty is largely outperformed by LASSO and SCAD that are capable of variable selection. The ridge regression is similarly outperformed by an elastic net penalty that shifts part of the weight from the ridge penalty to the LASSO penalty.

4.2. Bipolar Disorder

Bipolar disorder is a condition in which people go back and forth between mania, periods of a very good or irritable mood, and depression [28]. Figure 8 presents the AUC results as the number of candidate SNPs was increased. Unlike Crohn's disease, penalized regression does not always outperform MLE. Elastic net penalized regression and TLP perform best, although again the performance difference among the methods is small. As shown in Fig. 8(b), the inclusion of many SNPs failing to reach the genome-wide significance level does not diminish the discrimination strength of the penalized methods, and in fact ridge and elastic net regression and TLP better use these extra SNPs than both LASSO and SCAD to exceed or nearly exceed its performance achieved with only the few significant SNPs.

Next we include all SNPs in each penalized regression model and vary the tuning parameter λ (Fig. 9). Again it seems that, with a large candidate model containing a large number of predictors, ridge regression performs less well than the other three penalized methods, perhaps because of the former's inability for variable selection. LASSO and elastic net with α = 0.5 are the winners.