Comparison of Lasso and Stepwise Regression in Psychological Data

Stepwise Regression Using the F-Test

Traditionally, variable selection in stepwise regression is carried out based on null hypothesis testing, such as the F-test. In forward selection (i.e., selection starts with an empty model or just the intercept), for example, the p-value for each candidate predictor at each step t is calculated based on the partial F statistic as follows:

Each step selects the predictor with the smallest p-value. This iterative algorithm terminates when no remaining predictors meet the inclusion threshold (i.e., a prespecified $\alpha$ value). The magnitude of this $\alpha$ value determines the selection liberty: a higher $\alpha$ allows for the inclusion of more predictors but increases the false positive rate. Backward elimination does the opposite of forward selection. Starting with a full model that includes all variables, each step prunes a variable until no remaining predictor has a p-value larger than the prespecified inclusion threshold. Both-directional stepwise regression combines forward selection and backward elimination by either adding a variable or pruning a variable at each step; a variable that enters the model at a prior step can also be removed later. Regardless of the direction of variable selection, these methods aim to arrive at a final model that best fits the given data.

Stepwise Regression With Information-Based Fit Indices

Information-based fit indices penalize the number of predictors included in the stepwise regression to maximize model fit with the fewest predictors¹. The two most common indices, namely, AIC and BIC, differ in the degree of penalty they place on the number of predictors. AIC has a fixed penalty of 2, whereas the penalty level of BIC, ln(N), increases with the sample size. As ln(N) is typically larger than 2 in psychological research, BIC generally favors more parsimonious models than AIC. Specifically, at each step t, the difference in AIC and BIC is calculated as follows:

where $K$ is the number of predictors to be estimated in the model. Hence, $K_t-K_{t-1}$ is 1 for forward selection and -1 for backward elimination.

Lasso Regression

The Lasso regression aims to reduce overfitting with regularization, a general approach to reduce prediction error in unseen samples by finding a balance between bias and variance. In the context of regression, bias refers to the distance of the estimated coefficients from the true coefficients, and variance refers to the uncertainty of the coefficient estimates due to sampling variability. The two elements often compete with each other, a phenomenon known as the bias-variance tradeoff (Helwig, 2017; Liu & Rhemtulla, 2022; McNeish, 2015). Out-of-sample prediction errors are affected by both bias and variance. Hence, if the analytical goal is to optimize predictive performance, arriving at a model with minimal bias may not be ideal; instead, it may be helpful to enlarge bias for smaller variance. In regression, this is done by shrinking the coefficient estimates towards zero. Popular regularization methods include Ridge regression, Lasso regression, and Elastic-net methods. Among them, Lasso regression is arguably the most frequently used variable selection method (Tibshirani, 1996; Zou et al., 2007). In addition to minimizing the RSS, Lasso adds a penalty on the sum of the absolute values of the beta coefficients:

This constraint attenuates the coefficients of all predictors, with the degree of attenuation determined by the tuning parameter $\lambda$. As $\lambda$ increases, the coefficients of predictors with small effects on the outcome variable will be quickly shrunk to zero, whereas the coefficients of predictors with larger effects will remain non-zero.

As different $\lambda$ selects different sets of variables, an important task in Lasso regression is to select the optimal value of $\lambda$ to reach a desired level of model sparsity. A too-small $\lambda$ would include almost all variables. A too-large $\lambda$ can result in an under-fitted model with too few predictors. The information-based fit indices AIC or BIC (Zou et al., 2007) and the k-fold cross-validation method (Golub et al., 1979) can be used to determine an appropriate level of $\lambda$. Among these methods, k-fold cross-validation is used more frequently because information-based fit indices require the calculation of degrees of freedom, which is challenging (McNeish, 2015; Tibshirani & Taylor, 2012). K-fold cross-validation splits the dataset into k samples. Each sample takes turns validating the model trained by the rest of the k-1 samples. For each iteration of training and validation, a series of the $\lambda$ values are used to generate and compare the cross-validation mean-squared error (MSE)². The $\lambda$ value that minimizes the cross-validation MSE (hereinafter denoted as $\lambda$_min) is recommended for variable selection. Alternatively, others recommend using the largest $\lambda$ value that produces a cross-validation MSE within one standard error of the minimum cross-validation MSE (hereinafter denoted as $\lambda$_1se) for better model parsimony (Friedman et al., 2010; McNeish, 2015).

Lasso regression has been recently introduced to psychologists (Helwig, 2017; Jacobucci et al., 2019; McNeish, 2015), and has been gaining traction in empirical analyses to enhance the prediction of psychological behaviors (Haynos et al., 2021; Smith et al., 2020). Yet, most methodological studies assessing its performance focused on high-dimensional data (i.e., the number of predictors is far larger than the sample size; Bickel et al., 2009; Fan & Lv, 2010; Sirimongkolkasem & Drikvandi, 2019). In contrast, psychological datasets often contain more observations than variables. Although a few recent studies have found that Lasso is more accurate in predicting the hold-out sample than stepwise regression when applied to low-dimensional datasets (Ahrens et al., 2020; Greenwood et al., 2020; Xu et al., 2012), others have shown that this improvement in prediction is very small and limited to a few conditions (Hastie et al., 2020; Pavlou et al., 2016; Wester et al., 2022). The differences in data characteristics are the primary reason for these mixed findings. Among these studies, two investigated categorical response variables from a single empirical dataset (Greenwood et al., 2020; Xu et al., 2012), whereas others included continuous response variables in simulated datasets (Ahrens et al., 2020; Hastie et al., 2020; Wester et al., 2022). Moreover, their simulated data did not well represent characteristics of data commonly seen in psychological research. For example, Ahrens et al. (2020) and Hastie et al. (2020) examined datasets with 100 candidate predictors, which is rarely the case in psychological research. In Wester et al. (2022), which focuses on the selection of interactions to model treatment effect heterogeneity, relevant and irrelevant predictors were simulated to be independent. However, this assumption is unrealistic, given that correlated covariates are ubiquitous in psychological studies. Therefore, it is still unclear how this novel Lasso regression method compares to the more popular stepwise regression method when applied to psychological data. Given the increasing popularity of the Lasso method in psychology, more research is needed to better understand its performance in typical psychological data, especially its selection accuracy, which Lasso regression was not developed to optimize but is highly desired in psychological research to draw inferential conclusions.

To address this issue, we conducted a Monte Carlo simulation study to compare Lasso and stepwise regression across a representative range of data conditions that are typical in psychological studies. We organize the remainder of this article as follows. First, we describe our simulation study comparing Lasso and stepwise regression in out-of-sample prediction, within-sample selection accuracy, and model sparsity. Next, we demonstrate their differences in an empirical study that aims to identify risk factors of adolescent externalizing problems. Finally, we conclude with cautious notes and recommendations about the application of these methods in practice.

Simulation Design

Based on our literature review of recent psychological studies that involve variable selection methods (see Supplementary Material A in Zhou et al., 2024), we simulated data using a four-way factorial design. Specifically, we manipulated: 1) sample size (N), which can be 100, 200, 400, or 800; 2) candidate predictor size (p), which can be 5, 15, 25, 35, or 80³; 3) signal-to-noise ratio (SNR), defined as $\frac{Var\left(\widehat{y}\right)}{Var\left(y-\widehat{y}\right)}$ or $\frac{R^2}{1-R^2}$, can be 0.2, 0.5, 0.8, or 2; and 4) level of sparsity ($s$), defined as the percentage of non-informative predictors, which can be 20%, 40%, or 80%. In total, the factorial design yielded 240 (4×5×4×3) conditions, with 100 replications of each simulation condition. All data were simulated based on the multiple linear regression model without the intercept:$y=X\beta +ϵ$, where $X$is an $N\times p$matrix of candidate predictor values drawn from $N_p\left(0,\Sigma \right).$ $\Sigma$ is a $p\times p$ correlation matrix with the elements on the diagonal fixed to 1, and$$the other entries randomly and independently drawn from $Beta\left(3.5,3.5\right)$ bounded between -1 and 1. We used this distribution so that most correlations among the candidate predictors ranged between -0.5 and 0.5, but some candidate predictors were allowed to be highly correlated with each other. We ensured that the predictor correlation matrix was positive semi-definite by replacing a non-positive definite correlation matrix with the nearest positive-definite correlation matrix (Higham, 2002). In each replication, a unique correlation matrix $\Sigma$ was used to generate the candidate predictor values. $\beta ={\left({\beta }_1,\dots ,{\beta }_p\right)}^T$ contains regression coefficients with set $J_1=\left\{j,{\beta }_j=1\right\}$ representing the informative predictors and set $J_2=\left\{j, {\beta }_j=0 \right\}$ the non-informative predictors. $y={\left(y_1,\dots ,y_n\right)}^T$ is the outcome variable from $N_n\left(X\beta ,{\sigma }^2{I}_{\text{n}}\right)$ where ${\sigma }^2=\frac{Var\left(X\beta \right)}{SNR}$, and $I_{\text{n}}$ is the $n\times n$ identity matrix. Within each replication, we simulated a training set and a test set using the same parameter values. The test set was used to calculate the out-of-sample predictive accuracy, as described below.

Simulation Analysis

The stepAIC function in the MASS package was used (Venables & Ripley, 2002) to perform forward selection, backward elimination, and both-direction stepwise regression with three different selection criteria including the traditional F-statistic method with α = 0.15 (i.e., the threshold criteria compared to the p-value of each predictor) as recommended in the literature (Derksen & Keselman, 1992; Flack & Chang, 1987), and the information-based fit indices AIC and BIC. We then used the cv.glmnet function in the glmnet package (Friedman et al., 2010) to conduct Lasso regression. Five-fold cross-validation was used to optimize the tuning parameter $\lambda$. We performed Lasso with both $\lambda$_min and $\lambda$_1se (respectively denoted as Lasso.min and Lasso.1se below). Because Lasso is sensitive to the scale of the coefficient values during variable selection, all variables were standardized prior to the analysis (Hastie et al., 2009, Chapter 3)⁴.

Comparison Criteria

Analysis

We conducted a factorial repeated-measures analysis of variance (RM-ANOVA; Myers, 1979) to investigate the main effects of predictor size, sample size, SNR, sparsity, selection method, and their interaction effects on each comparison criterion. We evaluated those effect sizes through the generalized eta squared (${\eta }_G^2$) (Olejnik & Algina, 2003). Because there were many terms in the ANOVA models, we only reported results with at least a small effect size (i.e., ${\eta }_G^2\ge 0.01$; Cohen, 1988), and we focused on interaction effects that involve selection methods. All steps of simulation and analyses were conducted in R Version 4.0.0 (R Core Team, 2020), and the simulation code is in the Supplementary Material B of Zhou et al. (2024).

Preliminary Analyses

All three backward elimination methods (i.e., backward elimination with F-statistic, AIC, and BIC) produced far more complex models with significantly poorer predictive performance and selection accuracy. For example, the MSE of all backward elimination methods with small sample sizes (N < 200) were disproportionately high with values larger than 2000; the MCC of all backward elimination methods were almost 0 when predictor size was larger than 15; increase in sample size did not improve sensitivity as other methods did⁵. We also found that the forward selection methods performed very similarly to the both-directional stepwise regression methods. Hence, we only report the results of both-directional stepwise regression methods below⁶. The descriptive statistics for all methods are included in Supplementary Material C of Zhou et al. (2024), and Table 5 in Supplementary Material C of Zhou et al. (2024) details the ANOVA results for the both-directional stepwise regression methods and the Lasso methods.

Within-Sample Selection Accuracy

Sensitivity

Figure 1 depicts the levels of sensitivity across conditions and methods. Only sample size of 100 is shown here for demonstration because, 1) it is the closest to our empirical data sample size, and 2) sample size has the least influence on all measures without imposing any interaction effects with the analytical method.

The results of different sample sizes are included in the Supplementary Material C of Zhou et al. (2024) for interested readers. There were large main effects of predictor size (${\eta }^2$ = 0.65), sparsity (${\eta }^2$ = 0.4), selection method (${\eta }^2$ = 0.26), SNR (${\eta }^2$ = 0.19), and sample size (${\eta }^2$ = 0.15). There were also small interaction effects between SNR and selection method (${\eta }^2$ = 0.03), between predictor size and selection method (${\eta }^2$ = 0.03), and a three-way interaction effect between predictor size, sparsity, and selection method (${\eta }^2$ = 0.02). Specifically, sensitivity increased as sparsity, SNR, and sample size increased (see Table 2.2 in the Supplementary Material C of Zhou et al., 2024) but decreased as predictor size increased. The Lasso methods, especially Lasso.min had the greatest sensitivity across all conditions. Lasso.1se had higher sensitivity than stepwise regression unless SNR was as small as 0.2. This means that the Lasso methods were generally better at identifying the informative predictors than the stepwise regression methods unless the data were very noisy. Furthermore, the difference in sensitivity between Lasso and stepwise regression was minimal when only a few candidate predictors were considered ($p$ = 5), and only a few were informative (s = 0.8). This difference was larger when predictor size and the number of informative predictors increased.

Click to enlarge

Figure 1

Average Sensitivity Across Methods and Across Conditions When Sample Size Was 100

Specificity

Similar to sensitivity, all of the five factors of interest: selection method (${\eta }^2$ = 0.17), predictor size (${\eta }^2$ = 0.07), sparsity (${\eta }^2$ = 0.06), SNR (${\eta }^2$ = 0.04), and sample size (${\eta }^2$ = 0.02), had at least small effects on specificity. The interaction between predictor size and selection method yielded a small effect (${\eta }^2$ = 0.02). Figure 2 depicts the levels of specificity across conditions and methods.

Larger predictor sizes and sparsity levels, lower SNR, and smaller sample sizes (see Table 2.3 in the Supplementary Material C of Zhou et al., 2024) were associated with a higher level of specificity. Lasso.min had the lowest specificity across all conditions, and Lasso.1se outperformed stepwise regression methods only when SNR was lower than 0.5. This indicates that Lasso.min was more likely to classify non-informative predictors as significant and thus produced models with larger false positive or Type 1 error rate than all others. Lasso.1se was better but did not outmatch stepwise regression methods unless in conditions wherein only a small proportion of the outcome variance was explainable by the linear regression. Moreover, the difference in specificity across methods was magnified in conditions with only five predictors. For example, Lasso.min’s degree of specificity was almost twice as small as most other methods, and the specificity of Lasso.1se was the highest.

Click to enlarge

Figure 2

Average Specificity Across Methods and Across Conditions When Sample Size Was 100

MCC

Sparsity level (${\eta }^2$ = 0.3) and predictor size (${\eta }^2$ = 0.22) had a large effect size when predicting MCC. Selection method (${\eta }^2$ = 0.02) and the three-way interaction effect between predictor size, sparsity, and selection method (${\eta }^2$ = 0.02) exhibited small effects on MCC. Figure 3 depicts the levels of MCC across conditions and methods.

In general, MCC increased with higher sparsity levels and smaller candidate predictor sizes. This implies that selection is generally more accurate when dealing with a small pool of candidate variables of which only a few are truly informative. Across methods, although Lasso.min had the lowest MCC when predictor size was small ($p$ < 25) and only a few were informative (s = 0.8), Lasso regression methods had higher MCC across most other conditions. In particular, Lasso.1se had the highest MCC across all conditions. This indicates that Lasso regression methods, especially Lasso with a harsher shrinkage tuning parameter, were more likely to produce models with better selection accuracy than stepwise regression methods.

Click to enlarge

Figure 3

Average MCC Across Methods and Across Conditions When Sample Size Was 100

Model Sparsity

All the five main effects yielded large effect size: predictor size (${\eta }^2$ = 0.63), selection method (${\eta }^2$ = 0.47), sparsity (${\eta }^2$ = 0.26), SNR (${\eta }^2$ = 0.24), and sample size (${\eta }^2$ = 0.17). There were also two small-size interaction effects with the selection method: SNR and selection method (${\eta }^2$ = 0.04), and predictor size and selection method (${\eta }^2$ = 0.03). Figure 4 depicts differences in model sparsity across conditions and methods.

Model sparsity generally increased when predictor size and sparsity level increased; it decreased when sample size (see Table 2.5 in Supplementary Material C of Zhou et al., 2024) and SNR increased. This means that for all methods, a more parsimonious model was more likely to be obtained if the candidate predictor pool was large, a few predictors were truly relevant, the sample size was small, and the data were noisy. Between methods, Lasso.min produced the least parsimonious models across all conditions. On the other hand, stepwise regression with BIC produced more parsimonious models than all other methods in most conditions except when the SNR was 0.2 or when the predictor size was 5. In such cases, Lasso.1se produced models similar in size or more parsimonious than stepwise regression with BIC.

Click to enlarge

Figure 4

Average Model Sparsity Across Methods and Across Conditions When Sample Size Was 100

The Empirical Data

We demonstrate the variable selection differences between Lasso and stepwise regression in an empirical study aimed at predicting externalizing behavior (i.e., disruptive, aggressive, or delinquent actions that are directed outwardly) during adolescence. The data were drawn from the second timepoint of the Stanford Early Life Stress Study (Chahal et al., 2022; Gotlib et al., 2021), where externalizing behavior is assessed by the aggressive behavior and rule-breaking behavior scales of the Youth Self-Report (YSR; Achenbach, 2001), a widely used measure of behavioral problems in adolescents. Thirty-three externalizing-behavioral–relevant variables from the domains of pubertal development, sensitivity to stress and reward, emotional and behavioral problems, social support, physical and emotional neglect, emotional and affective regulation, early life stress severity, and demographic information of the child and parents, were pre-selected as predictors of externalizing behavior. Complete case analysis was conducted with a sample size of 141. A detailed description of the measures of the predictors and the sample is included in the Supplementary Material D of Zhou et al. (2024).

Analysis

To examine the predictive performance of the models, we split the full dataset into a training set with 80% of the sample and a test set with the remaining 20%. This random splitting procedure was repeated 1000 times to investigate selection variability due to sample variation. Methods were compared in terms of their average out-of-sample predictive accuracy (i.e., MSE), average model sparsity (i.e., the percentage of variables not selected across the 1000 iterations), and selection rate of each variable.

Empirical Study Results

Summary of Selection

Figure 5 shows the selection rate of each candidate predictor. Internalizing problems, impulsivity, affective reactivity, sex (being male), and sensitivity to punishment were selected most often by all methods. Among these methods, stepwise regression with BIC, stepwise regression with F statistic, and Lasso.1se tended to produce more parsimonious models than stepwise regression with AIC and Lasso.min. Moreover, they were more consistent in their selection, an indicator that these three methods were more robust to sample variability. It should be noted that, however, this discrepancy in selection consistency is confounded by model sparsity.

Click to enlarge

Figure 5

Selection Rate of Each Candidate Variable Across Methods

Cautious Notes and Recommendations in Practice

We recommend researchers to carefully consider their analytical goals and data characteristics before choosing a variable selection method. Lasso regression is generally more consistent in out-of-sample prediction, although its improvement in magnitude in typical psychological data is not as pronounced as it is in high-dimensional datasets. Lasso regression is also a competitive alternative to stepwise regression if the aim is to make accurate explanation for its better within-sample selection accuracy. However, neither of the two Lasso methods investigated here exhibited sufficiently high selection accuracy across all conditions. As shown in Figure 3, even for Lasso.1se, the one with the highest MCC, its MCC hardly exceeded 0.4 unless the predictor size was very small (p = 5) or the sparsity level was very high (s = 0.8). Because the sparsity level of the true model is unknown in empirical data, researchers may consider prescreening predictors based on theoretical knowledge or prior empirical evidence to narrow down the candidate pool to a size smaller than 15, if possible.

Given the tradeoff between power and sensitivity in variable selection, the ideal method should consider the analytical priority, such as to capture true effects or to avoid spurious findings. If statistical power is more of a concern, we recommend Lasso regression with the harsher penalty term (${\lambda }_{1se}$) for lower false positive rate, better overall selection accuracy, lower model complexity, and more consistent selection. On the other hand, when the priority is to reach lower false positive rate and to improve model interpretability, Lasso does not yield better results than stepwise regression with information-based fit indices, especially the BIC, particularly when the predictor size is small and only a small fraction of the predictors is truly relevant. To improve the false positive rate of the Lasso regression, we encourage researchers to consider more advanced derivatives such as the adaptive Lasso (Fan & Li, 2001; Zou, 2006) and the relaxed Lasso (Meinshausen, 2007). These methods further exclude less influential variables through their weighted degree of shrinkage that imposes a magnified penalty on non-informative predictors.

Following the machine learning literature, we recommend researchers to check the IRC assumption to evaluate the selection consistency when using Lasso (Meinshausen & Bühlmann, 2006; Zhao & Yu, 2006). It is important to note that the presence of multicollinearity—a common issue in psychological data where predictors are highly linearly dependent—can worsen the selection consistency of Lasso regression. This is because Lasso randomly selects one variable out of a group of highly correlated ones to reach a sparser model (Zou & Hastie, 2005). If researchers are interested in delivering a consistent model that contains all highly correlated predictors to include both interaction and main effects, derivatives of Lasso such as the elastic-net (Zou & Hastie, 2005) and the grouped Lasso (Yuan & Lin, 2006) are generally recommended. Another remedy for the unstable selection is to integrate Bayesian selection into variable selection. For example, the stochastic search variable selection (SSVS) assigns higher prior probability to more promising predictors (Bainter et al., 2020).