Least absolute shrinkage and selection operator type methods for the identification of serum biomarkers of overweight and obesity: simulation and application

Background The study of circulating biomarkers and their association with disease outcomes has become progressively complex due to advances in the measurement of these biomarkers through multiplex technologies. The Least Absolute Shrinkage and Selection Operator (LASSO) is a data analysis method that may be utilized for biomarker selection in these high dimensional data. However, it is unclear which LASSO-type method is preferable when considering data scenarios that may be present in serum biomarker research, such as high correlation between biomarkers, weak associations with the outcome, and sparse number of true signals. The goal of this study was to compare the LASSO to five LASSO-type methods given these scenarios. Methods A simulation study was performed to compare the LASSO, Adaptive LASSO, Elastic Net, Iterated LASSO, Bootstrap-Enhanced LASSO, and Weighted Fusion for the binary logistic regression model. The simulation study was designed to reflect the data structure of the population-based Tucson Epidemiological Study of Airway Obstructive Disease (TESAOD), specifically the sample size (N = 1000 for total population, 500 for sub-analyses), correlation of biomarkers (0.20, 0.50, 0.80), prevalence of overweight (40%) and obese (12%) outcomes, and the association of outcomes with standardized serum biomarker concentrations (log-odds ratio = 0.05–1.75). Each LASSO-type method was then applied to the TESAOD data of 306 overweight, 66 obese, and 463 normal-weight subjects with a panel of 86 serum biomarkers. Results Based on the simulation study, no method had an overall superior performance. The Weighted Fusion correctly identified more true signals, but incorrectly included more noise variables. The LASSO and Elastic Net correctly identified many true signals and excluded more noise variables. In the application study, biomarkers of overweight and obesity selected by all methods were Adiponectin, Apolipoprotein H, Calcitonin, CD14, Complement 3, C-reactive protein, Ferritin, Growth Hormone, Immunoglobulin M, Interleukin-18, Leptin, Monocyte Chemotactic Protein-1, Myoglobin, Sex Hormone Binding Globulin, Surfactant Protein D, and YKL-40. Conclusions For the data scenarios examined, choice of optimal LASSO-type method was data structure dependent and should be guided by the research objective. The LASSO-type methods identified biomarkers that have known associations with obesity and obesity related conditions.


Background
The study of circulating biomarkers and their association with disease outcomes has become progressively complex due to advances in biotechnologies available for the measurement of these biomarkers, including multiplex technologies. Although the availability of numerous biomarkers to study disease outcomes is highly promising, high dimensional biomarker data present statistical challenges. The Least Absolute Shrinkage and Selection Operator (LASSO) [1] is a popular high dimensional data analysis method that may be utilized for these biomarker data because it can simultaneously perform regularization and variable selection, which can improve both prediction accuracy and interpretation. This method, originally proposed for the linear regression model, minimizes the residual sum of squares, subject to the sum of the absolute value of the coefficients being less than a tuning parameter [1]. For the binary logistic regression model, the residual sum of squares is replaced by the negative log-likelihood [2]. If this tuning parameter is large, there is no effect on the estimated regression parameters. However, as the tuning parameter gets smaller, this may cause some coefficients to be shrunk towards zero or set to be zero. Still, there has been extensive and ongoing research towards the improvement of this method in obtaining a more sparse and consistent solution. Therefore, the key aim of this study was to evaluate five extensions that have been proposed to improve the sparsity and consistency of the original LASSO method. A simulation study was performed to compare variable selection properties of the original LASSO method with the Adaptive LASSO (AL), Elastic Net (EN), Iterated LASSO (IL), Bootstrap-Enhanced LASSO (BL), and Weighted Fusion (WF) for the binary logistic regression model.

LASSO
Details on penalized binomial logistic regression have been previously described [2]. Let x i = (x i1 , …, x ip ) ⊺ for i =1,…,N denote p-predictors for N observations. Assume that responses for the binary logistic regression model can take values G = 1,2. Then, where β 0 is the intercept and β = (β 1 , β 2 , …, β p ) ⊺ is a pvector of regression parameters. This implies The LASSO method then finds parameter values to minimize is the penalty function for the LASSO.

Adaptive LASSO
The LASSO method has shown to not always provide consistent variable selection. The LASSO penalizes all coefficients equally, even when the coefficients are large.
In contrast, the AL uses adaptive weights to penalize coefficients differently [3]. AL uses a weighted penalty, where w j ¼ 1 β j j j v ;β j is the maximum likelihood estimate and v > 0. The weighted penalty will allow variables with larger coefficients to receive smaller penalties and thus might provide a more consistent solution.

Elastic net
Zou and Hastie proposed the EN to address three issues related to the LASSO [4]. The first two issues relate to highly correlated variables in the n > p situation. For highly correlated variables, the LASSO tends to choose one variable and not the others. Also, predictive performance for ridge regression was empirically observed to be better than LASSO [4]. The third issue relates to the p > n situation in which the LASSO can at most select n variables. The penalty term for EN incorporates both the ridge penalty [5] λ ⋅ ∑ j = 1 p β j 2 and LASSO penalty λ ⋅ ∑ j = 1 p |β j |: where α is a value between 0 and 1. The EN penalty is equal to the LASSO penalty when α = 1 and to the ridge penalty when α = 0. The EN selects groups of correlated variables together, shares nice properties of both the LASSO and ridge regression, and can be considered for situations with p > n.

Iterated LASSO
The purpose of the IL is to consider the AL where the weights are based on LASSO estimates of the coefficients rather than maximum likelihood estimates [6]. An initial estimator is obtained by the IL to reduce the dimension of the model. The IL uses the LASSO to obtain an initial estimator and reduce the dimension of the model. The LASSO estimates are then used for the weighted penalty.

Bootstrap-enhanced LASSO
The bootstrap sample x* = (x 1 *, x 2 *, …, x n *) is obtained by randomly sampling the initial data points x 1 , x 2 , …, x n with replacement. The BL takes several bootstrapped replications of a sample and then considers the intersection of these estimates [7]. The motivation behind the BL is that if there existed several data sets from the same distribution, relevant variables would appear in all data sets and by running the LASSO for several bootstrapped replications would lead to a consistent model. As the BL may be too strict in intersecting models, it has been recommended to use a softened version of the BL (BL-S) [8]. In particular, a BL-S that considered variables that appeared in at least 90% of the bootstrap samples was shown to have better performance than the BL [8]. Similarly, for this study rather than considering a strict intersection, we a priori considered variables that appeared in at least 75% of the bootstrap samples (BL-75), consistent with other application studies [9,10].

Weighted fusion
The motivation behind WF is to improve the EN by utilizing additional information from the correlation structure. In particular, the WF utilizes information from correlated variables by using correlation driven weights to penalize for the pairwise differences of these coefficients [11]. The penalty term is defined as where tuning parameters λ 1 ≥ 0, λ 2 ≥ 0, and J is a correlation driven penalty function, x j is the sample correlation between predictors, γ > 0 is a tuning parameter, and s ij is the sign of ρ ij . The correlation driven weights encourage correlated variables to be considered together.

Current study
The simulation study utilized the data structure of the population-based Tucson Epidemiological Study of Airway Obstructive Disease (TESAOD). This real dataset presents data scenarios that are likely to be encountered in serum biomarker research. Examples of such scenarios may include high correlation between biomarkers, weak to moderate associations with the outcome, and sparse number of true signals. The current study seeks to compare the LASSO and LASSO-type methods for scenarios similar to the TESAOD data with the intent that findings may be useful for similar biomarker studies. In the application study, a comparison of methods including the application of each LASSO and LASSO-type method to the TESAOD data was performed for the identification of serum biomarkers of overweight and obesity.

Simulation study
A simulation study was performed to compare the performance of the LASSO, AL, EN, IL, BL-75, and WF methods for the binary logistic regression model. A total of five scenarios were considered for the overweight and obese outcomes as shown in Table 1. Biomarkers for the simulation study were generated from the multivariate normal distribution using the R package 'mnormt' (version 1.5-3) [12] and the outcomes were generated from the binomial distribution using R package 'stats' (version 3.1.0) [13]. Tuning parameters were estimated by using 10-fold cross validation with deviance loss. For EN, the alpha parameters were also chosen by cross validation  [11]. The data structure of TESAOD was utilized in the simulation study design. A total of 879 subjects had 86 measured serum biomarkers, thus the choice of N = 1000 with 100 biomarkers was made. Simulations for N = 500 were also completed as in the case of stratified analyses or a smaller study population. The first four scenarios considered five biomarkers with non-zero association with the outcome (or true signals) out of the 100 biomarkers, while the fifth scenario considered 20 biomarkers truly associated with the outcome. Correlation of serum biomarkers ranged from <0.01 to 0.70 for TESAOD, thus the correlation levels of 0.20, 0.50, and 0.80 were chosen. In the TESAOD dataset, 40% of subjects were overweight and 12% were obese. Thus in the simulation study the proportion of overweight was set at 40% and that of obesity at 12%. For each scenario a total of 1000 simulated datasets were generated. Overall, five scenarios, two outcomes, and two sample sizes were considered for a total of 20 simulation studies.
To evaluate the variable selection properties of each method, we first assessed if the solution provided by each method was sparse. Second, we evaluated how well each method identified the predictors with correct non-zero and zero coefficients. In addition, we also evaluated the area under the Receiver Operating Characteristic curve (AUC) on an independent validation data set to evaluate the discriminatory predictive performance of each method. The independent validation data sets were generated for each of the 1000 simulations and for each scenario with similar data structure as the training data set. Depending on the scenario, the validation data sets had sample size of either 1000 or 500 per validation data set.

Application study Study population
TESAOD is a population-based prospective cohort study initiated in 1972 in Tucson, Arizona [14]. At study initiation, 3805 participants from a stratified cluster of 1655 white Tucson households between the ages of 6 to 95 were enrolled. At study baseline, participants completed a standardized questionnaire, had height and weight measured by a study nurse, and had a sample of their blood collected. For the current study, serum biomarkers were measured on 879 non-Hispanic white subjects who were between the ages of 21 to 70 years at the 1972 baseline survey and for whom a baseline serum sample with sufficient volume was available. The University of Arizona Institutional Review Board approved the TESAOD study.
A total of 86 biomarkers were considered for the application study, of which 81 were measured by Myriad-RBM and five were measured locally at the ARC. Of the 87 biomarkers measured by Myriad-RBM, six were not considered. Factor VII, Insulin, and Prostate Specific Antigen, Free were dropped as we only had measurements for 322 of the 879 subjects. Biomarkers that had greater than 90% undetectable values (i.e., Interleukin-2, Interleukin-12p70, and Lymphotactin) were also not considered. Missing data were present for Alpha-2 Macroglobulin (3 subjects), CC16 (2 subjects), YKL-40 (4 subjects), Eotaxin (1 subject), and SPD (3 subjects). Missing values were replaced with the median value for the given biomarker. Biomarker measurements fell into one of the following categories: undetectable values, normal values, and high values ( Table 10 in the Appendix). Biomarkers that had less than 15% of undetectable values were considered as continuous. These biomarkers were evaluated for normality and were log transformed to obtain approximate normality when appropriate. Their concentrations were analyzed in statistical analyses as standardized values. Samples used to measure the biomarker concentrations were run using different batches in the laboratory. Variability in measurement may be introduced by not running samples in one batch. To adjust for any batch effects, the ComBat function as part of the R package 'sva' (version 3.10.0) was used [16,17]. Biomarkers that had between 15 and 50% undetectable values were categorized at their median value. Biomarkers that had between 50 and 90% undetectable values were categorized at their detection limit.

Statistical analysis
For the TESAOD study data with the panel of 86 biomarkers, each LASSO and LASSO-type method was applied to two separate binary logistic regression models, one comparing the 306 overweight subjects with the 463 normal-weight subjects, and the other comparing the 66 obese subjects with the 463 normal-weight subjects. Furthermore, we additionally applied each LASSO-type method to a binary logistic regression model comparing the 372 overweight or obese subjects to the 463 normal weight subjects. In addition, separate univariable binary logistic regression models were performed to verify that each biomarker was independently associated with each outcome. For the BL-75, we report the median values for the variables that appeared in at least 75% of the bootstrap samples.
The analyses of simulated and real data were performed using R versions 2.15.1 and 3.0.2 and Stata version 12.0 (Statacorp LP, College Station, TX, USA). The R packages glmnet (version 1.9-8) [2,18], lqa (version 1.0-3) [19], and ROCR (version 1.0-5) [20,21] were used. Table 2 shows sparsity of results for all methods in each of the five scenarios. The only method that did not provide the most sparse solution for either the overweight or obese outcome across all scenarios was the EN. When considering scenarios with weak associations of biomarkers with the outcome (scenarios 1 and 2) and scenarios with highly correlated variables (scenarios 2 and 4), the IL provided the most sparse solution for the overweight outcome and for the obese outcome, the LASSO, IL, BL-75, and WF methods provided the most sparse solution given a particular scenario. For scenario 5 that considers a moderate number of true signals with moderate correlation, the IL provided the most sparse solution for the overweight outcome and the BL-75 provided the most sparse solution for the obese outcome. Tables 3 and 4 show each method's ability to correctly identify the non-zero and zero coefficients for each of the five scenarios. Considering N = 1000 with 100 biomarkers (Table 3), for the overweight outcome, the WF outperformed the other methods in identifying the nonzero coefficients, but performed the worst in identifying the zero-coefficients. The IL identified more of the zerocoefficients correctly. For the obese outcome, one method did not clearly outperform the other methods. Results for simulations with N = 500 with 100 biomarkers are shown in Table 4. For the overweight outcome, the WF tended to choose more of the non-zero coefficients correctly and the IL tended to identify more of the zero-coefficients correctly. For the obese outcome, the WF outperformed the others at identifying the correct non-zero coefficients, but tended to incorrectly include more noise variables. Table 5 shows median and range of AUC results. Discriminatory predictive performance increased as the association with the outcome increased and also with larger sample size. In general, when considering scenarios with either weak associations of biomarkers with the outcome (scenarios 1 and 2) or scenarios with highly correlated variables (scenarios 2 and 4), the LASSO method showed good predictive ability over the other methods. The LASSO, EN, and WF performed comparably and outperformed the other methods. The BL-75 generally outperformed the AL. When considering scenario 5 with a moderate number of true signals with moderate correlation, the BL-75 showed the best predictive ability for the overweight outcome. In most scenarios, the IL performed poorly.

Application results
A total of 86 biomarkers were considered for the application study. We identified the biomarkers that were consistently chosen across the six LASSO and LASSO-type methods and had a significant univariable association with   Table 6). The combination of both outcomes was also considered and 16 biomarkers were identified, namely Adiponectin, ApoH, Calcitonin, sCD14, C3, CRP, Ferritin, GH, IgM, IL-18, Leptin, MCP-1, Myoglobin, SHBG, SPD, and YKL-40 (Table 7).
Results for all biomarkers, not limited to those chosen across all 6 LASSO and LASSO-type methods, can be found in Tables 11-13 in the Appendix.

Discussion
In the simulation study, we compared the variable selection properties of the LASSO and five LASSO-type methods for the binary logistic regression model and did find certain situations in which one method outperformed the others. In general, when we considered scenarios with weak associations of biomarkers with the outcome and scenarios with high correlation between biomarkers, the IL tended to provide the most sparse solution, but had poor discriminatory predictive performance. The WF tended to correctly identify more of the true signals, but also incorrectly included more noise variables. Still, we were not able to identify one method that had an overall superior performance over the others. In general, our simulation set-up considered much smaller effects compared to those studied in the original methodological papers that proposed the LASSO [1], AL [3], EN [4], IL [6], BL [7], and WF [11] methods. This more realistic setting for biomarker research may contribute to why we did not see clear improvements of one method over the other.
Similar to our results, when comparing AL results to the LASSO, and considering both large and small effect sizes, Zou found that there was not one single method that consistently outperformed the others [3]. They found that in low sample size scenarios the LASSO performed the best with a low signal to noise ratio (SNR), while the AL outperformed the LASSO when high SNR was present [3]. In general, the LASSO was able to identify more of the non-zero coefficients as compared to   Both the LASSO and AL methods have been shown to have good prediction accuracy [3], however the LASSO did show better discrimination than the AL in our simulation study. The EN has been shown to outperform the LASSO when collinearity is present. However, the EN typically chooses more variables than the LASSO, creating a less sparse solution [4]. For scenarios with high correlation (scenarios 2, 4, and 5), our results show that the EN correctly identified the true signals either comparably or better than the LASSO (Tables 3 and 4). As expected, the EN included more noise variables for all scenarios. Similarly, in our simulation study the LASSO provided a more sparse solution as compared to the EN ( Table 2). The IL has been shown to provide a more sparse solution than the LASSO [6]. We confirmed in our simulation study that the IL tended to provide the most sparse solution as compared to the others, but it also demonstrated poor predictive ability. With sparse data-generating models, the BL has been shown to outperform the LASSO [8]. However, the BL has also been shown to be too stringent [8]. In the present study, in order to minimize the exclusion of any potentially important biomarkers, we a priori considered the BL-75. We found that when considering a moderate number of true signals with moderate correlation, the BL-75 showed the best predictive ability. WF has been shown to outperform the LASSO and EN [11]. Similar to the EN, WF also tends to over select variables and thus creates a less sparse solution. In our simulation study, we found that the WF correctly identified more true signals, but incorrectly included more noise variables. Given these biomarker data, choice of optimal LASSOtype method was dependent on the characteristics of how the data were generated and should be guided by the research objectives. For objectives that aim to identify the maximal number of true signals, the WF was most optimal and the IL and AL the least. While identifying the greatest number of true signals, the WF also included more noise variables. The LASSO and EN performed well in the identification of many true signals and exclusion of more noise variables. For objectives that aim to maximize prediction, the LASSO, EN, and WF would also be optimal. While we found that no method had a clear overall advantage over the others, the IL was outperformed by the other methods in both variable selection and prediction. Additionally, while the BL showed the best predictive ability when considering a moderate number of true signals with moderate correlation, it was outperformed in variable selection. Given that the original LASSO method was not outperformed by the other methods, this method would be the most ideal method since it is the most direct and efficient method to implement. However, in general we recommend that the choice of optimal LASSO-type method should be guided by the underlying scientific question and by the research objectives. When applying the methods to the TESAOD dataset, we considered three different outcomes, namely overweight, obese, and overweight or obese. We had expected similar biomarkers for both conditions and possibly stronger associations with obesity. We chose to combine overweight and obesity in a separate analysis as this would provide the highest power for most biomarkers as well as provide a summary measure. Considering overweight and obesity in separate models would allow us to estimate biomarker association levels specific to each condition. Not all biomarkers were chosen consistently across all methods when considering the overweight and obese outcomes separately. ApoH and sCD14 were consistently chosen across all methods for the overweight outcome and across all methods except for the IL for the obese outcome. Similarly, Adiponectin, C3, Ferritin, IgM and Myoglobin were consistently chosen across all methods for the overweight outcome and across all methods for the obese outcome except for AL, IL, and BL-75. Results from the simulation study suggest that the IL might not always choose the true signal and that the LASSO, EN, and WF might be more likely to identify the true signal as compared to AL and BL-75. Other biomarkers such as Calcitonin and GH were consistently chosen for the overweight outcome, but never chosen for the obese outcome. In contrast, MCP-1 and vWF were never chosen for the overweight outcome, but always chosen for the obese outcome. These differences may in part be due to the methodology and different sample sizes in the overweight and obese categories, but they also might be due to biological differences between the two outcomes.
Overall, in addition to the well-known effects of Adiponectin, CRP, and Leptin, the LASSO methods identified multiple biomarkers that have been reported to be associated with obesity and/or obesity related conditions  and that were largely classified into either a group of hormones or hormone related proteins (Adiponectin, Calcitonin, GH, Leptin, SHBG), or into a group of positive acute phase reactants and other biomarkers of inflammation (C3, CRP, Ferritin, IL-18, MCP-1, vWF, and YKL-40). The LASSO-type methods shrink coefficients and set other coefficients to 0, thus producing biased estimates. Of note, Tables 6 and 7 show biased estimates and the magnitude of the estimation differences could be noted. The tuning parameter chosen by cross validation affects the amount of shrinkage and this tuning parameter may differ between LASSO-type methods.

Strengths and limitations
A strength of this study is that real-world data were used in the development of the simulation study parameters. In particular, the scenarios that were considered represent a more realistic setting that might be present in high dimensional serum biomarker research such as sparse number of true signals, weak to moderate association of the biomarkers with the outcomes, and high correlation between biomarkers. In the application study we confirmed potentially important biomarkers of overweight and obesity using results that were consistent across six LASSO and LASSO-type methods.
A limitation to this study is that although we found that in general the IL tended to provide the most sparse solution and the WF tended to correctly identify the most number of true signals, the choice of optimal LASSO method is data structure dependent and results from this study may not be generalizable to other biomarker studies. Furthermore, for the BL method, rather than consider a strict intersection, we considered the BL-75. We acknowledge that this may not fully optimize results and a different frequency threshold may outperform the BL-75. In addition, the primary goal of this research was to study the binary logistic regression model. Particularly, we were interested in how the LASSO-type models would compare when considering a common (40% overweight) or less common (12% obese) outcome. Considering both outcomes together as a single ordinal or as a multinomial response might have been more efficient than to consider them as separate binary responses and performance of LASSO-type approaches for such outcome measures will be a future research topic.
In this paper we focus on situations where the number of variables p is smaller than the sample size N. As shown in Fan and Fan [46], Fan and Liv [47], and Fan, Samworth, and Wu [48], when p > N, performance of the LASSO-type methods is inferior to that of certain two-stage methods, which first apply a screening procedure to reduce the number of important variables below the sample size, then apply methods like LASSO to the selected subset of variables. Comparison of such twostage methods for the p > N situation is beyond the scope of this paper and will be the topic of future investigation.

Conclusions
For the data scenarios examined, the LASSO method was overall not outperformed by any of the other methods. Choice of optimal LASSO-type method was dependent on characteristics of how the data were generated. In general, these characteristics are unknown and can only to some extent be estimated. Nevertheless, choice of optimal LASSO-type method should be guided by knowledge of the underlying scientific question and by the research objectives. The LASSO-type methods were able to identify biomarkers that were known to be associated with obesity and obesity related conditions, demonstrating the promise of such methods in future investigations.