- Introduction
Variable selection is important for high-dimensional data analysis in many research areas such as biology, signal processing, and collaborative filtering. For example, microarray experiments allow one to measure thousands of variables (genes, proteins) simultaneously. The data sets generated by these experiments are generally very large in terms of the number of predictors ( ) and often small in terms of the number of biological samples ( ). In regression analysis, this problem is often termed the “large and small problem” ( ) and presents a major barrier to traditional statistical methods.
With the development of computer and data collection technologies, the database sizes continue to grow and various statistical methodologies have been developed over the past several decades to cope with the challenges presented by these data. In particular, there are major challenges in parameter estimation, model, and variable selection. Several regression methods have been proposed for fitting multiple regression models, especially for the case when where the least-squares method could not be used.
Tibshirani proposed Lasso (Least Absolute Shrinkage and Selection Operator) [1], which minimizes the residual sum of squares subject to an -norm constraint. The Lasso penalty results in some coefficients being estimated to be completely zero, thus performing estimation and variable selection simultaneously. Following the seminal paper of Tibshirani [1], various extensions of the Lasso were developed, for example, the adaptive Lasso [2], Smoothly Clipped Absolute Deviation (SCAD) [3], etc.
Quantile regression, introduced by Koenker and Bassett [4], could be used when an estimate of the various quantiles (such as the median) of a conditional distribution is of interest. This allows one to see and compare how some quantiles of the response variable may be more affected by some predictor variables than other quantiles.
Some methods have combined regularized and robust regression methods to perform variable selection in high-dimensional data with outliers. For example, Rosset and Zhu [5] proposed the Huber Lasso method which combines Huber’s criterion loss with a Lasso penalty. The LAD-adaptive Lasso method is proposed by Wang et al. [6], combining the idea of Least Absolute Deviance (LAD) and adaptive Lasso. Lambert-Lacroix and Zwald [7] developed Huber’s Criterion with an adaptive Lasso which combines Huber’s loss function and adaptive Lasso penalty.
Fujisawa and Eguchi [8] introduce the gamma divergence for regression. It measures the difference between two conditional probability density functions. Arnold and Tibshirani [9] implemented the dual algorithm and their implementation is available in the R package genLasso. Taddy [10] proposed the gamma Lasso (GL) algorithm which can be seen as a computationally more attractive, multi-convex relaxation of best variable selection. Yi and Huang [11] developed an algorithm, called Semismooth Newton Coordinate Descent (SNCD), to obtain better efficiency and scalability for computing the solution paths of penalized quantile regression. Qin et .al [12] proposed a method called Maximum Tangent Likelihood Estimation (MTE). Christidis et al. [13] introduced the Split Regularized Regression (SRR) method which can be seen as a computationally more attractive, multi-convex relaxation of best-split selection. Zhu et .al [14] proposed Whitening Lasso (WLasso) to remove the correlations by applying a whitening transformation to the data before using the generalized Lasso criterion designed by Tibshirani and Taylor [15].
In the next section, we will give an overview of some group variable selection methods in linear regression.
2. Material and methods
We start from the standard model for multiple linear regression to describe the regression regularization methods. Let the data and the design matrix denoted by , . . , , the general linear model is usually written as
(1)
Here are the regression coefficients are the random errors, are the regressors for observation and The ordinary least squares (OLS) method estimates by minimizing the residual squared error, i.e. .
In general, OLS tends to give estimators with low biases but high variances, and better prediction accuracy can usually be obtained by lowering the variance with a little increased bias.
2.1 Lasso Regression
Tibshirani [1] proposed the Lasso penalty, a regularization technique for simultaneous estimation and variable selection for large data sets. The Lasso estimate is defined by:
(2)
An equivalent form of the Lasso is,
(3)
or
(4)
lambda is the parameter deciding the weight on minimizing the RSS compared to the penalty term is the sum of the absolute value of coefficients.
The Lasso minimizes the residual sum of squares subject to the sum of the absolute value of the coefficients being less than a constant. In other words, Lasso is a regression shrinkage method typically used in models with a large number of variables but relatively few observations. The main purpose of Lasso is to perform variable selection while fitting the regression line to the data. This is done by shrinking certain coefficients but in addition setting some of the coefficients also to zero. Lasso performs a regularization by adding a penalty to the objective under optimization. This penalty is the sum of the absolute value of coefficients and determines which coefficients to shrink and how much.
2.2 Adaptive Lasso
Zou [2] proposed a new version of the Lasso, which is called the adaptive Lasso. The penalized least squares with the adaptive Lasso are defined as
(5)
Instead of simply using the absolute value of the parameters as the penalization, adaptive weights are added for penalizing different coefficients differently. Zou (2006) suggested the use of stated weights, , where comes from minimizing the OLS or Lasso and is a user-chosen constant. The choice of is very important and Zou [2] suggested using OLS while can is chosen by K-fold cross-validation. The adaptive Lasso selects the true set of nonzero coefficients with probability tending to one.
2.3 Lad Lasso
Wang et al. [6] developed combined methods from Least Absolute Deviation (LAD) regression that is useful for robust regression, and also Lasso which is a popular choice for shrinkage estimation and variable selection, becoming Lad Lasso. The LAD-Lasso can be written as (Wang et al., [6]).
(6)
As can be seen, the LAD- criterion combines the LAD criterion and the Lasso penalty, and hence the resulting estimator is expected to be robust against outliers and also to enjoy a sparse representation.
2.4 Adaptive Lad Lasso
The Lad estimator is more robust than the OLS estimator, especially when datasets are subject to heavy-tailed errors or outliers. Lasso is a popular choice for shrinkage estimation. Adaptive Lad-Lasso is combining the two classical ideas to put forward a robust detection method to estimate change points in the mean-shift model. The basic idea is to convert the change point estimation problem into a variable selection problem with a penalty. The Adaptive Lad-Lasso can be written as (Lambert-Lacroix and Zwald, [7])
(7)
where is a known weights vector. In this model, the estimator is robust to outliers because the squared loss has been replaced by the -loss.
- Split Regularized Regression (SRR)
Christidis et al. [13] recently introduced the Split Regularized Regression (SRR) method which can be seen as a computationally more attractive, multi-convex relaxation of best-split selection. In high-dimensional regression, the proposed method builds an ensemble of models by splitting the set of covariates into different but possibly overlapping groups. A penalty term is introduced to encourage diversity between groups, and model stacking is used to generate accurate predictions.
The SRR-supervised variable clustering problem is to simultaneously estimate models, one for each cluster, then predict the target based on the average coefficient vector across clusters. SRR estimates multiple sparse coefficient vectors and encourages these vectors to be as diverse as possible.
Although SRR does not explicitly search for variable clusters, they may be inferred from the coefficient vectors: a variable belongs to a given cluster if its coefficient in the corresponding coefficient vector is nonzero. This approach can identify overlapping clusters and does not force coefficients of variables in the same cluster to tend towards the same value. The objective function of SRR is
(8)
The matrix is essentially a cluster membership matrix, where variable belongs to cluster if . Variables can belong to multiple clusters, but hard clusters can also be defined such that . Maximal diversity is achieved when the rows of contain only one nonzero element and thus each variable belongs to a single cluster (i.e. . The final vector of regression coefficients used for prediction is an average across all vectors :
.
Like many of the coefficient-grouping methods, cluster assignment is more of a side effect than a principal objective of SRR. Moreover, the best solutions in terms of prediction error tend to be complex models with extensive cluster overlap. As a result, variables are assigned to clusters with less certainty, and the clustering performance is less optimal.
2.6 DLasso
Haselimashhadi and Vinciotti [16] proposed a new penalty term that is capable of producing similar results to other well-known penalty functions in the context of regularized. The new penalty is differentiable and this penalty opens up the possibility of using it in many contexts where differentiability plays a key role. For example, a differentiable objective function could lead to more effcient implementations of parameter estimation procedures for certain models or to improved model selection criteria by a more accurate estimation of the bias term. The method is implemented in the R package DLASSO freely available from CRAN, http://CRAN.R-project.org/package=DLASSO.
- Simulation Study
In this section, we compare some regularized regression methods in low-dimensional with sparse and non-sparse coefficients , and high-dimensional with sparse coefficients ( , settings. For the sparse settings, we use a classical simulation setting, e.g. Yu et al. [17] and Li et al.[18], where . We draw the independent variables from a multivariate normal distribution, . The pairwise covariance between and is set to be . For the error , we choose a range of distributions in order to test the robustness of the methods to departures from normality. In particular, we consider the following cases , We design a mixture normal distribution with large outliers, similar to Lambert-Lacroix and Zwald [7], by drawing of the data from distribution and from a distribution, Laplace distribution, mixture of two Laplace distributions, t-distribution with 3 ) degrees of freedom and Gamma . Under all these cases, we compare the regularized regression methods described in the previous section, namely adaptive Lad Lasso (Xu and Ying, [19]; Lambert-Lacroix and Zwald, [7]), Split Regularized Regression (SRR) and DLasso (AIC, GIC, BIC, CGV). For the adaptive Lad Lasso we adapt some of the functions in the R package,for the SRR method, we use the R package and for the DLasso (AIC, GIC, BIC, CGV) methods, we use the R package .For the correlation , we experiment both with and . For the values we consider three cases:
- , which corresponds to the very sparse case with structures in the predictors.
- , which corresponds to the sparse case with structures in the predictors.
- for all , which corresponds to a dense case.
3.1 Simulation 1: low-dimensional with very sparse coefficients (Case 1)
In this section, we consider low-dimensional data with very sparse coefficients set with and . Table 1A, Table 1B, and Figure1 report the results of the simulation. We consider both the case of low correlation ( ) and that of high correlation ( ) of the predictors. The top panels report the median model error over iterations (similar results for the mean error), with the model error computed by , where are the estimated parameters and the sample covariance. The bottom panels report the true positives which are the number of correctly found non-zero coefficients. Here three correspond to the case of all non-zero coefficients being correctly detected.
Our results show that: the DLasso (GIC, BIC) methods do not perform well when the predictors are highly correlated; the adaptive Lad Lasso and the Split Regularized Regression (SRR) methods outperform all other methods for most error distributions.
Table1A: Average Median Model Error over 500 replications for the case: and values as in simulation1, Best method indicated in bold.
|
|
|
|
|
|
|
|
|
|
0.061
|
0.063
|
0.058
|
0.059
|
0.074
|
0.046
|
|
|
0.132
|
0.192
|
0.128
|
0.126
|
0.134
|
0.070
|
|
|
0.139
|
0.222
|
0.140
|
0.133
|
0.146
|
0.038
|
|
|
0.117
|
0.176
|
0.116
|
0.113
|
0.115
|
0.041
|
|
|
0.170
|
0.321
|
0.187
|
0.160
|
0.189
|
0.050
|
|
|
0.976
|
1.898
|
1.846
|
1.003
|
0.218
|
0.107
|
Table1B: Average Median Model Error over 500 replications for the case: and values as in simulation1, Best method indicated in bold.
|
|
|
|
|
|
|
|
|
|
0.084
|
0.089
|
0.081
|
0.081
|
0.059
|
0.056
|
|
|
0.157
|
0.217
|
0.159
|
0.149
|
0.106
|
0.122
|
|
|
0.185
|
0.230
|
0.191
|
0.179
|
0.121
|
0.055
|
|
|
0.154
|
0.191
|
0.141
|
0.137
|
0.097
|
0.052
|
|
|
0.252
|
0.314
|
0.259
|
0.233
|
0.172
|
0.089
|
|
|
1.363
|
1.956
|
1.912
|
1.730
|
0.179
|
0.163
|
Figure1: Comparison of regularized regression methods under different error distributions, for low (left) and high (right) correlated predictors. The top panels plot the median model error over 500 replications for simulation1 and the bottom panels the average true positives when and .
3.2 Simulation 2: high-dimensional with very sparse coefficients (Case 1)
We consider a similar setting to simulation 3.1 but with different sample size and several predictors. In particular, we consider a high-dimensional example with very sparse coefficients with and .Table 2A, Table 2B and Figure 2 report the median model error over replications, with the model error computed in the same way as in Figure1.
Table 2A: Average Median Model Error over 500 replications for the case: and values as in simulation 2, Best method indicated in bold.
|
|
|
|
|
|
|
|
|
|
0.946
|
0.274
|
0.252
|
0.306
|
0.291
|
0.081
|
|
|
1.717
|
0.799
|
0.382
|
0.378
|
0.454
|
0.134
|
|
|
1.928
|
1.063
|
0.573
|
0.528
|
0.490
|
0.085
|
|
|
1.530
|
0.684
|
0.416
|
0.474
|
0.448
|
0.069
|
|
|
2.281
|
1.319
|
0.732
|
0.563
|
0.628
|
0.123
|
|
|
2.799
|
2.634
|
2.413
|
2.031
|
0.923
|
2.039
|
Table2B: Average Median Model Error over 500 replications for the case: and values as in simulation 2, Best method indicated in bold.
|
|
|
|
|
|
|
|
|
|
0.628
|
0.234
|
0.184
|
0.203
|
0.178
|
0.335
|
|
|
1.373
|
0.849
|
0.370
|
0.329
|
0.330
|
0.299
|
|
|
1.447
|
1.007
|
0.410
|
0.338
|
0.363
|
0.141
|
|
|
1.087
|
0.668
|
0.287
|
0.269
|
0.297
|
0.096
|
|
|
2.004
|
1.411
|
0.673
|
0.450
|
0.448
|
0.277
|
|
|
2.822
|
2.662
|
2.523
|
1.972
|
0.526
|
0.421
|
Figure 2: Comparison of regularized regression methods under different error distributions, for low (left) and high (right) correlated predictors. The top panels plot the median model error over 500 replications for simulation 2 and the bottom panels the average true positives when and .
The results support the performance of the methods: DLasso (AIC, GIC) do not perform well when the predictors are highly correlated, and the adaptive Lad Lasso method outperforms all other methods as departures from normality increase. This is particularly evident in the case and .
3.3 Simulation 3: low- dimensional with non-sparse coefficients (Case 2)
To investigate the performance of variable selection methods, we set up a new simulation where we have as in case 2, that is a sparse situation. Table 3A, Table 3B and Figure 3 report the median model error over replications for the cases and
Table3.A: Average Median Model Error over 500 replications for the case: and values as in simulation 3, Best method indicated in bold.
|
|
|
|
|
|
|
|
|
|
0.260
|
0.228
|
0.221
|
0.222
|
0.159
|
2.788
|
|
|
0.388
|
0.785
|
0.355
|
0.360
|
0.291
|
2.260
|
|
|
0.400
|
0.851
|
0.366
|
0.365
|
0.302
|
2.776
|
|
|
0.389
|
0.632
|
0.339
|
0.347
|
0.250
|
2.562
|
|
|
0.521
|
1.148
|
0.523
|
0.472
|
0.409
|
2.429
|
|
|
1.537
|
3.618
|
3.566
|
1.657
|
0.450
|
2.400
|
Table3B: Average Median Model Error over 500 replications for the case: and values as in simulation3, Best method indicated in bold.
|
|
|
|
|
|
|
|
|
|
0.757
|
0.428
|
0.460
|
0.513
|
0.117
|
0.366
|
|
|
0.968
|
0.819
|
0.744
|
0.826
|
0.203
|
0.712
|
|
|
1.091
|
0.865
|
0.787
|
0.847
|
0.209
|
0.496
|
|
|
1.036
|
0.741
|
0.684
|
0.738
|
0.188
|
0.488
|
|
|
1.133
|
1.107
|
0.977
|
0.933
|
0.275
|
0.426
|
|
|
2.944
|
3.564
|
3.526
|
2.678
|
0.352
|
0.610
|
Figure 3: Comparison of regularized regression methods under different error distributions, for low (left) and high (right) correlated predictors. The top panels plot the median model error over 500 replications for simulation 3 and the bottom panels the average true positives when and .
From the results in Table 3A, Table 3B and Figure 3, our simulation study confirms that the SRR outperforms all other methods as departures from normality increase. This is particularly evident in the case when the predictors are highly correlated.
3.4 Simulation 4: high-dimensional with non-sparse coefficients (Case 2)
To investigate the performance of variable selection methods, we set up a new simulation where we have as in case 2, that is a sparse situation. Table 4A, Table 4B and Figure 4 report the median model error over replications for the cases and .
Table 4A: Average Median Model Error over 500 replications for the case: and values as in simulation 4, Best method indicated in bold.
|
|
|
|
|
|
|
|
|
|
0.942
|
0.410
|
0.417
|
0.507
|
0.452
|
1.515
|
|
|
1.739
|
1.000
|
0.731
|
0.826
|
0.782
|
2.378
|
|
|
1.893
|
1.107
|
0.793
|
0.863
|
0.775
|
2.108
|
|
|
1.607
|
0.892
|
0.693
|
0.830
|
0.729
|
1.648
|
|
|
2.282
|
1.466
|
0.953
|
0.786
|
0.839
|
2.549
|
|
|
2.864
|
2.694
|
2.578
|
3.347
|
1.174
|
2.328
|
Table 4B: Average Median Model Error over 500 replications for the case: and values as in simulation 4, Best method indicated in bold.
|
|
|
|
|
|
|
|
|
|
0.701
|
0.264
|
0.236
|
0.247
|
0.237
|
0.459
|
|
|
1.428
|
0.901
|
0.438
|
0.409
|
0.388
|
0.542
|
|
|
1.514
|
0.919
|
0.491
|
0.449
|
0.424
|
0.772
|
|
|
1.298
|
0.749
|
0.428
|
0.410
|
0.379
|
0.666
|
|
|
2.038
|
1.324
|
0.642
|
0.500
|
0.496
|
0.469
|
|
|
2.852
|
2.799
|
2.742
|
1.949
|
0.658
|
0.543
|
Figure 4: Comparison of regularized regression methods under different error distributions, for low (left) and high (right) correlated predictors. The top panels plot the median model error over 500 replications for simulation 4 and the bottom panels the average true positives when and .
From the results in Table 4A and Table 4B and Figure 4 our simulation study confirms the performances of the adaptive Lad Lasso and the SRR methods outperform all other methods for most error distributions, Furthermore, the results show how DLasso (AIC) is the worst performing method in the case of departure from normality especially when the predictors are highly correlated.
3.5 Simulation 5: low- dimensional with non-sparse coefficients (Case 3)
To investigate the performance of variable selection methods, we set up a new simulation where we have as in case 3, that is a non-sparse situation. Table5A, Table 5B and Figure 5 report the median model error over replications for the cases and .
Table5.A: Average Median Model Error over 500 replications for the case: and values as in simulation 5, Best method indicated in bold.
|
|
|
|
|
|
|
|
|
|
0.387
|
0.139
|
0.304
|
0.386
|
0.085
|
0.219
|
|
|
0.456
|
0.230
|
0.292
|
0.456
|
0.124
|
0.256
|
|
|
0.446
|
0.245
|
0.333
|
0.446
|
0.132
|
0.243
|
|
|
0.430
|
0.217
|
0.287
|
0.429
|
0.122
|
0.249
|
|
|
0.377
|
0.310
|
0.344
|
0.377
|
0.163
|
0.221
|
|
|
0.428
|
1.955
|
1.891
|
0.371
|
0.198
|
0.269
|
Table5B: Average Median Model Error over 500 replications for the case: and values as in simulation5, Best method indicated in bold.
|
|
|
|
|
|
|
|
|
|
0.922
|
0.123
|
0.131
|
0.149
|
0.040
|
0.094
|
|
|
1.027
|
0.216
|
0.215
|
0.238
|
0.066
|
0.125
|
|
|
0.740
|
0.243
|
0.246
|
0.316
|
0.067
|
0.091
|
|
|
0.976
|
0.197
|
0.202
|
0.224
|
0.059
|
0.090
|
|
|
0.768
|
0.307
|
0.308
|
0.374
|
0.083
|
0.099
|
|
|
0.564
|
1.064
|
0.917
|
0.570
|
0.095
|
0.141
|
Figure 5: Comparison of regularized regression methods under different error distributions, for low (left) and high (right) correlated predictors. The median model error is plotted over 500 replications for simulation 5 when and .
From the results in Figure 5, Table 5A and Table 5B our simulation study confirms that the SRR outperforms all other methods as departures from normality increase. This is particularly evident in the case when the predictors are highly correlated.
3.6 Simulation 6: high-dimensional with non-sparse coefficients (Case 3)
To investigate the performance of variable selection methods, we set up a new simulation where we have as in case 3, that is a non-sparse situation. Table 6A, Table 6B and Figure 6 report the median model error over replications for the cases and .
Table 6.A: Average Median Model Error over 500 replications for the case: and values as in simulation 6, Best method indicated in bold.
|
|
|
|
|
|
|
|
|
|
0.949
|
0.686
|
0.828
|
1.780
|
0.779
|
1.675
|
|
|
1.720
|
1.160
|
1.082
|
2.043
|
1.170
|
2.220
|
|
|
1.827
|
1.229
|
1.041
|
1.677
|
1.175
|
1.848
|
|
|
1.610
|
1.077
|
0.977
|
1.794
|
1.265
|
1.723
|
|
|
2.304
|
1.534
|
1.217
|
1.536
|
1.234
|
1.584
|
|
|
2.838
|
2.427
|
2.220
|
1.982
|
1.732
|
2.187
|
Table 6B: Average Median Model Error over 500 replications for the case: and values as in simulation 6, Best method indicated in bold.
|
|
|
|
|
|
|
|
|
|
3.211
|
0.512
|
0.569
|
0.794
|
0.428
|
7.648
|
|
|
1.983
|
0.863
|
0.752
|
0.878
|
0.674
|
9.493
|
|
|
2.289
|
0.956
|
0.816
|
0.942
|
0.703
|
10.072
|
|
|
2.050
|
0.774
|
0.760
|
0.879
|
0.633
|
11.435
|
|
|
2.749
|
1.399
|
0.962
|
1.059
|
0.844
|
7.312
|
|
|
2.745
|
2.676
|
2.591
|
2.264
|
1.022
|
7.942
|
Figure 6: Comparison of regularized regression methods under different error distributions, for low (left) and high (right) correlated predictors. The median model error is plotted over 500 replications for simulation 6 when and .
From the results in Table 6A and Table 6B and Figure 6 our simulation study confirms the performances of the SRR method outperform all other methods as departures from normality increase, Furthermore, the results show how adaptive Lad Lasso is the worst performing method in the case of departure from normality especially when the predictors are highly correlated.