# Stochastic Restricted LASSO-Type Estimator in the Linear Regression Model

Authors: Manickavasagar Kayanan; Pushpakanthie Wijekoon
Journal: Journal of Probability and Statistics
 (2020)

---

## Abstract

Among several variable selection methods, LASSO is the most desirable estimation procedure for handling regularization and variable selection simultaneously in the high-dimensional linear regression models when multicollinearity exists among the predictor variables. Since LASSO is unstable under high multicollinearity, the elastic-net (Enet) estimator has been used to overcome this issue. According to the literature, the estimation of regression parameters can be improved by adding prior information about regression coefficients to the model, which is available in the form of exact or stochastic linear restrictions. In this article, we proposed a stochastic restricted LASSO-type estimator (SRLASSO) by incorporating stochastic linear restrictions. Furthermore, we compared the performance of SRLASSO with LASSO and Enet in root mean square error (RMSE) criterion and mean absolute prediction error (MAPE) criterion based on a Monte Carlo simulation study. Finally, a real-world example was used to demonstrate the performance of SRLASSO.


---

## Body

## 1. Introduction

Let us consider the linear regression model(1)y=Xβ+ϵ,where y is the n × 1 vector of observations on the dependent variable, X is the n × p matrix of observations on the nonstochastic predictor variables, β is a p × 1 vectors of unknown coefficients, and ϵ is the n × 1 vector of random error terms, which is independent and identically normally distributed with the mean zero and common variance σ2, that is, E(ϵ) = 0 and E(ϵϵ′) = Ω = σ2I.It is well-known that ordinary least square estimator (OLSE) is a best linear unbiased estimator for model (1), and it is defined as(2)β^OLSE=X′X−1X′y.Furthermore, researchers [1, 2] have shown that parameter estimation is improved when prior information for regression coefficients is available, which can be in the form of exact linear restrictions or stochastic linear restrictions. Let us assume that there exists prior information on β in the form of stochastic linear restriction as(3)ϕ=Rβ+v,where ϕ is the q × 1 vector, R is the q × p matrix with rank q, and v is the q × 1 vector of disturbances, such that E(v) = 0, D(v) = E(vv′) = Ψ = σ2W (W is positive definite) and E(vϵ′) = 0. Note that equation (3) will be the exact linear restriction when v = 0.Theil and Goldberger [2] proposed the mixed regression estimator (MRE) by combining the models (1) and (3), and it is defined as(4)β^MRE=X′X+R′W−1R−1X′y+R′W−1ϕ.According to the literature, it was noted that OLSE and MRE are unstable when the number of predictors is high. In this case, the variable selection methods such as forward selection, backward selection, and step-wise selection have been used. However, these methods are also unstable when multicollinearity exists among predictor variables. As a remedial solution to this problem, Tibshirani [3] proposed the least absolute shrinkage and selection operator (LASSO) by considering the model (1) to handle both multicollinearity and variable selection simultaneously in the high-dimensional linear regression model. The LASSO estimator is defined as(5)β^LASSO=argminβy−Xβ′y−Xβsubject to∑j=1pβj≤t,where t ≥ 0 is a turning parameter. The solutions of LASSO can be obtained using either the standard quadratic programming technique or least angle regression (LARS) [4] algorithm. According to Zou et al. [5, 6], LASSO is unstable when high multicollinearity exists among predictor variables. Therefore, they proposed elastic-net (Enet) estimator as an alternative to LASSO to handle this issue. The Enet estimator is defined as(6)β^Enet=argminβy−Xβ′y−Xβ+k∑j=1pβj2subject to∑j=1pβj≤t.The Enet solutions can be obtained using LARS-EN algorithm, which is the modified version of LARS algorithm.Norouzirad et al. [7] and Tuaç and Arslan [8] attempted to combine LASSO with the exact linear restriction, and their work did not well define how the exact restriction is incorporated since there is no analytical solution for LASSO. In this article, we proposed stochastic restricted LASSO-type estimator (SRLASSO) by combining LASSO and stochastic restrictions. Furthermore, we compared the performance of SRLASSO with LASSO and Enet in root mean square error (RMSE) criterion and mean absolute prediction error (MAPE) criterion using a Monte Carlo simulation study and a real-world example. The structure of the rest of the article is as follows: Section 2 describes SRLASSO and algorithm to find SRLASSO solutions, Section 3 shows the performance of SRLASSO, and Section 4 concludes the article and references provided at the end of the paper.


## 2. Stochastic Restricted LASSO-Type Estimator (SRLASSO)

By considering equation (3) as an additional constrain, we define stochastic restricted LASSO-type estimator (SRLASSO) as(7)β^SRLASSO=argminβy−Xβ′y−Xβsubject to∑j=1pβj≤t,Rβ=ϕ−v.We can view this as a quadratic optimization problem. Here, we have 2p + q constraints. However, this method is not suitable in a practical situation if p is large. So, we propose a stochastic restricted LARS (SRLARS) algorithm, which is the modified version of the LARS algorithm, to find SRLASSO solutions. In SRLARS, we consolidate MRE with LARS.


### 2.1. Stochastic Restricted LARS (SRLARS)

Standardize the predictor variablesX to have a mean of zero and a standard deviation of one, and response variable y to have a mean zero. Let residuals of the models (1) and (3) be r=ϵ^ and τ=v^, respectively.Step 1.
Start withβ^=0, r0 = y, and τ0 = ϕ.Step 2.
Find the predictorXj1 most correlated with r0 as follows:(8)Xj1=maxjCorXj,r0;j=1,2,…,p.
Letβ^ji be the regression coefficient of Xji. Then, increase the estimate of β^j1 from 0 toward MRE until any other predictor Xj2 has a high correlation with the current residual as Xj1 does. At this point, SRLARS proceeds in the equiangular direction between the two predictors Xj1 and Xj2 instead of continuing in the direction based on Xj1.
In a similar way,ith variable Xji eventually earns its way into the active set, and then SRLARS proceeds in the equiangular direction between Xj1, Xj2, …, Xji. Continue adding variables to the active set in this way moving in the direction defined by the least angle direction. In the intermediate steps, the coefficient estimates are updated using the following formula:(9)β^ji=β^ji−1+αiui,where αi is a value between 0 and 1 which represents how far the estimate moves in the direction before another variable enters the model and the direction changes again, and ui is the equiangular vector.
The directionui is calculated using the following formula based on MRE:(10)ui=EiEi′X′X+R′W−1REi−1Ei′X′ri−1+R′W−1τi−1,where Ei is the matrix with column (ej1, ej2, …, eji) and ej be the jth standard unit vector in ℝp, which has the index of variables selected in each subsequent step.
Then,αi is calculated as follows:(11)αi=minαji+,αji−,αji∗∈0,1,where(12)αji±=Corri−1,Xji±Corri−1,XjCorri−1,Xji±CorXui,Xj,for any j such that β^ji−1=0 and(13)αji∗=−β^ji−1ui,
for anyj such that β^ji−1≠0.Step 3.
Ifαi=αji∗, then Ei is the matrix formed by removing the column ej from Ei−1. Then ri and τi related to the current step is calculated as(14)ri=ri−1−αiXui,τi=τi−1−αiRui,and then move to the next step where ji+1 is the value of j such that αi=αji+ or αi=αji− or αji∗.Step 4.
Proceed Step2 until αi = 1.


### 2.2. Properties of SRLARS

SRLARS algorithm sequentially updates the SRLASSO estimates. It requiresO(m3 + pm2) operations, where m is the number of steps. The prediction performance of SRLARS is evaluated using the RMSE criterion and MAPE criterion, which is described in Section 3. According to Efron et al. [4], the conventional tuning parameter is t=∑j=1pβ^j, and the suitable value of turning parameter t for the particular problem is selected using K-fold cross-validation.


### 2.3. Selection of Prior Information

According to Nagar and Kakwani [9], we can define the prior information as follows: Let β1 be a vector of some selected q elements of β and β2 be the rest of elements. Assume that b is the known unbiased estimates of β1. By using the “two sigma rule,” now we can write the range of β1 as b ± 2SE(b). Based on that, we can set the expressions of equation (3) as r^=b, R^=10…0…001…0…0..…...00…1…0q×p, β=β1β2, and σ2W^=SEb10…00SEb2…0..….00…SEblq×q.


## 2.1. Stochastic Restricted LARS (SRLARS)

Standardize the predictor variablesX to have a mean of zero and a standard deviation of one, and response variable y to have a mean zero. Let residuals of the models (1) and (3) be r=ϵ^ and τ=v^, respectively.Step 1.
Start withβ^=0, r0 = y, and τ0 = ϕ.Step 2.
Find the predictorXj1 most correlated with r0 as follows:(8)Xj1=maxjCorXj,r0;j=1,2,…,p.
Letβ^ji be the regression coefficient of Xji. Then, increase the estimate of β^j1 from 0 toward MRE until any other predictor Xj2 has a high correlation with the current residual as Xj1 does. At this point, SRLARS proceeds in the equiangular direction between the two predictors Xj1 and Xj2 instead of continuing in the direction based on Xj1.
In a similar way,ith variable Xji eventually earns its way into the active set, and then SRLARS proceeds in the equiangular direction between Xj1, Xj2, …, Xji. Continue adding variables to the active set in this way moving in the direction defined by the least angle direction. In the intermediate steps, the coefficient estimates are updated using the following formula:(9)β^ji=β^ji−1+αiui,where αi is a value between 0 and 1 which represents how far the estimate moves in the direction before another variable enters the model and the direction changes again, and ui is the equiangular vector.
The directionui is calculated using the following formula based on MRE:(10)ui=EiEi′X′X+R′W−1REi−1Ei′X′ri−1+R′W−1τi−1,where Ei is the matrix with column (ej1, ej2, …, eji) and ej be the jth standard unit vector in ℝp, which has the index of variables selected in each subsequent step.
Then,αi is calculated as follows:(11)αi=minαji+,αji−,αji∗∈0,1,where(12)αji±=Corri−1,Xji±Corri−1,XjCorri−1,Xji±CorXui,Xj,for any j such that β^ji−1=0 and(13)αji∗=−β^ji−1ui,
for anyj such that β^ji−1≠0.Step 3.
Ifαi=αji∗, then Ei is the matrix formed by removing the column ej from Ei−1. Then ri and τi related to the current step is calculated as(14)ri=ri−1−αiXui,τi=τi−1−αiRui,and then move to the next step where ji+1 is the value of j such that αi=αji+ or αi=αji− or αji∗.Step 4.
Proceed Step2 until αi = 1.


## 2.2. Properties of SRLARS

SRLARS algorithm sequentially updates the SRLASSO estimates. It requiresO(m3 + pm2) operations, where m is the number of steps. The prediction performance of SRLARS is evaluated using the RMSE criterion and MAPE criterion, which is described in Section 3. According to Efron et al. [4], the conventional tuning parameter is t=∑j=1pβ^j, and the suitable value of turning parameter t for the particular problem is selected using K-fold cross-validation.


## 2.3. Selection of Prior Information

According to Nagar and Kakwani [9], we can define the prior information as follows: Let β1 be a vector of some selected q elements of β and β2 be the rest of elements. Assume that b is the known unbiased estimates of β1. By using the “two sigma rule,” now we can write the range of β1 as b ± 2SE(b). Based on that, we can set the expressions of equation (3) as r^=b, R^=10…0…001…0…0..…...00…1…0q×p, β=β1β2, and σ2W^=SEb10…00SEb2…0..….00…SEblq×q.


## 3. Performance of SRLASSO

SRLASSO is compared with LASSO and Enet using the RMSE criterion and MAPE criterion, which are the expected prediction errors of the algorithms, and are defined as(15)RMSEβ^=1nynew−Xnewβ^′ynew−Xnewβ^,MAPEβ^=1n∑ynew−Xnewβ^,where (ynew, Xnew) denotes the new data which are not used to obtain the parameter estimates, n is the number of new observations, and β^ is the estimated value of β using the respective algorithm. A Monte Carlo simulation study and a real-world example are used for the comparison.


### 3.1. Simulation Study

According to McDonald and Galarneau [10], first we generate the predictor variables by using the following formula:(16)xi,j=1−ρ2zi,j+ρzi,m+1;i=1,2,…,n,j=1,2,…,m,where zi,j is an independent standard normal pseudo random number and ρ is the theoretical correlation between any two explanatory variables.In this study, we have used a linear regression model of 100 observations and 20 predictors. A dependent variable is generated by using the following equation:(17)yi=β1xi,1+β2xi,2+⋯+β5xi,20+ϵi;i=1,2,…,100,where ϵi is a normal pseudo random number with a mean zero and common variance σ2.We chooseβ = (β1, β2, …, β20) as the normalized eigenvector corresponding to the largest eigenvalue of X′X for which β′β = 1. To define the prior information according to Section 2.3, we assume that OLSE estimates of the first four elements of β are unbiased, which are the estimates of b. To investigate the effects of different degrees of multicollinearity on the estimators, we choose ρ = (0.5, 0.7, 0.9), which represents weak, moderated, and high multicollinearity. For the analysis, we have simulated 50 data sets consisting of 50 observations to fit the model and 50 observations to calculate the RMSE and MAPE. The cross-validated RMSE and MAPE of the estimators are displayed in Figure 1 and Figure 2, respectively. The median cross-validated RMSE and MAPE of the estimators are displayed in Table 1.Figure 1
Cross-validated RMSE values of the estimators when (a)ρ = 0.5, (b) ρ = 0.7, and (c) ρ = 0.9.
(a)
(b)
(c)Figure 2
Cross-validated MAPE values of the estimators when (a)ρ = 0.5, (b) ρ = 0.7, and (c) ρ = 0.9.
(a)
(b)
(c)Table 1
Median cross-validated RMSE and MAPE values.
Estimators
RMSE
MAPE
t
k
Selected variables
ρ = 0.5
LASSO
3.440
2.830
7.60
—
16
Enet
3.454
2.854
7.14
0.01
17
SRLASSO
3.280
2.63
6.16
—
15
ρ = 0.7
LASSO
3.512
2.804
8.66
—
15
Enet
3.609
2.879
8.95
0.01
17
SRLASSO
3.322
2.676
6.43
—
15
ρ = 0.9
LASSO
3.502
2.769
13.48
—
15
Enet
3.466
2.814
14.47
0.01
16
SRLASSO
3.318
2.732
11.66
—
15From Figures1 and 2 and Table 1, we can observe that SRLASSO always shows better performance compared to LASSO and Enet in both RMSE criterion and MAPE criterion under all degrees of multicollinearity.


### 3.2. Real-World Example

As a numerical example, the well-known Prostate Cancer Data [11] was used to compare the performance of SRLASSO. This data set is attached with “lasso2” R package. In the Prostate Cancer Data, the predictors are the following eight clinical measures: log cancer volume (lcavol), log prostate weight (lweight), age, log of the amount of benign prostatic hyperplasia (lbph), seminal vesicle invasion (svi), log capsular penetration (lcp), Gleason score (gleason), and percentage Gleason score 4 or 5 (pgg45). The response is the log of prostate specific antigen (lpsa), and the dataset has 97 observations. The variance inflation factor (VIF) values of the predictor variables of the dataset are 3.09, 2.97, 2.47, 2.05, 1.95, 1.37, 1.36, and 1.32, and the condition number is 243, which shows evidence of multicollinearity among the predictor variables. Stamey et al. [11] have examined the correlation between the level of prostate specific antigen with those eight clinical measures. Furthermore, Tibshirani [3] and Tibshirani et al. [4] have used this data to examine the performance of LASSO and LARS algorithm, respectively. We have used 67 observations to fit the model and 30 observations to calculate RMSE and MAPE. We assume that OLSE estimates of the first three regression coefficients of Prostate Cancer Data are unbiased and we defined the prior information for this data based on Section 2.3. The cross-validated RMSE and MAPE of the estimators are displayed in Table 2, and coefficient paths of each estimator are displayed in Figure 3.Table 2
Cross-validated RMSE and MAPE values of Prostate Cancer Data.
Estimators
RMSE
MAPE
t
k
LASSO
0.8567
0.6496
1.2263
—
Enet
0.8063
0.6010
1.2897
0.01
SRLASSO
0.7833
0.5803
1.0415
—Figure 3
Coefficient paths of the (a) LASSO, (b) Enet, and (c) SRLASSO versust=∑j=1pβ^j for the Prostate Cancer Data.
(a)
(b)
(c)From Table2, we can observe that SRLASSO outperforms LASSO and Enet on Prostate Cancer Data in both RMSE criterion and MAPE criterion. Furthermore, we can note that the selection of variables is different for each estimator by comparing Figures 3(a)–3(c).


## 3.1. Simulation Study

According to McDonald and Galarneau [10], first we generate the predictor variables by using the following formula:(16)xi,j=1−ρ2zi,j+ρzi,m+1;i=1,2,…,n,j=1,2,…,m,where zi,j is an independent standard normal pseudo random number and ρ is the theoretical correlation between any two explanatory variables.In this study, we have used a linear regression model of 100 observations and 20 predictors. A dependent variable is generated by using the following equation:(17)yi=β1xi,1+β2xi,2+⋯+β5xi,20+ϵi;i=1,2,…,100,where ϵi is a normal pseudo random number with a mean zero and common variance σ2.We chooseβ = (β1, β2, …, β20) as the normalized eigenvector corresponding to the largest eigenvalue of X′X for which β′β = 1. To define the prior information according to Section 2.3, we assume that OLSE estimates of the first four elements of β are unbiased, which are the estimates of b. To investigate the effects of different degrees of multicollinearity on the estimators, we choose ρ = (0.5, 0.7, 0.9), which represents weak, moderated, and high multicollinearity. For the analysis, we have simulated 50 data sets consisting of 50 observations to fit the model and 50 observations to calculate the RMSE and MAPE. The cross-validated RMSE and MAPE of the estimators are displayed in Figure 1 and Figure 2, respectively. The median cross-validated RMSE and MAPE of the estimators are displayed in Table 1.Figure 1
Cross-validated RMSE values of the estimators when (a)ρ = 0.5, (b) ρ = 0.7, and (c) ρ = 0.9.
(a)
(b)
(c)Figure 2
Cross-validated MAPE values of the estimators when (a)ρ = 0.5, (b) ρ = 0.7, and (c) ρ = 0.9.
(a)
(b)
(c)Table 1
Median cross-validated RMSE and MAPE values.
Estimators
RMSE
MAPE
t
k
Selected variables
ρ = 0.5
LASSO
3.440
2.830
7.60
—
16
Enet
3.454
2.854
7.14
0.01
17
SRLASSO
3.280
2.63
6.16
—
15
ρ = 0.7
LASSO
3.512
2.804
8.66
—
15
Enet
3.609
2.879
8.95
0.01
17
SRLASSO
3.322
2.676
6.43
—
15
ρ = 0.9
LASSO
3.502
2.769
13.48
—
15
Enet
3.466
2.814
14.47
0.01
16
SRLASSO
3.318
2.732
11.66
—
15From Figures1 and 2 and Table 1, we can observe that SRLASSO always shows better performance compared to LASSO and Enet in both RMSE criterion and MAPE criterion under all degrees of multicollinearity.


## 3.2. Real-World Example

As a numerical example, the well-known Prostate Cancer Data [11] was used to compare the performance of SRLASSO. This data set is attached with “lasso2” R package. In the Prostate Cancer Data, the predictors are the following eight clinical measures: log cancer volume (lcavol), log prostate weight (lweight), age, log of the amount of benign prostatic hyperplasia (lbph), seminal vesicle invasion (svi), log capsular penetration (lcp), Gleason score (gleason), and percentage Gleason score 4 or 5 (pgg45). The response is the log of prostate specific antigen (lpsa), and the dataset has 97 observations. The variance inflation factor (VIF) values of the predictor variables of the dataset are 3.09, 2.97, 2.47, 2.05, 1.95, 1.37, 1.36, and 1.32, and the condition number is 243, which shows evidence of multicollinearity among the predictor variables. Stamey et al. [11] have examined the correlation between the level of prostate specific antigen with those eight clinical measures. Furthermore, Tibshirani [3] and Tibshirani et al. [4] have used this data to examine the performance of LASSO and LARS algorithm, respectively. We have used 67 observations to fit the model and 30 observations to calculate RMSE and MAPE. We assume that OLSE estimates of the first three regression coefficients of Prostate Cancer Data are unbiased and we defined the prior information for this data based on Section 2.3. The cross-validated RMSE and MAPE of the estimators are displayed in Table 2, and coefficient paths of each estimator are displayed in Figure 3.Table 2
Cross-validated RMSE and MAPE values of Prostate Cancer Data.
Estimators
RMSE
MAPE
t
k
LASSO
0.8567
0.6496
1.2263
—
Enet
0.8063
0.6010
1.2897
0.01
SRLASSO
0.7833
0.5803
1.0415
—Figure 3
Coefficient paths of the (a) LASSO, (b) Enet, and (c) SRLASSO versust=∑j=1pβ^j for the Prostate Cancer Data.
(a)
(b)
(c)From Table2, we can observe that SRLASSO outperforms LASSO and Enet on Prostate Cancer Data in both RMSE criterion and MAPE criterion. Furthermore, we can note that the selection of variables is different for each estimator by comparing Figures 3(a)–3(c).


## 4. Conclusions

This study clearly showed that SRLASSO does a better performance than LASSO and Enet in both RMSE criterion and MAPE criterion when multicollinearity exists among the predictor variables. Therefore, SRLASSO can be used as an alternative estimator of LASSO and Enet if prior information is accessible on the regression coefficients. The proposed SRLARS algorithm can be used to obtain SRLASSO solutions.



---

*Source: 7352097-2020-03-30.xml*