A comparative study of different variable selection methods based on numerical simulation and empirical analysis

Lasso method and basic characteristics

Tibshirani (1996) introduced a technique for selecting variables. The fundamental concept underlying Lasso involves the imposition of penalty factors to restrict the L₁ norm of the estimator β, building upon the framework of Ordinary Least Square (OLS) estimation (Tibshirani, 1996).

As for the data set with p predictive variables and n predictive variable-response variable pairs ${\left\{\left(x_i,y_i\right)\right\}}_{i=1}^n$, we could use Lasso to find an estimate of $\hat{\beta }$ that can better fit the data through minimization of $RSS\left(\hat{\beta }\right)$. (4)${\displaystyle RSS\left(\hat{\beta }\right)=\underset{\hat{\beta }}{\underbrace{argmin}}\sum _{i=1}^n{\left(y_i-\sum _{j=1}^p{\beta }_j{x}_{ij}\right)}^2+\lambda \sum _{j=1}^p\left|{\beta }_j\right|}$ λ ≥ 0 is the Penalized Parameter, $\lambda {\sum }_{j=1}^p\left|{\beta }_j\right|$ is compression constraints. (5)${\displaystyle \hat{\beta }=argmin\left\{{∥y-x\beta ∥}_2^2\right\}s.t∥\beta ∥\le t.}$

t ≥ 0 is a tuning Parameter, which controls the intensity of compression. If the Parameter obtained by the least squares is estimated as β⁰, Lasso can be compression as long as $t<{\sum }_{j=1}^p\left|\widehat{{\beta }_j^0}\right|$. In addition, for some models with small absolute values, the coefficient is compressed to zero. Therefore, the inequality $∥\beta ∥\le t$ effectively limits the parameter space and makes the final model explicable.

It can be proved that the relation between t and λ (Wu, 2015) is (6)${\displaystyle t=\sum _{j=1}^{p}sign\left(\widehat{{\beta }_j^0}\right)\cdot \widehat{{\beta }_j^0}-p\frac{\lambda }{2}.}$ The estimator $\hat{\beta }$ can be computed through numerical approximation techniques, such as coordinate descent and minimum angle regression.

The fundamental attributes of Lasso entail pursuing the sparse representation of the model, essentially. The aforementioned procedure is executed through the optimization of a function problem that comprises a combination of a loss and a penalty term. In cases where a set of predictive variables exhibit high correlation, Lasso regression may select one variable while disregarding the others, leading to unstable outcomes. This indicates that collinearity cannot be effectively and accurately addressed.

Adaptive Lasso method and basic characteristics

It can be seen from Lasso analysis that Lasso penalty function is the same for all the estimators in β = (β₁, β₂, …⋯, β_p). However, if the penalty function can be made to give a smaller penalty to the larger estimators, thus effectively reducing the model’s bias while ensuring the sparsity of the model. Therefore, Zou (2006) proposed Adaptive Lasso, in which the punishment parameter can be adjusted according to the size of the estimator, instead of λ being a fixed punishment in Lasso (Liang & Feng, 2016). Before the Adaptive Lasso in $∥{\beta }_j∥$ increase weight coefficient ${\hat{w}}_j=\frac{1}{{\left|{\hat{\beta }}_j\right|}^{\gamma }},\gamma$ constant is greater than zero, Zou advised to choose a $\sqrt{n}$ consistent estimator satisfy, which uses least squares estimator $\hat{\beta }$ as the initial estimator, ${\hat{\beta }}_{\mathrm{j}}={\hat{\beta }}^{\mathrm{OLS}}$, Adaptive Lasso is in the form of (7)${\displaystyle \hat{\beta }=\arg \underset{\beta }{min}{∥y-X\beta ∥}_2^2+\lambda \sum _{j=1}^p{\hat{w}}_j{∥{\beta }_j∥}_1.}$

By formula Eq. (7), when γ ≥ 1, if some components of the initial estimate of $\hat{\beta }$ are large (such as ${\hat{\beta }}_j$), the weight of the corresponding punishment item weight coefficient ${\hat{w}}_j=\frac{1}{{\left|{\hat{\beta }}_j\right|}^{\gamma }}$ relative to other components, so the Adaptive - Lasso has the following features: (1) The utilization of a penalty function serves to decrease the bias of the Lasso estimator, while simultaneously guaranteeing the accuracy of the estimator.

(2) It also proves that Adaptive Lasso can satisfy the following Oracle properties with an appropriate λ value:

(a) The consistency of variable selection:

$\left\{j,{\hat{\beta }}_j\ne 0\right\}=\left\{j,{\beta }_j\ne 0\right\}\triangleq S_0$. (b) Asymptotic normality:

$\sqrt{n}\left(\hat{\beta }-{\beta }_{S_0}\right)\stackrel{d}{\to }N\left(0,{\sigma }^2\right)$.

* S₀ is the active set of predictive variables

$S_0=\left\{j:{\beta }_j^0\ne 0,j=1,2,\dots \dots ,p\right\}$,

$\hat{S}\left(\lambda \right)=\left\{j:{\hat{\beta }}_j\left(\lambda \right)\ne 0,j=1,2,\dots \dots ,p\right\}$,

$\hat{S}\left(\lambda \right)$ Is a nonzero coefficient subscript set of parameters estimated by Lasso method.

According to the CV method to get ${\hat{\lambda }}_{CV}$, eventually get $\hat{S}\left({\hat{\lambda }}_{CV}\right)$ which has a high probability to include S ₀ and $\left|\hat{S}\left({\hat{\lambda }}_{CV}\right)\le min\left(n,p\right)\right|$.

(3) The Adaptive Lasso method conforms to the three characteristics of the penalty function introduced by Fan and Li, namely impartiality, parsimony, and continuity. These attributes are an advancement over the Lasso method. However, the Adaptive Lasso method still encounters challenges in effectively addressing collinearity.

Elastic Net method and basic features

The Ridge regression technique is distinguished by its ability to evenly allocate weights to associated characteristic variables. Conversely, the Lasso method is not capable of effectively addressing collinearity. In response, it introduced the Elastic Net approach, which combines the strengths of both Ridge and Lasso regression techniques. Specifically, Elastic Net is a convex combination of Ridge and Lasso regression (Rodgers & Nicewander, 1988). (8)${\displaystyle \hat{\beta }=arg\underset{\beta \in R^p}{min}\left\{{∥y-X\beta ∥}_2^2+{\lambda }_1\sum _{j=1}^p\left|{\beta }_j\right|+{\lambda }_2\sum _{j=1}^p{\beta }_j^2\right\}}$

The penalty term is $P_{\lambda }\left(\left|{\beta }_j\right|\right)={\lambda }_1\left|{\beta }_j\right|+{\lambda }_2{\left|{\beta }_j\right|}^2$, λ₁, λ₂ are two nonnegative punishment parameters, order $a=\frac{{\lambda }_2}{{\lambda }_1+{\lambda }_2}$, (9)${\displaystyle \hat{\beta }=arg\underset{\beta \in R^p}{min}\left\{{∥y-X\beta ∥}_2^2+\left(1-a\right)\sum _{j=1}^p\left|{\beta }_j\right|+a\sum _{j=1}^p{\beta }_j^2\right\}}$

Penalty term $P_{\lambda }\left(\left|{\beta }_j\right|\right)=\left(1-a\right){\sum }_{j=1}^p\left|{\beta }_j\right|+a{\sum }_{j=1}^p{\beta }_j^2$ at zero point is the guide, for ∀a > 0 is strictly convex.

To: $X_{\left(n+p\right)\times p}^{\ast }=\frac{1}{\sqrt{1+{\lambda }_2}}\left(\genfrac{}{}{0.0pt}{}{X}{\sqrt{{\lambda }_2}I}\right)$,  $Y_{n+p}^{\ast }=\left(\genfrac{}{}{0.0pt}{}{Y}{0}\right)$,

$\gamma =\frac{{\lambda }_1}{\sqrt{1+{\lambda }_2}}$,   ${\beta }^{\ast }=\sqrt{1+{\lambda }_2}\beta$,

then

${\hat{\beta }}^{\ast }=arg{min}_{{\beta }^{\ast }\in R^p}\left\{{∥y^{\ast }-X^{\ast }{\beta }^{\ast }∥}_2^2+\gamma {\sum }_{j=1}^p\left|{\beta }_j^{\ast }\right|\right\}$.

Therefore, the Elastic Net problem can be converted to the Lasso problem. The optimal Elastic Net solution ${\hat{\beta }}^{\ast }$ can be found on $\hat{\beta }=\frac{1}{\sqrt{1+{\lambda }_2}}{\hat{\beta }}^{\ast }$. Because I_p×p is full rank, considering the composition of $X_{\left(n+p\right)\times p}^{\ast }$ and it has p columns, the Lasso matrix $X_{\left(n+p\right)\times p}^{\ast }$ (matrix rank p), the applied matrix is X_(n+p)×p. Therefore, Elastic Net can select up to p variables which solved Lasso can only select a maximum of n (n < p) variables. According to the findings, the Elastic Net technique is capable of selecting group variables, meaning that it can effectively select multiple highly correlated independent variables. This addresses the issue commonly encountered with Lasso regression methods, which tend to select only one variable and disregard the others. This approach is appropriate for certain application scenarios that require analysis of the correlation between dependent variables and predictive variables exhibiting group characteristics. Examples of such scenarios include variable screening and prediction of gene expression profile data.

SCAD method and basic features

In their study, Fan & Li (2001) introduced a non-convex penalty function with the aim of attaining unbiasedness of the estimator β. This penalty function ensures that the penalty on the coefficient decreases as the coefficient estimator increases, thereby ensuring the approximate unbiasedness of large coefficients. (10)${\displaystyle \hat{\beta }=arg\underset{\beta ϵR^p}{min}\left\{{∥y-X\beta ∥}_2^2+\sum _{j=1}^p{P}_{\lambda }\left(\left|{\beta }_j\right|\right)\right\},}$

$P_{\lambda }\left(\left|{\beta }_j\right|\right)P$ is the penalty term. The proposed method is a continuously differentiable penalty function, called smooth cohesion absolute deviation penalty SCAD (Smoothly Clipped Absolute Deviation Penalty). (11)${\displaystyle P\left({\beta }_j{|}_{\lambda ,a}\right)=\left\{\begin{array}{c}\lambda \left|{\beta }_j\right|\left|{\beta }_j\right|\le \lambda ,\hfill \\ \frac{2a\lambda \left|{\beta }_j\right|-{\left|{\beta }_j\right|}^2-{\lambda }^2}{2\left(a-1\right)}\lambda <\left|{\beta }_j\right|<a\lambda ,\hfill \\ \frac{\left(a+1\right){\lambda }^2}{2}\left|{\beta }_j\right|\ge a\lambda ,\hfill \end{array}\right.}$

Take the derivative of the penalty function SCAD, get (12)${\displaystyle P^{{}^{\prime }}\left({\beta }_j{|}_{\lambda ,a}\right)=\left\{\begin{array}{c}sgn\left({\beta }_j\right)\lambda \left|{\beta }_j\right|\le \lambda \hfill \\ sgn\left({\beta }_j\right)\left(\alpha \lambda -\left|{\beta }_j\right|\right)/\left(\alpha -1\right)\lambda <\left|{\beta }_j\right|<a\lambda \hfill \\ 0\left|{\beta }_j\right|\ge a\lambda \hfill \end{array}\right.}$

Thus, its coefficient estimate speed of punished decreases with the coefficient estimators $\left|{\beta }_j\right|$ increases, namely when the β_j absolute value is larger, SCAD penalty is constant, when the β_j absolute value is small, the smaller the absolute value of the coefficient of the degree of compression than the LASSO method is bigger and therefore better able to induce the sparse structure, t is convenient for screening and obtaining a sparse subset of variables. Therefore, SCAD guarantees the unbiasedness of large coefficients. Not differentiable at the origin, but continuously differentiable at (−∞,0) ∪ (0,+ ∞), a local quadratic approximation is used to obtain a local optimal solution while ensuring sparsity and continuity. Fan, Huang, and Kim respectively proved that N-SCAD (conventional), H-SCAD (high dimensional), and UH-SCAD (super high dimensional) have Oracle properties under certain hypothetical conditions, and have the advantage over Adaptive Lasso in that it does not require prior access to the consistent estimator of $\sqrt{n}$ to ${\hat{\beta }}^0$. It is important to note that the local optimal solution may not always coincide with the global optimal solution. Numerical methods often yield a particular local solution, thereby creating a disparity between theoretical and practical applications. Consequently, research in this area remains a prominent topic of interest.

Numerical simulation and comparative analysis

This study conducts a comparative analysis of the penalty function of Lasso, Adaptive Lasso, Elastic-Net, and SCAD variable selection methods under varying conditions through numerical simulation. The study also examines the applicable application scenarios of these methods.

The model’s performance was assessed using four indicators, namely: (1) the consistency of model variable selection, (2) the mean square error, (3) the efficiency of the model algorithm, and (4) the stability of the model. Regarding the evaluation of the consistency of model variable selection, the false discovery rate and false exclusion rate methods as well as the calculation of the probability of predictive variables entering the model were not utilized in our approach (Jing, 2013; Wood, 2018; Xu, 2022; Min, 2022; Chu, 2018; Hastie, Tibshirani & Wainwright, 2015; Bai & Zhao, 2005). However, we choose variable selection compatibility of the model to be found through calculating the Angle between uniformed estimated coefficient vector $\hat{\beta }$ and the uniformed true coefficient vector β $\text{Angle}=\frac{180}{\pi }arccos\left({\beta }^T\hat{\beta }\right)$ reflection coefficient estimation accuracy, namely consistency ($\left\{j:\widehat{{\beta }_j}\ne 0\right\}=\left\{j:{\beta }_j\ne 0\right\}\triangleq S_0$; S₀ is the active set of predictive variables). The Angle is closer to zero, the higher the degree of compatibility between β and $\hat{\beta }$, the stronger the model interpretation ability (geometric interpretation of Pearson correlation coefficient formula) (Liang, Feng & Song, 2016). $Err=\frac{1}{n}\sum {\left({\hat{y}}_0-y_0\right)}^2$, the mean square error of the predicted variable obtained from the training set which can be carried to the model, the mean square error of the response variable ${\hat{y}}_0$ and the actual response value y₀ in the test set can reflect the prediction accuracy of the model. If the predicted response value is close to the actual response value, ΔErr will be small; the stability of the model was assessed by utilizing the boxplot tool to examine the outliers of Angle and Err across various variable selection methods. This approach provided insight into the robustness of the established model. The evaluation of the model algorithm’s efficiency is based on the duration of time taken by the computer to perform a regression fitting.

Four indicators were used to evaluate the model performance:

Angle:consistency of model variable selection;

—- Evaluate the explanatory power of the model

Err:mean square error;

—- Evaluate the prediction accuracy of the model

Outliers:Evaluate the stability of model

Time:Evaluate the efficiency of model algorithm

Design of simulated data set

Our data are generated from the following linear random effects model (13)${\displaystyle Y=X\beta +Z\delta +\varepsilon }$ y_n×1 is the response variable, X_n×p is the design matrix of the prediction variable, Z_n×q is the design matrix of the random effect variable, which is constructed in the same way as X, q ≤ p.β_p×1 is a nonrandom parameter vector called deterministic effect, δ_q×1 is a random parameter vector called random effect, ɛ_n×1 is a random noise vector. δ follows the normal distribution with zero mean and G covariance matrix. δ ∼ N(0, G), ɛ follows normal distribution, ɛ ∼ N(0, R), R = σ²I_n , $cov\left(\delta ,\varepsilon \right)=0$. In order to design the simulation data set, the function mvrnorm() in MASS package is firstly used to design the n × p dimensional multivariate normal simulation matrix about the prediction variable X, X_n×p is the prediction variable, let ${\beta }_{p\times 1}=\left(1,1.2,3.2,4.1,-5,\underset{45}{\underbrace{0,\dots 0}}\right)$- is the coefficient vector, and y_n×1 is the response variable. The model discussed is the sparse model, ɛ_n×1 is the random noise vector, which is independently identically distributed and follows the standard normal distribution, ${\varepsilon }_{n\times 1}\sim N\left(0,I_n\right)$, X and ɛ are independent. The numerical simulation experiment was repeated for 300 times, and the data in the X matrix was updated each run. In this study, the mean value of the prediction variables in each column of the X matrix was designed to be a random variable mu ∼N (0, I_n) subject to normal distribution instead of constant, to simulate the internal fluctuation of the prediction variables and the non-homogeneity of the random factors Zδ between individuals. The fitted model has linear random effects. The extent of collinearity exhibited by the predictor variable X is strongly associated with the correlation coefficient that exists between them. In order to assess collinearity, experimental designs are implemented that incorporate correlation coefficients. A correlation coefficient of 1 corresponds to complete collinearity, while a coefficient of 0 indicates the absence of collinearity. As the correlation coefficient approaches 1, the degree of collinearity between variables increases. The dataset was partitioned into two subsets, namely the training set and the test set. Four variable selection techniques, namely Lasso, Elastic-Net, SCAD, and Adaptive Lasso, were employed to perform regression modeling on the dataset. The evaluation model was based on the average value of four indicators obtained from conducting a numerical simulation experiment for a total of 300 times.

The numerical simulation is carried out in the following two scenarios.

(1) n ≫ p, where the number of observed samples n is much larger than the conventional data set with the number of predictive variables P, n:p = 2500:50; The correlation coefficient between the predictive variables was cor = 0.2/0.4/0.6/0.8/0.98.

(2) high dimensional data set which satisfying $\mathit{log}\left(p\right)=n^a\left(0<a<1\right)$

n:p =40:50; Correlation coefficient cor =0.2/0.4/0.6/0.8/0.98

We use R to analyze and complete numerical simulation, mainly using the following software packages that MASS, car, psych, glmnet, ncvreg, msgps, ISLR etc. We select parameters in software analysis:

For SCAD, we used ncvreg(), where the gamma parameter selection is 3.7 (i.e., the punishment function where a is recommended by Fan); alpha = 1 (MCP/SCAD penalty); When n<p, lambda.min is set to 0.05 instead of 0.001;

For Adaptive Lasso, we used msgps(): where gamma parameter selection is 1 (i.e., γ in the penalty function); Initial estimators, as suggested by Zou when n ≫ p, select the minimum square estimator that meets the requirements of $\sqrt{n}$ convergence ${\hat{\beta }}_{\mathrm{j}}={\hat{\beta }}^{\mathrm{OLS}}$ as the initial estimator, corresponding to the parameter Lambda = 0; If n<p, Lambda = 0.001 as the initial estimate of ${\hat{\beta }}_{\mathrm{j}}={\hat{\beta }}^{ridge}$.

For Elastic-net, we used msgps() where alpha =0.5;

For Lasso, we used msgps() where alpha defaults to 1;

For SCAD, we used local quadratic approximation to obtain the locally optimal solution. We can obtain the Lambda of Lasso, Adaptive Lasso and Elastic-net that used Generalized-Cross-Validationg(GCV) to solve the optimal value and obtain the corresponding coefficient value $\hat{\beta }$.

Numerical simulation result analysis

(1) First examine n:p =2500:50 (number of samples on the training set: number of predictors); Number of trial repeats m =300, the correlation coefficient is cor =0.2/0. 4/0.6/0.8/0.98

Numerical simulation obtains Angle, Err boxplot (Fig. 1), the time data (Table 1) and Angle median-correlation coefficient, Err median ∼correlation coefficient relationship curve (Fig. 2).

(2) Comparative analysis. In the context of a large sample size (n ≫p) and without accounting for model stability, this study examines the accuracy of model predictions and the coincidence of model coefficients when the correlation coefficient is less than 0.8. According to the results, SCAD exhibited the highest performance, trailed by Adaptive-Lasso and Lasso. The reason for this is that SCAD is capable of attaining an approximately unbiased estimation of large coefficients through the design of its penalty function. Theoretically, Adaptive-Lasso can improve variable selectivity by weighting corrections to the Lasso penalty terms, However, the results of this numerical simulation experiment show that the two perform basically the same in terms of variable selectivity (For Adaptive Lasso iterative solution of the minimum loss function, Zou (2006) proposed that to obtain satisfactory Oracle properties which required an initial estimator ${\hat{\beta }}^0$ with fully satisfies the $\sqrt{n}$ consistency. It is difficult to achieved whaterver using coordinate descent method or gradient descent method for iterative solution. The models fitted in these four methods all perform well in terms of prediction accuracy (Err ≤1.005), where the irreducible error ɛ∼N(0,1) is an upper bound on the prediction accuracy of y, maximum value is 1 (i.e., ΔErr ≤ 0.005). In cases where a group of predictors exhibit strong correlation (i.e., cor>0.8), the model’s variable selection compatibility Angle experiences a rapid decline. However, the model’s Err remains relatively small and is largely unaffected by the increasing impact of the correlation coefficient. Hence, the model constructed amidst the presence of pronounced collinearity among the predictor variables is apt for forecasting the outcome based on the predictor variables. However, it is not conducive for making inferences regarding the association between each predictor variable and the response variable.

The evaluation of model stability and algorithm efficiency can be performed through the utilization of Angle, Err boxplot, and regression fitting timetable Time. The convex penalty function possesses favourable stability properties that can enhance the performance of models constructed through the utilization of Lasso, Elastic-Net, and Adaptive-Lasso techniques.

The Lasso algorithm is characterized by its simplicity and efficiency, which results in a notably high speed of model fitting. Thus, the Lasso method is highly appropriate for modeling in scenarios where the number of observations is greater than the number of predictor variables, and the degree of collinearity among the predictor variables is minimal. Elastic-Net exhibits superior stability under conditions of high correlation coefficients, rendering it appropriate for application scenarios characterized by predictive variables possessing group attributes. The complexity of SCAD’s penalty functions, owing to their non-convex nature, results in a model form that is more intricate than that of Lasso, Elastic-Net, and Adaptive Lasso. Additionally, the iterative algorithm employed in SCAD runs notably slower, and the instability of the local quadratic approximation algorithm is evidenced by the Angle and Err in the boxplot with the highest number of outliers. Furthermore, models fitted using this method exhibit poor stability. (3) Next, consider the n:p =40:50 (number of samples on the training set: number of predictors) high-dimensional dataset; The number of test replicates m =300 and the correlation coefficient was cor =0.2/0. 4/0.6/0.8/0.98. The study conducted a numerical simulation and generated several graphical and tabular outputs, including a boxplot of Angle and Err (Fig. 3), a table of time data (Table 2), and a curve depicting the relationship between the median of Angle-correlation coefficient and Err-correlation coefficient (Fig. 4). (4) Comparative analysis. Upon analyzing the Angle and Err indicators, it can be concluded that SCAD exhibits superior performance. However, when compared to the state of affairs pertaining to a large sample data, the Angle indicator exhibits a significant increase by several to dozens of times. The interpretability of models constructed by the four classical methods is notably weak when dealing with high-dimensional data. The Lasso, SCAD, and Adaptive-Lasso variable selection methods have produced models with prediction errors of less than 3. These models may be suitable for application scenarios where high prediction accuracy is not a requirement. In cases where the correlation coefficient between variables exceeds 0.8, the respective indicators of the model experience a rapid decline, thereby rendering the model’s guiding significance null and void.

The unsuitability of the Lasso-like classical variable selection method, which is appropriate for conventional data, for modeling high-dimensional data has been demonstrated through numerical simulation. Over the past decade, there has been a surge in research on high-dimensional models in the field of statistics both domestically and internationally. Numerous modeling theories and application practices have surfaced for various forms of high-dimensional data (Niu & Li, 2021). These have proven to be a potent means of managing vast quantities of information during the age of big data.

Comparative experiment and result discussion

The evaluation of the quality of the fitted model for Lasso, Elastic-Net, Adaptive-Lasso, and SCAD variable selection methods, as per literature, involves the computation of the average of correctly selected predictive variables (“cc”), the average of correctly eliminated unrelated predictive variables (“cd”), and the accuracy of coefficient estimation $\hat{\beta }$ (Rodgers & Nicewander, 1988). The evaluation of the efficacy of the fitted model in literature is conducted through a comparison of the test mean square error (PMSE), the count of non-zero coefficients that are erroneously identified as zero (“IC”), and the count of zero coefficients that are correctly identified as zero (“C”) (Yanbiao, 2011). The literature employs the use of SIS as a means of reducing the dimensionality of high-dimensional data with varying correlation numbers. Subsequently, the fitting process is carried out using Lasso, Adaptive-Lasso, Elastic-Net, and SCAD methods. The resulting regression coefficients $\hat{\beta }$ of each model are then compared to assess the efficacy of the fitting model. All of the findings indicate that the most effective approach for SCAD fitting involves the utilization of non-convex penalty function regression. From a new perspective, this study proposes to calculate the $\text{Angle}=\frac{180}{\pi }arccos\left({\beta }^{\mathrm{T}}\hat{\beta }\right)$ between the unitized estimation coefficient vector $\hat{\beta }$ and the unitized true coefficient vector βto evaluate the consistency of variable selection of the fitted model, and analyze with the help of boxplot tool the distribution and stability of Angle. In order to demonstrate the efficacy and practicality of the proposed approach, a comparative study is undertaken between the model evaluation technique presented in this manuscript and conventional analytical methods, under identical conditions.

The linear mixed random effects model is described as Y  = Xβ + Zδ + ɛ. This study discusses the sparse model. The coefficient vector of the real model is set as $\beta =\left(1,1.2,-3.2,4.1,-5,\underset{45}{\underbrace{0,\dots ,0}}\right)$, in which the number of non-zero coefficients is 5, and the number of zero coefficients is 45. The study conducted numerical simulation experiments in the context of large sample and high dimension. The experiments involved varying correlation coefficients of 0.2, 0.6, and 0.98, and were repeated 300 times. Please refer to Tables 3 and 4 for the obtained outcomes.

Cor: correlation coefficient between predictive variables.

Angle(Median, Outliers): median and number of outliers are given to evaluate respectively the explanatory power and stability of the model

Err (Median , Outliers):(Mean square error of the model)

median and number of outliers are given to evaluate respectively the prediction accuracy and stability of the model

CC: the mean of the number of non-zero coefficients selected

CD: the mean of the number of zero coefficients selected

CEMR: the average error of coefficient estimation $\hat{\beta }$

Table 3 presents the four indicators proposed by this study, namely Angle Median, Angle Outliers, Err Median, and Err Outliers, along with the three commonly used indicators in the research field pertaining to correctly selected predictive variables (represented by “CC”), correctly eliminated irrelevant variables (represented by “CD”), and the average error of coefficient estimation $\hat{\beta }$ (represented by “CEMR”), under the condition of a large sample (n ≫ p).The numerical simulation results indicate that: (1) Of the four variable selection methods, namely Lasso, Elastic-Net, Adaptive-Lasso, and SCAD, the latter exhibits superior performance. Specifically, SCAD yields the smallest median angle, which reflects the consistency of variable selection, and achieves a 100% elimination rate for the “cd” value, indicating the removal of irrelevant variables. The model’s prediction accuracy, as represented by the median of Err, is at its lowest point, while the estimation error of coefficient $\hat{\beta }$ is minimized. Elastic-Net, akin to Lasso and Adaptive-Lasso, exhibits limited capacity to effectively shrink non-relevant variables. (2) The performance metrics of Lasso, Elastic-Net, Adaptive-Lasso, and SCAD, which are four variable selection techniques, exhibit a decline with the increase in correlation “cor” among the predictive variables. However, SCAD continues to demonstrate superior performance among the four methods.

(3) The Angle and Err outliers presented in this study can serve as a means to delineate the stability of the model. Empirical evidence confirms the theoretical assertion that the SCAD model’s stability is compromised by the non-convexity of its penalty function and the instability of the local quadratic approximation algorithm.

(4) The metric used to evaluate the effectiveness of a model’s performance. The present study utilizes median values for Angle and Err, as opposed to mean values. This methodological choice serves to mitigate the influence of outliers on statistical analyses and to enhance the precision of the fitted model’s performance evaluation.

Table 4 presents the experimental outcomes of the numerical simulation comparison of two assessment techniques in high-dimensional settings. The analytical results presented in Table 4 are equivalent to those displayed in Table 3 when dealing with sizable sample sizes. Table 4 demonstrates a decline in the indicators of the models acquired through both evaluation methods. When the correlation coefficient (cor) among the predictor variables is significantly high, the performance indicators of each model tend to deteriorate rapidly. This could imply that the model has lost its ability to provide meaningful guidance. The present comparative analysis demonstrates that the outcomes of the two assessment techniques utilized for the four variable selection methods are in agreement. This finding substantiates the efficacy and suitability of the model evaluation method proposed in this study. Additionally, the stability attributes of the model fitted by the four variable selection methods are ascertained by scrutinizing the outliers of Angle and Err indices. In the statistical analysis of Angle and Err, which serve as performance indicators for model inference and prediction ability, the median value is utilized in lieu of the mean value to enhance the precision of the analysis outcomes.

Elastic-Net

The parsing path of the Elastic-Net penalty function and the cross-validation mean squared error graph are depicted in Fig. 6. The R programming language’s in-built ten-fold cross-validation function was utilized to determine the optimal adjustment parameter λ value of 0.107, which corresponds to the smallest training mean square error. Subsequently, an Elastic Net regression model was fitted using this value, and the resulting regression coefficients are presented in Table 6.

The model was validated on the test set and the mean square error of the index, which serves as a measure of the model’s accuracy, was calculated to be Err =0.926. The study aimed to perform variable selection and identify the factors that significantly affect systolic blood pressure, including tobacco consumption, adipose tissue concentration, obesity, alcohol consumption, and age. The quantitative relationship between these factors was established, with age and adipose tissue concentration found to have the most significant impact.

Lasso

The parsing path of the Lasso penalty function and the cross-validation mean squared error graph are depicted in Fig. 7. The R language’s in-built ten-fold cross-validation function was utilized to determine the optimal adjustment parameter λ = 0.0536, which yielded the smallest training mean square error. Subsequently, a Lasso regression model was fitted using this parameter, and the resulting regression coefficients are presented in Table 7.

The model should be validated on the test set and the mean square error of the index, which serves as an indicator of the model’s accuracy, should be calculated. The value of Err, which is equal to 0.929, should be obtained. Additionally, variable selection should be performed. In contrast to the Elastic-Net approach, the removal of the obesity term and the augmentation of the impact factor of adipose tissue concentration on systolic blood pressure not only streamlines the model, but also possesses a logical basis, as evidenced by the correlation coefficient of 0.72 obtained from the corr.test() function in the psych package when analyzing the relationship between the obesity term and adipose tissue concentration. The process of parsing exhibits a comparable pattern, and demonstrates a higher rate of compression as the coefficient experiences an increase in lambda.

SCAD

The parsing path of the SCAD penalty function is illustrated in Fig. 8. The region that is shaded in the diagram signifies the attainment of a local optimum by SCAD within that particular area, which consequently leads to the emergence of an unstable solution defect. Utilize the R programming language’s inherent cross-validation function, specifically the ten-fold cross-validation method, to derive the adjustment parameter λ = 0.149 at the point of minimal training mean square error. Subsequently, generate the SCAD fitting regression model. Table 8 displays the regression coefficients that correspond to the given data.

The model should be validated on the test set, and the mean square error of the index that reflects the model’s accuracy, Err =0.946, should be calculated. Additionally, variable selection should be performed. The study identified and isolated five key factors that exhibited a significant influence on systolic blood pressure. Additionally, the quantitative relationship between these factors was elucidated. Nevertheless, the solution provided by SCAD exhibits instability, and the assigned coefficient of ldl is negative, a result that lacks plausibility.

Adaptive Lasso

The parsing path of the Adaptive Lasso penalty function is illustrated in Fig. 9. The msgps package provides adjustment parameters for various variable selection criteria, including AIC, BIC, GCV, and Cp. The present study has opted for the GCV (ms.tuning = 2.86) approach to derive an Adaptive Lasso regression model with a good fit. Table 9 displays the regression coefficients that correspond to the given data.

The model is validated on the test set, and the mean square error of the index is calculated to reflect the accuracy of the model, resulting in an error of 0.92. This indicates that variable selection has been successfully achieved. Table 10 demonstrates the presented information. The concentration of adipose tissue, age, tobacco and alcohol consumption are identified as significant factors that exert a greater influence on systolic blood pressure. Further, the quantitative relationship between these factors is refined for better understanding.

Result analysis: (1) The Lasso technique is effective in achieving variable selection and identifying the factors that exert a significant influence on systolic blood pressure. In contrast to the Elastic-Net approach, the exclusion of the obesity variable and the augmentation of the determinants of adipose tissue concentration on systolic blood pressure are observed. The model can be deemed as straightforward and rational, as evidenced by the correlation coefficient of 0.72 between the obesity term and the adipose tissue concentration. The process of parsing exhibits a likeness as well, and experiences expedited compression as the coefficient ascends in tandem with lambda. (2) The instability of the solution provided by SCAD is evidenced by the negative coefficient assigned to ldl, which is deemed implausible. (3) Similar to Lasso, the Adaptive Lasso method effectively eliminated adipose tissue concentration, age, tobacco, and alcohol as variables that exert a significant influence on systolic blood pressure. Moreover, the Adaptive Lasso approach yielded the lowest mean square error.

The utilization of Elastic-Net, Lasso, Adaptive Lasso, and SCAD has been employed in fitting the random effect of the coefficient estimate $\hat{\beta }$ depicted in Fig. 10, as well as the random effect of the mean square error illustrated in Fig. 11, based on the given design. To examine and evaluate the stochastic impact of $\hat{\beta }$, it is possible to observe. The model developed through the utilization of Adaptive Lasso and SCAD exhibits a marginally superior capacity for variable selection in comparison to Elastic-Net and Lasso. The model coefficients are subject to significant fluctuations under the influence of random effects, whereas Elastic-Net and Lasso exhibit relatively stable behavior. Both the mean square error and stability of each model are similarly impacted by random factors. Thus, during the practical implementation phase, the most suitable technique for variable selection can be chosen to align with the model requirements.