Lasso 模型

根据Lasso模型，实现了lasso回归的R函数代码，使用了Bostan2与Boston3数据集，最后将lambda.min与lambda.lse的Lasso回归、经过lasso筛选features的重新线性回归、线性回归、lambda.min与lambda.lse的岭回归和PCR模型进行了比较

Part I: Lasso Implementation

Firstly, load the data.

df <- read.table("https://liangfgithub.github.io/F21/Coding2_myData.csv", 
                 header = TRUE,
                 sep = ",")

And below is the one_var_lasso function

one_var_lasso = function(r, x, lam) {
  xx = sum(x^2)
  xr = sum(r * x)
  b = (abs(xr) - lam/2)/xx
  b = sign(xr) * ifelse(b > 0, b, 0)
  return(b)
}

Below is MyLasso function

MyLasso = function(X, y, lam.seq, maxit = 500) {
  
  # X: n-by-p design matrix without the intercept 
  # y: n-by-1 response vector 
  # lam.seq: sequence of lambda values (arranged from large to small)
  # maxit: number of updates for each lambda 
  # Center/Scale X
  # Center y
  
  n = length(y)
  p = dim(X)[2]
  nlam = length(lam.seq)
  lam.seq = sort(lam.seq , decreasing = TRUE)
  
  
  ##############################
  # YOUR CODE: 
  # Record the corresponding means and scales
  # For example, 
  y.mean = mean(y)
  X.mean= colMeans(X)
  #print(X.mean)
  #X.std=sqrt(colMeans(t(t(X)-X.mean)^2))
  
  X.std=apply((t(t(X)-X.mean)^2),2,sum)
  X.std=sqrt(X.std/(n-1))

  Xs = scale(X, center = TRUE, scale = TRUE)
  ##############################
  
  # Initilize coef vector b and residual vector r
  b = rep(0, p)
  r = y
  B = matrix(nrow = nlam, ncol = p + 1)
  
  # Triple nested loop
  for (m in 1:nlam) {
    lam = 2 * n * lam.seq[m]
    for (step in 1:maxit) {
      for (j in 1:p) {
        r = r + (Xs[, j] * b[j])
        b[j] = one_var_lasso(r, Xs[, j], lam)
        r = r - Xs[, j] * b[j]
      }
    }
    B[m, ] = c(0, b/X.std)
    
    B[m,1]=y.mean-b %*% (X.mean/X.std)
    
    #B[m,1]=y.mean-sum((X.mean)*b)
  }
  
  return(t(B))
}

X = as.matrix(df[, -14])
y = df$Y

lam.seq = exp(seq(-1, -8, length.out = 80))
myout = MyLasso(X, y, lam.seq, maxit = 100) 
rownames(myout) = c("Intercept", colnames(X)) 

x.index = log(lam.seq)
beta = myout[-1, ]  # beta is a 13-by-80 matrix
matplot(x.index, t(beta),
        xlim = c(min(x.index), max(x.index)),
        lty = 1,
        xlab = "Log Lambda",
        ylab = "Coefficients",
        type="l", 
        lwd = 1)

library(glmnet)

## 载入需要的程辑包：Matrix

## Loaded glmnet 4.1-2

lasso.fit = glmnet(X, y, alpha = 1, lambda = lam.seq)
# coef(lasso.fit)
write.csv(as.matrix(coef(lasso.fit)), file = "Coding2_lasso_coefs.csv", 
          row.names = FALSE)

So, if I check the maximum of the difference value of mylasso output and lasso.fit output

max(abs(coef(lasso.fit) - myout))

## [1] 0.00456362

We can look the plot from glmnet

plot(lasso.fit, xvar = "lambda")

Part II: Simulation Study of BostonData2

library(glmnet) 
library(pls)

## 
## 载入程辑包：'pls'

## The following object is masked from 'package:stats':
## 
##     loadings

set.seed(4358) 

myData = read.csv("https://liangfgithub.github.io/F21/BostonData2.csv")
myData = myData[, -1]

Repeat the following simulation 50 times: In each iteration, randomly split the data into two parts, 75% for training and 25% for testing. So, firstly, we should split the data set.

X = data.matrix(myData[,-1])  
Y = myData[,1]  

T = 50
n = length(Y)
ntest = round(n * 0.25)  # test set size
ntrain = n - ntest  # training set size
all.test.id = matrix(0, ntest, T)  # 
for(t in 1:T){
  all.test.id[, t] = sample(1:n, ntest)
}
#save(all.test.id, file="alltestID.RData")
#test.id = all.test.id[,1] 
#MSPE=rep(0,7)
MSPE =matrix(0,7,T)
names(MSPE) = c("Full", "R_min", "R_1se", "L_min", "L_1se", "L_Refit", "PCR")

Full Model

for (i in c(1:50)){
full.model = lm(Y ~ ., data = myData[-all.test.id[,i],])
Ytest.pred = predict(full.model, newdata = myData[all.test.id[,i],])
MSPE[1,i] = mean((myData$Y[all.test.id[,i]] - Ytest.pred)^2)
}

Ridge Model

mylasso.lambda.seq = exp(seq(-10, -2, length.out = 100))
for (i in c(1:50)){
cv.out = cv.glmnet(X[-all.test.id[,i], ], Y[-all.test.id[,i]], alpha = 0, 
                   lambda = mylasso.lambda.seq)
best.lam = cv.out$lambda.min
Ytest.pred = predict(cv.out, s = best.lam, newx = X[all.test.id[,i], ])
MSPE[2,i] = mean((Y[all.test.id[,i]] - Ytest.pred)^2)

best.lam = cv.out$lambda.1se
Ytest.pred = predict(cv.out, s = best.lam, newx = X[all.test.id[,i], ])
MSPE[3,i] = mean((Y[all.test.id[,i]] - Ytest.pred)^2)
}

Lasso

for (i in c(1:50)){
cv.out = cv.glmnet(X[-all.test.id[,i], ], Y[-all.test.id[,i]], alpha = 1)
best.lam = cv.out$lambda.min
Ytest.pred = predict(cv.out, s = best.lam, newx = X[all.test.id[,i], ])
MSPE[4,i] = mean((Y[all.test.id[,i]] - Ytest.pred)^2)

best.lam = cv.out$lambda.1se
Ytest.pred = predict(cv.out, s = best.lam, newx = X[all.test.id[,i], ])
MSPE[5,i] = mean((Y[all.test.id[,i]] - Ytest.pred)^2)

mylasso.coef = predict(cv.out, s = best.lam, type = "coefficients")
var.sel = row.names(mylasso.coef)[which(mylasso.coef != 0)[-1]]
mylasso.refit = lm(Y ~ ., myData[-all.test.id[,i], c("Y", var.sel)])
Ytest.pred = predict(mylasso.refit, newdata = myData[all.test.id[,i], ])
MSPE[6,i] = mean((Ytest.pred - Y[all.test.id[,i]])^2)
}

PCR

for (i in c(1:50)){
mypcr = pcr(Y ~ ., data= myData[-all.test.id[,i], ], validation="CV")
CVerr = RMSEP(mypcr)$val[1, , ]
adjCVerr = RMSEP(mypcr)$val[2, , ]
best.ncomp = which.min(CVerr) - 1 
if (best.ncomp==0) {
    Ytest.pred = mean(myData$Y[-all.test.id[,i]])
  } else {
    Ytest.pred = predict(mypcr, myData[all.test.id[,i],], ncomp=best.ncomp)
  }

MSPE[7,i] = mean((Ytest.pred - myData$Y[all.test.id[,i]])^2)
}

Overview for BostonData2

library(ggplot2)
par(pin = c(5,4))
boxplot(t(MSPE),main="Errors for Seven methods",col=terrain.colors(7))
#axis(2,seq(0,1,0.1),seq(0,1,0.1))
legend("topright", inset=0,   c("Full","Ridge_min","Ridge_lse","Lasso_min","Lasso_lse","L_refit","PCR"), fill=terrain.colors(7),bty='n')

So, We can find from the figure above that ridge regression using lambda.min is the best model.

Part II: Simulation Study of BostonData3

Firstly, load the data set, which has many noise data

myData3 = read.csv("https://liangfgithub.github.io/F21/BostonData3.csv")
myData3=myData3[,-1]

X = data.matrix(myData3[,-1])  
Y = myData3[,1]  

T = 50
n = length(Y)
ntest = round(n * 0.25)  # test set size
ntrain = n - ntest  # training set size
all.test.id = matrix(0, ntest, T)  # 
for(t in 1:T){
  all.test.id[, t] = sample(1:n, ntest)
}

MSPE3 =matrix(0,6,T)
names(MSPE3) = c( "R_min", "R_1se", "L_min", "L_1se", "L_Refit", "PCR")

Ridge Model

mylasso.lambda.seq = exp(seq(-10, -2, length.out = 100))
for (i in c(1:50)){
cv.out = cv.glmnet(X[-all.test.id[,i], ], Y[-all.test.id[,i]], alpha = 0, 
                   lambda = mylasso.lambda.seq)
best.lam = cv.out$lambda.min
Ytest.pred = predict(cv.out, s = best.lam, newx = X[all.test.id[,i], ])
MSPE3[1,i] = mean((Y[all.test.id[,i]] - Ytest.pred)^2)

best.lam = cv.out$lambda.1se
Ytest.pred = predict(cv.out, s = best.lam, newx = X[all.test.id[,i], ])
MSPE3[2,i] = mean((Y[all.test.id[,i]] - Ytest.pred)^2)
}

Lasso

for (i in c(1:50)){
cv.out = cv.glmnet(X[-all.test.id[,i], ], Y[-all.test.id[,i]], alpha = 1)
best.lam = cv.out$lambda.min
Ytest.pred = predict(cv.out, s = best.lam, newx = X[all.test.id[,i], ])
MSPE3[3,i] = mean((Y[all.test.id[,i]] - Ytest.pred)^2)

best.lam = cv.out$lambda.1se
Ytest.pred = predict(cv.out, s = best.lam, newx = X[all.test.id[,i], ])
MSPE3[4,i] = mean((Y[all.test.id[,i]] - Ytest.pred)^2)

mylasso.coef = predict(cv.out, s = best.lam, type = "coefficients")
var.sel = row.names(mylasso.coef)[which(mylasso.coef != 0)[-1]]
mylasso.refit = lm(Y ~ ., myData3[-all.test.id[,i], c("Y", var.sel)])
Ytest.pred = predict(mylasso.refit, newdata = myData3[all.test.id[,i], ])
MSPE3[5,i] = mean((Ytest.pred - Y[all.test.id[,i]])^2)
}

PCR

for (i in c(1:50)){
mypcr = pcr(Y ~ ., data= myData3[-all.test.id[,i], ], validation="CV")
CVerr = RMSEP(mypcr)$val[1, , ]
adjCVerr = RMSEP(mypcr)$val[2, , ]
best.ncomp = which.min(CVerr) - 1 
if (best.ncomp==0) {
    Ytest.pred = mean(myData3$Y[-all.test.id[,i]])
  } else {
    Ytest.pred = predict(mypcr, myData3[all.test.id[,i],], ncomp=best.ncomp)
  }

MSPE3[6,i] = mean((Ytest.pred - myData3$Y[all.test.id[,i]])^2)
}

Overview for BostonData3

library(ggplot2)
boxplot(t(MSPE3),main="Errors for Seven methods",col=terrain.colors(6))
#axis(2,seq(0,1,0.1),seq(0,1,0.1))
legend("topright", inset=0,   c("Ridge_min","Ridge_lse","Lasso_min","Lasso_lse","L_refit","PCR"), fill=terrain.colors(6),bty='n')

From the figure above, we can find that PCR is the worst model and ridge models are relatively worse than the rest models.

Actually, it’s easy to explain, because there are lots of noise factors, so if the data set form a new matrix, it will also contain noise factors.

Because lasso model is well-suited for models showing high levels of muticollinearity, and the coefficients will be zero for those noise factors.