Comparing Feature Selection Methods in Clinical Data Modeling: LASSO, Ridge, Elastic Net, SCAD, Boruta, Best Subset, and Stepwise Regression

Comparing Feature Selection Methods in Clinical Data Modeling: LASSO, Ridge, Elastic Net, SCAD, Boruta, Best Subset, and Stepwise Regression Yumeng Wang1, Iulia Cristina Tudor1* 1Corcept Therapeutics, Redwood city, CA 94065 *Corresponding author: [email protected] Abstract Accurate variable selection is essential for predictive modeling in clinical research, yet guidance on method choice under varying data conditions is limited. We systematically evaluated stepwise regression, best subset selection, Boruta, and four regularization-based methods (LASSO, Ridge, Elastic Net, SCAD) using simulated datasets and a prospective study of endogenous hypercortisolism in patients with difficult-to-control type 2 diabetes. Key metrics included predictive performance, parameter estimation, and feature selection accuracy. Stepwise and best subset methods performed well in low-dimensional, low-collinearity settings, offering interpretable and parsimonious models without hyperparameter tuning. In medium-to-high dimensional or highly correlated settings, regularized methods yielded more stable and accurate predictions, while Boruta identified a complementary set of features overlooked by linear models. These results provide practical guidance for choosing variable selection methods and suggest that combining complementary approaches can enhance both feature discovery and interpretability in clinical datasets. Key Words: LASSO, stepwise regression, simulation, hypercortisolism. 1. Introduction Feature selection is a critical step in clinical data modeling, yet there remains a lack of practical guidance on when to use specific methods. Each approach has its own strengths and limitations, making the choice highly context-dependent. Clinical datasets often present unique challenges: a large number of candidate predictors, many of which may be irrelevant or noisy; relatively small sample sizes due to the high cost of patient recruitment; and strong correlations among variables. These factors increase the risk of overfitting and poor generalizability if not properly addressed. Moreover, interpretability is essential in clinical research, which explains why linear regression remains widely used despite the availability of more complex modeling techniques. Taken together, these characteristics underscore the importance of robust feature selection strategies in building reliable and clinically meaningful models. A commonly used classical approach for variable selection is stepwise regression, which iteratively adds or removes predictors based on information criteria or hypothesis tests. While stepwise methods have been widely applied in clinical research for their simplicity and interpretability, they are known to suffer from several limitations in high-dimensional settings. These include model instability, and a propensity to overfit, particularly when the number of features (p) is large relative to the sample size (n) (Derksen & Keselman, 1992; Harrell, 2015). Another commonly used classical approach is best subset regression, which attempts to identify the combination of predictors that provides the best fit according to a criterion such as AIC, BIC, or adjusted R2. Unlike stepwise procedures, best subset regression conducts an exhaustive search across all possible predictor combinations, which guarantees that the identified model is globally optimal under the chosen metric. The main advantages of this approach are its interpretability and transparency. However, it quickly becomes computationally prohibitive as the number of predictors grows, and like stepwise regression, it is prone to instability and overfitting in small-sample, highdimensional settings (Miller, 2002). To address these challenges, regularization methods have emerged as powerful alternatives by imposing penalties on model complexity. The Least Absolute Shrinkage and Selection Operator (LASSO) (Tibshirani, 1996) applies an L1 penalty, encouraging sparsity by shrinking some coefficients exactly to zero and thereby performing variable selection. Ridge regression (Hoerl & Kennard, 1970), which uses an L2 penalty, shrinks coefficients toward zero but retains all predictors, offering improved stability when predictors are highly correlated, as correlated features tend to shrink together. The Elastic Net (Zou & Hastie, 2005) combines L1 and L2 penalties, providing a balance between sparsity and stability, and is particularly effective when groups of correlated features are present. In addition, nonconvex penalties such as the Smoothly Clipped Absolute Deviation (SCAD) (Fan & Li, 2001) have been proposed to overcome the estimation bias introduced by LASSO, especially for predictors with large true effects, while still producing sparse and interpretable models. Beyond regression-based methods, feature selection approaches leveraging tree-based models have gained popularity in clinical and biomedical research. Methods such as Boruta (Kursa & Rudnicki, 2010) build upon Random Forests to assess the importance of each variable in predicting the outcome. Boruta iteratively compares the importance of real features to that of permuted features, selecting only those that consistently show stronger predictive power. The advantage of tree-based approaches lies in their ability to naturally capture nonlinearities and interactions among predictors, without the need for strong modeling assumptions. However, they often sacrifice interpretability compared to linear models. Taken together, these diverse methods highlight the trade-offs between interpretability, stability, and predictive performance that are central to feature selection in clinical data modeling. Despite the wide array of available approaches, there remains limited practical guidance for clinical researchers on when to use which method. This study aims to fill that gap by systematically comparing stepwise regression, best subset regression, LASSO, Ridge, Elastic Net, SCAD, and Boruta, using both simulated data and real clinical datasets. 1.1 Stepwise Regression via F-Test-Based Feature Selection Variable selection in stepwise regression is based on null hypothesis testing, most commonly through the F-test. In forward selection, the procedure begins with an empty model (or only the intercept). At each step t, the improvement in model fit from adding a candidate predictor is assessed using the partial F-statistic. The predictor with the smallest p-value is added, provided it falls below a prespecified inclusion threshold α. A larger α increases the likelihood of including more predictors but also raises the false positive rate. In contrast, backward elimination starts with the full model and sequentially removes the predictor with the largest p-value, continuing until no variable exceeds the elimination threshold. The bidirectional (stepwise) procedure combines both approaches, allowing variables to be added or removed at each step, so that predictors previously included can later be excluded. The algorithm terminates once no further predictors meet the entry or removal criteria, yielding a final model determined by sequential hypothesis testing. 1.2 Best subset regression In best subset regression, the goal is to identify the model of size 𝒌𝒌∈{𝟏𝟏,⋯,𝒑𝒑} that minimizes a chosen selection criterion, such as AIC, BIC, Mallows’Cp, or adjusted R2. A naïve implementation requires fitting all 2p possible models, which is computationally infeasible when p is large. To address this, the “leaps and bounds” algorithm (Furnival & Wilson, 1974) was developed, which prunes the search space by discarding subsets that cannot possibly improve the criterion value. This approach substantially reduces the computational burden while still guaranteeing that the optimal subset is identified for each model size. Practical implementations are widely available, including the leaps package in R for linear regression, and the glmulti package, which generalizes best subset selection to generalized linear models such as logistic regression (Calcagno & de Mazancourt, 2010). Despite its advantages in terms of interpretability and optimality under the chosen criterion, best subset regression can remain computationally demanding in moderate-to-high dimensions and is sensitive to collinearity among predictors. 1.3 LASSO LASSO adds an ℓ1 penalty to the ordinary least squares (OLS) objective: 𝜷𝜷�=𝒂𝒂𝒂𝒂𝒂𝒂𝒎𝒎𝒎𝒎𝒎𝒎 𝜷𝜷�𝟏𝟏 𝟐𝟐𝒎𝒎��𝒚𝒚𝒎𝒎−𝒙𝒙𝒎𝒎𝑻𝑻𝜷𝜷�𝟐𝟐 𝒎𝒎 𝒎𝒎=𝟏𝟏 +𝝀𝝀��𝜷𝜷𝒋𝒋� 𝒑𝒑 𝒋𝒋=𝟏𝟏 � where λ ≥ 0 is a tuning parameter that controls the strength of regularization. The ℓ1 penalty encourages sparsity, shrinking some coefficients exactly to zero, thereby performing variable selection. As λ increases, more coefficients are shrunk to zero, yielding simpler models; when λ=0, the solution reduces to ordinary least squares. In practice, λ is typically chosen via K-fold crossvalidation, minimizing prediction error on held-out folds. Efficient algorithms make LASSO computationally feasible even when p ≫ n, which has contributed to its widespread use in clinical applications. 1.4 Ridge Ridge regression uses an ℓ2 penalty instead: 𝜷𝜷�=𝒂𝒂𝒂𝒂𝒂𝒂𝒎𝒎𝒎𝒎𝒎𝒎 𝜷𝜷�𝟏𝟏 𝟐𝟐𝒎𝒎��𝒚𝒚𝒎𝒎−𝒙𝒙𝒎𝒎𝑻𝑻𝜷𝜷�𝟐𝟐 𝒎𝒎 𝒎𝒎=𝟏𝟏 +𝝀𝝀�𝝆𝝆𝒋𝒋𝟐𝟐 𝒑𝒑 𝒋𝒋=𝟏𝟏 � Unlike LASSO, Ridge does not produce sparse models. All coefficients are shrunk toward zero but remain nonzero. This is beneficial when all variables contribute meaningfully to the outcome, such as in medical diagnosis models where biomarkers collectively improve disease prediction. However, it does not perform variable selection in the strict sense. 1.5 Elastic Net Elastic Net combines the strengths of LASSO and Ridge by incorporating both ℓ1 and ℓ2 penalties: 𝜷𝜷�=𝒂𝒂𝒂𝒂𝒂𝒂𝒎𝒎𝒎𝒎𝒎𝒎 𝜷𝜷�𝟏𝟏 𝟐𝟐𝒎𝒎��𝒚𝒚𝒎𝒎−𝒙𝒙𝒎𝒎𝑻𝑻𝜷𝜷�𝟐𝟐 𝒎𝒎 𝒎𝒎=𝟏𝟏 +𝝀𝝀[𝜶𝜶��𝜷𝜷𝒋𝒋� 𝒑𝒑 𝒋𝒋=𝟏𝟏 +(𝟏𝟏−𝜶𝜶)�𝝆𝝆𝒋𝒋𝟐𝟐 𝒑𝒑 𝒋𝒋=𝟏𝟏 ]� where α ∈ [0,1] controls the mixing between LASSO (α=1) and Ridge (α=0). The dual regularization provides a flexible compromise between sparsity and stability. Tuning involves selecting both λ and α, usually via a two-dimensional grid search with cross-validation. 1.6 SCAD The Smoothly Clipped Absolute Deviation (SCAD) penalty was designed to overcome the bias of LASSO for large coefficients while still yielding sparse solutions. The SCAD penalty is defined by a piecewise function that behaves like LASSO near zero (encouraging sparsity), but flattens out for larger coefficients, thereby providing unbiased estimates for large coefficients. The SCAD estimator solves: 𝜷𝜷�=𝒂𝒂𝒂𝒂𝒂𝒂𝒎𝒎𝒎𝒎𝒎𝒎 𝜷𝜷�𝟏𝟏 𝟐𝟐𝒎𝒎��𝒚𝒚𝒎𝒎−𝒙𝒙𝒎𝒎𝑻𝑻𝜷𝜷�𝟐𝟐 𝒎𝒎 𝒎𝒎=𝟏𝟏 +�𝒑𝒑𝝀𝝀(�𝜷𝜷𝒋𝒋�) 𝒑𝒑 𝒋𝒋=𝟏𝟏 � where 𝒑𝒑𝝀𝝀 (⋅) is the SCAD penalty function depending on λ and a tuning parameter a (commonly set to 3.7). Compared to LASSO and Ridge, SCAD’s nonconvexity makes the optimization more challenging than for convex penalties. 1.7 Boruta algorithm The Boruta algorithm (Kursa & Rudnicki, 2010) is a wrapper method built around Random Forests to perform all-relevant feature selection. Unlike LASSO, Ridge, or SCAD, which impose penalties in regression models, Boruta leverages the variable importance measure from Random Forests. At each iteration, a Random Forest model is fitted, and the importance score of each real feature is compared against the maximum importance among the permutated features (shadows). Features consistently outperforming the shadows are deemed “confirmed”, while those consistently underperforming are “rejected”. Features with intermediate performance are left “tentative”. The main advantage of Boruta is its ability to detect weak but relevant features, though at the cost of higher computational burden compared to single-model methods. To summarize, the methods discussed above differ in their underlying selection mechanisms, computational efficiency, and behavior under multicollinearity or high-dimensional settings. Table 1 provides a side-by-side comparison of these seven approaches, highlighting their key features, advantages, and limitations. Table 1: Comparison of key features of the selected methods Method Type Feature Selection Approach Multicollinearity Interpretability Scalability Stepwise Regression Linear (GLM) Greedy (Both directions) Poorly High Moderate Best Subset Linear (GLM) Exhaustive Search Poorly High Slow LASSO Linear (GLM) Shrinkage (L1 penalty) Good High Good Ridge Linear (GLM) Shrinkage (L2 penalty) Good High Good Elastic Net Linear (GLM) Shrinkage (L1 + L2 penalty) Good High Good SCAD Linear (GLM) Shrinkage (Nonconvex penalty) Good High Moderate Boruta Treebased Importance from Random Forest Generally robust Moderate Moderate In this paper, we compare stepwise regression, best subset, LASSO, ridge regression, elastic net, and Boruta across both simulated datasets and real-world clinical data. Our evaluation focuses on predictive performance, parameter estimates, and feature selection accuracy under varying levels of sample size, predictor dimensionality, signal-to-noise ratio, sparsity, and collinearity. The remainder of the article is organized as follows. In Section 2, we present a simulation study designed to systematically compare the performance of these methods across controlled scenarios. Section 3 applies the methods to an empirical study aimed at identifying risk factors for hypercortisolism, illustrating their practical differences in a real clinical context. We conclude in Section 4 with a discussion of key findings and practical considerations for applying these techniques in clinical research. 2. Simulation Study 2.1 Simulation Design To systematically assess model performance under diverse clinical data scenarios, we conducted a comprehensive simulation study based on a five-factor factorial design. The parameters were selected to mimic a wide range of conditions commonly encountered in clinical research. Specifically, we varied the following factors: 1. Sample size (N): Values of 100, 500, and 1000 were considered to represent small, medium, and large study sizes typical of clinical settings. 2. Number of candidate predictors (p): We examined predictor sets of size 15, 30, 50, 100, and 1000, spanning low-, moderate-, and high-dimensional settings. 3. Signal-to-noise ratio (SNR): let’s denote the model as: 𝒚𝒚𝟎𝟎 =𝒇𝒇(𝒙𝒙𝟎𝟎)+ 𝜺𝜺𝟎𝟎 where 𝒙𝒙𝟎𝟎 ∈ Rp is the independent variable and 𝒚𝒚𝟎𝟎 ∈ R is the response variable. Then, the SNR is defined as: 𝑺𝑺𝑺𝑺𝑺𝑺=𝒗𝒗𝒂𝒂𝒂𝒂(𝒇𝒇(𝒙𝒙𝟎𝟎)) 𝒗𝒗𝒂𝒂𝒂𝒂(𝜺𝜺𝟎𝟎) SNR values ranged from 0.05 to 4.5, corresponding to 2–80% of variance explained, with SNR = 1 representing 50% explained variance. 4. Sparsity level: Sparsity was defined as the proportion of truly nonzero coefficients in the underlying model. To approximate realistic conditions, five predictors were assigned coefficients of 1, while the remaining coefficients decayed exponentially toward zero: 𝜷𝜷𝟎𝟎𝒎𝒎=𝟎𝟎.𝟓𝟓𝒎𝒎−𝒔𝒔, for 𝒎𝒎=6, … ,p. A threshold of 1e–2 distinguished true coefficients from negligible ones. By varying the total number of predictors, this setup captured differing assumptions about the number of relevant features. 5. Collinearity level: Predictors were generated from a multivariate normal distribution with first order autoregressive correlation structures. Correlation values of 0.2, 0.5, and 0.8 were used to represent weak, moderate, and strong collinearity, respectively. Each simulation scenario was replicated 50 times to estimate the mean and standard deviation of each evaluation metric. For each replication, we generated both a training set and a test set of equal size under the same parameter configurations. The training set was used for model fitting, while the test set served to evaluate predictive accuracy, as described below. 2.2 Evaluation Metrics To evaluate and compare the performance of each regression method across the simulated scenarios, we considered three aspects: predictive performance, parameter estimates and variable selection accuracy. These metrics were chosen to reflect the utility of the fitted model for prediction, the estimated coefficient values, and its ability to recover the true underlying structure of the data. • Predictive performance was measured using the normalized root mean squared error (NRMSE) on an independent test set: 𝑁𝑁𝑁𝑁𝑁𝑁𝑁𝑁𝑁𝑁=�1 𝑛𝑛𝑡𝑡𝑡𝑡𝑡𝑡𝑡𝑡∗𝑠𝑠𝑠𝑠(𝑦𝑦𝑖𝑖)�(𝑦𝑦𝚤𝚤 �−𝑦𝑦𝑖𝑖)2 𝑛𝑛𝑡𝑡𝑡𝑡𝑡𝑡𝑡𝑡 𝑖𝑖=1 where 𝑦𝑦𝑖𝑖 are the observed responses, 𝑦𝑦𝚤𝚤 � is the predicted responses, 𝑠𝑠𝑠𝑠(𝑦𝑦𝑖𝑖) is the standard deviation of the outcome in the test set, and 𝑛𝑛𝑡𝑡𝑡𝑡𝑡𝑡𝑡𝑡 is the number of observations in the test set. Normalization by outcome variability allows comparability across scenarios. • Parameter estimates were assessed by relative test error (RTE) (Hastie, Tibshirani, & Tibshirani, 2017), which measures the expected test error relative to the Bayes error rate: 𝑁𝑁𝑅𝑅𝑁𝑁(𝛽𝛽󰆹) = 𝑬𝑬�𝑦𝑦0−𝑥𝑥0𝑇𝑇𝛽𝛽󰆹�2 𝜎𝜎2 =(𝛽𝛽󰆹−𝛽𝛽0)𝑇𝑇 𝛴𝛴(𝛽𝛽󰆹−𝛽𝛽0) + 𝜎𝜎2 𝜎𝜎2 A perfect score of 1 indicates that the model exactly recovers the true coefficients, whereas the null score is given by 𝛽𝛽0𝑇𝑇𝛴𝛴𝛽𝛽0+𝜎𝜎2 𝜎𝜎2= 𝑁𝑁𝑁𝑁𝑁𝑁+ 1. • Variable selection accuracy was evaluated by comparing the estimated support set (nonzero coefficients) against the ground truth. We computed the Matthews Correlation Coefficient (MCC) (Matthews, 1975), which provides a balanced summary of all four classification outcomes (TP, TN, FP, FN): 𝑁𝑁𝑀𝑀𝑀𝑀= 𝑅𝑅𝑇𝑇∗𝑅𝑅𝑁𝑁−𝐹𝐹𝑇𝑇∗𝐹𝐹𝑁𝑁 �(𝑅𝑅𝑇𝑇+𝐹𝐹𝑇𝑇)(𝑅𝑅𝑇𝑇+𝐹𝐹𝑁𝑁)(𝑅𝑅𝑁𝑁+𝐹𝐹𝑇𝑇)(𝑅𝑅𝑁𝑁+𝐹𝐹𝑁𝑁) TP (true positives) are correctly selected nonzero coefficients, TN (true negatives) are correctly excluded zero coefficients, FP (false positives) are zero coefficients incorrectly selected, and FN (false negatives) are nonzero coefficients missed. The MCC ranges from -1 (complete disagreement) to 1 (perfect agreement), with 0 indicating random performance. • Positive variable selection accuracy was quantified using the F1 score (van Rijsbergen, 1979), which balances sensitivity and precision in identifying true predictors: 𝐹𝐹1 = 2∗𝑇𝑇𝑃𝑃𝑃𝑃𝑃𝑃𝑃𝑃𝑠𝑠𝑃𝑃𝑃𝑃𝑛𝑛∗𝑁𝑁𝑃𝑃𝑃𝑃𝑅𝑅𝑅𝑅𝑅𝑅 𝑇𝑇𝑃𝑃𝑃𝑃𝑃𝑃𝑃𝑃𝑠𝑠𝑃𝑃𝑃𝑃𝑛𝑛∗𝑁𝑁𝑃𝑃𝑃𝑃𝑅𝑅𝑅𝑅𝑅𝑅 , where 𝑇𝑇𝑃𝑃𝑃𝑃𝑃𝑃𝑃𝑃𝑠𝑠𝑃𝑃𝑃𝑃𝑛𝑛=𝑇𝑇𝑇𝑇 𝑇𝑇𝑇𝑇+𝐹𝐹𝑇𝑇 and 𝑁𝑁𝑃𝑃𝑃𝑃𝑅𝑅𝑅𝑅𝑅𝑅=𝑇𝑇𝑇𝑇 𝑇𝑇𝑇𝑇+𝐹𝐹𝐹𝐹 A key limitation of F1 score is its bias toward larger models: adding more variables typically increases recall (reducing false negatives) and thus inflates F1, even if many irrelevant predictors are included. To address this, we applied an adjusted F1 score that directly penalizes model size: 𝐴𝐴𝑠𝑠𝐴𝐴𝐴𝐴𝑠𝑠𝐴𝐴𝑃𝑃𝑠𝑠 𝐹𝐹1 = 𝐹𝐹1∗(1 −𝑁𝑁𝐴𝐴𝑁𝑁𝑁𝑁𝑃𝑃𝑃𝑃 𝑃𝑃𝑜𝑜 𝑜𝑜𝑃𝑃𝑅𝑅𝐴𝐴𝐴𝐴𝑃𝑃𝑃𝑃𝑠𝑠 𝑠𝑠𝑃𝑃𝑅𝑅𝑃𝑃𝑃𝑃𝐴𝐴𝑃𝑃𝑠𝑠 𝑁𝑁𝐴𝐴𝑁𝑁𝑁𝑁𝑃𝑃𝑃𝑃 𝑃𝑃𝑜𝑜 𝑜𝑜𝑃𝑃𝑅𝑅𝐴𝐴𝐴𝐴𝑃𝑃𝑃𝑃𝑠𝑠 ) This adjustment reduces scores for unnecessarily large models, providing a fairer comparison across methods with different levels of sparsity. 2.3 Parameter tuning Stepwise regression was performed with a p-value cutoff of 0.10 for variable entry and 0.05 for retention. For the best subset regression model, the Bayesian Information Criterion (BIC) was used to determine the optimal model size. All cross-validation procedures were conducted using 10-fold cross-validation. The LASSO was tuned using cross-validation to select the penalty parameter λ. Two versions of the LASSO model were considered: one using λmin, the value that minimizes the cross-validated prediction error, and another using λ1se, the largest value of λ within one standard error of the minimum. These are referred to as LASSO (min) and LASSO (1se), respectively, in the following sections. Ridge regression was tuned in a similar manner, with λ selected according to the λ1se rule. Elastic net regression was fit with equal weighting of the ℓ1 and ℓ2 penalties (α=0.5), and λ tuning was performed following the λ1se rule. For SCAD, the concavity parameter a was set to the commonly used default value of 3.7, and λ was chosen as λmin. In the Boruta procedure, only features confirmed as important were retained as selected predictors. The final random forest model was then built and tuned through a grid search over five candidate values of the number of predictors considered at each split, with cross-validation used to select the optimal configuration. For stepwise regression and best subset selection, computational burden increased substantially as the feature space grew. To address this, we implemented predictor prefiltering based on the marginal correlation between each predictor and the response. Specifically, when the number of predictors was ≥ 1000, stepwise regression was restricted to the top 100 features. For best subset selection, when the number of predictors was ≥ 100, the analysis was restricted to the top 50 features. 2.4 Simulation Study Results 2.4.1 Predictive Performance Predictive performance was evaluated using normalized RMSE, with results summarized in Figure 1. Each bar represents a method, where scenarios with p < 100 are considered low-dimensional and those with SNR < 0.5 are considered low-signal settings. Dots indicate model size, mapped to the secondary y-axis on the right. Model size for Ridge regression is not shown, since all coefficients remain nonzero and its model size is equal to the full feature space. As expected, performance improved with higher SNR across all methods, yielding smaller RMSE values. In low-dimensional settings (small number of predictors), differences between stepwise regression, best subset, and regularization-based methods were minimal, and became negligible as sample size increased. Neither SNR level nor predictor correlation substantially altered this pattern. For LASSO, we compared two tuning strategies: λmin which minimizes cross-validation error and results in larger models, and λ1se, which applies stronger regularization and produces more parsimonious models. In most cases, λmin achieved better predictive accuracy, while λ1se remained smaller in size, offering a trade-off between accuracy and parsimony. In high-dimensional settings, the picture was markedly different. Regularization methods substantially outperformed stepwise and best subset approaches, particularly when sample size was small. These performance gaps diminished as sample size increased. Ridge regression, however, exhibited notable overfitting when SNR was high, since it shrinks coefficients but does not eliminate irrelevant predictors, leading to inflated prediction error. Tree-based models displayed distinct behavior. Random forests performed comparably well under low SNR, but under high SNR they underperformed relative to stepwise regression. This suggests that their inherent bias, resulting from bagging and random splits, can hinder performance when the underlying signal is purely linear. Interestingly, predictor collinearity improved the predictive accuracy of random forests: when multiple correlated predictors carried the true signal, the probability of capturing it at a split increased, thereby enhancing model performance. Overall, regularization-based methods showed the most stable predictive performance across simulation scenarios, particularly in high-dimensional settings. While SNR had little impact on the relative ranking of methods, collinearity modestly benefited tree-based approaches. These results suggest that regularization methods are generally preferable for clinical datasets with many predictors, whereas stepwise regression and best subset remain attractive options in settings with large samples and small feature spaces, as they require no hyperparameter tuning. Figure 1. Predictive performance measured by normalized RMSE across simulation settings. Bars represent methods, with low-dimensional settings defined as p < 100 and low-SNR settings as SNR < 0.5. Dots indicate model size (right y-axis). 2.4.2 Parameter Estimates Parameter estimates were evaluated using relative test error, which quantifies how closely the estimated regression coefficients align with the true coefficients. A score of 1 indicates perfect recovery of the true coefficients, and lower scores reflect smaller estimation error. Figure 2 summarizes the results, with the dotted black curve representing the null model. Boruta is excluded from this metric because tree-based methods do not yield explicit coefficient estimates. The relative test error patterns are primarily determined by the sample size-to-predictor ratio. When p = 15 (Figure 2a), methods cluster into three groups. The first group is LASSO (min), which consistently achieves the lowest relative test error across conditions. The second group includes stepwise regression, best subset, and SCAD. These methods perform well when correlation is low and SNR is high, but their performance worsens as predictor correlation increases. The third group is LASSO (1se), Ridge, and Elastic Net. These methods perform reasonably when SNR is low but deteriorate sharply as SNR increases, making them the least accurate in this setting. As sample size increases, the performance of the first and second groups converges, while the third group remains suboptimal. The difference between LASSO (min) and LASSO (1se) highlights the critical role of tuning regularization parameters. Although the curves appear well separated in Figure 2a, the actual differences are modest, on the order of 0.05 on the y-axis. When p = 50 (Figure 2b), stepwise regression and best subset are the worst-performing methods at small sample sizes such as N = 100. With larger sample sizes, their accuracy improves and approaches that of LASSO (min). At the same time, Elastic Net, Ridge, and LASSO (1se) appear relatively worse as N grows. The overall pattern resembles the case with p = 15, except that stepwise regression and best subset never fully reach the accuracy of LASSO (min). The case with p = 30 exhibits a pattern similar to that observed for p = 50 and is therefore not shown. When p = 1000 (Figures 2c-d), Ridge regression dominates the plot with very high relative test errors, compressing the lines of other methods. After removing Ridge from the display, it becomes clear that stepwise regression and best subset perform poorly in high dimensions, although their accuracy improves toward that of the regularized methods as SNR and sample size increase. Overall relative test error decreases with larger sample sizes, showing that sufficient data can offset the challenges of high dimensionality. Across these scenarios, LASSO (min) and SCAD provide the most accurate coefficient estimates. SCAD performs best when correlation is low, whereas LASSO (min) excels when correlation is high. The case with p = 100 exhibits a similar pattern to p = 1000. In summary, model accuracy is primarily determined by the sample size to predictor ratio. High SNR improves the performance of stepwise regression and best subset but has limited benefit for regularized methods. Predictor correlation destabilizes stepwise regression and best subset while leaving regularization-based approaches relatively unaffected. LASSO (min) consistently produces the most accurate coefficient estimates, although in settings with large sample size relative to the number of predictors, stepwise regression and best subset achieve similar accuracy without requiring hyperparameter tuning. Figure 6. Heatmap of feature selection across methods. Each row represents a feature and each column a method. Methods are ordered by hierarchical clustering. Figure 6 presents the feature selection heatmap, where each row represents a feature and each column a method, with methods ordered by hierarchical clustering. Four feature clusters, outlined by black boxes, highlight groups of methods with similar selection patterns. LASSO (min) and SCAD produced the largest models, reflecting weaker regularization, whereas LASSO (1se) and Elastic Net selected more similar and parsimonious sets of features under stronger regularization. Stepwise regression and best subset identified largely overlapping features, with stepwise regression selecting more. Boruta, which performed best in simulation studies, displayed the most distinct selection pattern. To explore these differences in detail, Table 3 contrasts features selected by stepwise regression (Buse et al., 2025) versus Boruta. Stepwise regression uniquely identified prior use of SGLT2 inhibitors, analgesics, tirzepatide, maximum dose of GLP-1 receptor agonists, and BMI, features primarily related to hyperglycemia. In contrast, Boruta uniquely selected prior use of diuretics and beta-blockers, features more reflective of hypertension. Given the strong correlation between hyperglycemiaand hypertension-related variables, these differences suggest that distinct modeling approaches may emphasize different aspects of patient clinical profiles. Stepwise regression provides a more straightforward interpretation, as hyperglycemia is directly linked to poorly controlled diabetes, but Boruta raises the possibility that hypertension-related factors may also play an important role in this population. Table 3. Features Selected by Stepwise Regression Compared with Boruta Analysis Feature selected by Stepwise regression (Buse et al., 2025) Feature selected by Boruta Times selected in 1000 replicates Ethnicity = Hispanic/Latino 1000 Number of Antihypertensive Classes 999 Age 807 Fibrates 846 Diuretics 988 Region 909 Beta Blocking Agents 595 SGLT2 inhibitor 471 Analgesics 203 Tirzepatide 171 BMI >= 30 34 Took Maximum Dose of Any GLP-1 3 To better visualize frequently selected features (defined as appearing in at least 50% of replicates) across methods, an UpSet plot was generated (Figure 7). The LASSO (min) model produced the largest set, more than twice as large as the other methods, with 5 unique features and 5 features shared only with SCAD. As noted earlier, four features were consistently identified by all methods. Best subset yielded the most parsimonious model, selecting only these four features, while LASSO (1se) added one additional predictor (SGLT2 inhibitor). Boruta, stepwise regression, and Elastic Net produced models of comparable size, largely selecting subsets of LASSO (min) features. Overall, there was substantial overlap in selected features, though the relative importance ranking varied across methods. Figure 7. UpSet plot of features selected ≥50% of the time across different methods. In summary, the real-data analyses reinforce our simulation findings: while multiple methods achieve strong predictive performance, they differ in the balance between parsimony and stability. Stepwise regression provides interpretable, compact models without the need for hyperparameter tuning and proves to be a practical and effective option in this setting. Meanwhile, Boruta identifies complementary signals that may capture alternative clinical pathways, offering additional insights beyond those of traditional linear approaches. 4. Discussion In this study, we systematically evaluated the performance of stepwise regression, best subset, Boruta, and four regularization-based methods, LASSO, Ridge, SCAD and Elastic Net, across simulated scenarios and real-world clinical data from the CATALYST study. Our aim was to assess predictive performance, feature selection stability, and interpretability under varying sample sizes, feature dimensionality, sparsity patterns, and collinearity levels. From a practical standpoint, our results suggest tailoring method choice to data characteristics may be warranted. When the number of predictors is small and collinearity is low, stepwise and best subset methods offer good predictive performance and variable selection accuracy, with the added advantage of interpretability and no need for hyperparameter tuning. In medium-to-large feature spaces or when collinearity is high, regularized methods, particularly Elastic Net and LASSO (1se), provide more stable and parsimonious models. Boruta is especially useful in correlated settings, where its ability to exploit redundancy among predictors yields strong feature selection performance. By contrast, LASSO (min) and SCAD tend to over-select features, sacrificing interpretability and selection accuracy despite competitive raw predictive performance. Application to the CATALYST study confirmed these patterns. LASSO (min) and Ridge achieved the highest predictive accuracy, while stepwise regression ranked second but yielded a much more parsimonious model. LASSO and stepwise regression largely converged on a core set of predictors for hypercortisolism, whereas Boruta identified additional features that were missed by the linear models. These findings have important implications. First, they underscore the fact that there is no “one size fits all” model. Classic methods such as stepwise regression demonstrate good predictive power when the feature number is small and feature interaction is low. Regularized methods are more stable under high feature dimensions, but hyperparameter tuning is critical, and small differences in these parameters may lead to very different conclusions and models. Second, our results highlight the value of combining complementary approaches. Methods such as Boruta may reveal alternative or correlated predictors that linear models under-emphasize, offering a broader perspective on potential clinical pathways. There are several limitations to note. Our simulation settings, while extensive, cannot capture all complexities of real-world datasets. For instance, nonlinear relationships are common in clinical settings, yet our simulation framework was restricted to linear structures. Furthermore, our evaluation emphasized linear modeling strategies, while more flexible approaches such as gradient boosting or neural networks may deliver additional predictive gains, albeit at the expense of interpretability. In conclusion, this study highlights the practical trade-offs among variable selection methods in clinical research. These insights can help guide the choice of analytic strategy in future discovery and risk prediction efforts, balancing predictive accuracy, model stability, and interpretability in complex clinical populations. Acknowledgements The authors thank Daniel Einhorn for his thoughtful advice and contributions, which have strengthened this work. The authors are also grateful to Nina Pashova for her support in conducting the analyses. In addition, the authors acknowledge the Corcept biometrics team for their insightful comments and discussions, which improved the clarity and rigor of this work. References Buse, J. B., et al. (2025). Prevalence of hypercortisolism in difficult-to-control type 2 diabetes. Diabetes Care, dc242841. Calcagno, V., & de Mazancourt, C. (2010). glmulti: An R package for easy automated model selection with (generalized) linear models. Journal of Statistical Software, 34(12), 1–29. Derksen, S., & Keselman, H. J. (1992). Backward, forward and stepwise automated subset selection algorithms: Frequency of obtaining authentic and noise variables. British Journal of Mathematical and Statistical Psychology, 45(2), 265–282. Fan, J., & Li, R. (2001). Variable selection via nonconcave penalized likelihood and its oracle properties. Journal of the American Statistical Association, 96(456), 1348–1360. Furnival, G. M., & Wilson, R. W. (1974). Regression by leaps and bounds. Technometrics, 16(4), 499–511. Harrell, F. E. (2015). Regression modeling strategies: With applications to linear models, logistic and ordinal regression, and survival analysis (2nd ed.). Springer. Hastie, T., Tibshirani, R., & Tibshirani, R. J. (2017). Extended comparisons of best subset selection, forward stepwise selection, and the lasso. arXiv preprint arXiv:1707.08692. Hoerl, A. E., & Kennard, R. W. (1970). Ridge regression: Biased estimation for nonorthogonal problems. Technometrics, 12(1), 55–67. Kursa, M. B., & Rudnicki, W. R. (2010). Feature selection with the Boruta package. Journal of Statistical Software, 36(11), 1–13. Matthews, B. W. (1975). Comparison of the predicted and observed secondary structure of T4 phage lysozyme. Biochimica et Biophysica Acta, 405(2), 442–451. Miller, A. (2002). Subset selection in regression (2nd ed.). Chapman and Hall/CRC. Tibshirani, R. (1996). Regression shrinkage and selection via the Lasso. Journal of the Royal Statistical Society: Series B (Methodological), 58(1), 267–288. van Rijsbergen, C. J. (1979). Information Retrieval (2nd ed.). Zou, H., & Hastie, T. (2005). Regularization and variable selection via the Elastic Net. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 67(2), 301–320.