利用glmnet函数实现岭回归和Lasso

#Hitters（棒球）数据集介绍

library(ISLR)

names(Hitters)

#[1] "AtBat"     "Hits"      "HmRun"   

# [4] "Runs"      "RBI"       "Walks"   

#[7] "Years"     "CAtBat"    "CHits"   

#[10] "CHmRun"    "CRuns"     "CRBI"     

#[13] "CWalks"    "League"    "Division"

#[16] "PutOuts"   "Assists"   "Errors"   

#[19] "Salary"    "NewLeague"

sum(is.na(Hitters$Salary))

Hitters=na.omit(Hitters)    #先剔除缺失值

####岭回归和Lasso####

x=model.matrix(Salary~.,Hitters)[,-1]  #删除了截距项

y=Hitters$Salary

library(glmnet)

grid=10^seq(10,-2,length=100)#10^-2至10^10

ridge.mod=glmnet(x,y,alpha=0,lambda=grid)

grid[50]   #11497.57

grid[60]   #705.4802

dim(coef(ridge.mod))

ridge.mod$lambda[50]  #lamuda=11498时的系数向量

ridge.mod$lambda[60]   #lamuda=705时的系数向量

#说明，比较上述两个系数向量，可知lamuda值小，岭回归的系数更大

#预测lamuda值为50的岭回归系数

predict(ridge.mod,s=50,type="coefficients")[1:20,] #将系数矩阵转化为向量

#划分训练集和测试集

set.seed(1)

train=sample(1:nrow(x),nrow(x)/2)

test=(-train)

y.test=y[test]

ridge.mod=glmnet(x[train,],y[train],alpha= 0,lambda = grid,thresh = 1e-12)

ridge.pred=predict(ridge.mod,s=4,newx=x[test,])

mean((ridge.pred-y.test)^2)  # lamuda=4时,MSE=101036.8

mean((mean(y[train])-y.test)^2)   #只含截距项模型的测试误差193253.1

ridge.pred=predict(ridge.mod,s=1e10,newx= x[test,])

mean((ridge.pred-y.test)^2)  # lamuda=10^10时,MSE=193253.1

#当lamuda很大时，模型的测试误差与只含截距项模型的测试误差相当

set.seed(1)

cv.out=cv.glmnet(x[train,],y[train],alpha=0)

par(mfrow=c(1,1))

plot(cv.out)

                              

bestlam=cv.out$lambda.min  #交叉验证确定最优的lamuda=211.7416

ridge.pred=predict(ridge.mod,s=bestlam,newx= x[test,])

mean((ridge.pred-y.test)^2)  # lamuda=212时,96015.51

out=glmnet(x,y,alpha=0)

predict(out,type="coefficients",s=bestlam)[1:20,]

#岭回归无法将系数变量压缩至0

###

#最小二乘即lamda=0的岭回归

ridge.pred=predict(ridge.mod,s=0,newx= x[test,],exact = TRUE)

mean((ridge.pred-y.test)^2)  #114783

lm.fit=lm(y~x,subset= train)

lasso.mod=glmnet(x[train,],y[train],alpha=1,lambda=grid)

plot(lasso.mod)

set.seed(1)

cv.out=cv.glmnet(x[train,],y[train],alpha=1)

plot(cv.out)

bestlam=cv.out$lambda.min  #交叉验证确定最优的lamuda=35.32063

lasso.pred=predict(lasso.mod,s=bestlam,newx= x[test,])

mean((lasso.pred-y.test)^2)   # lamuda=35时, MSE=100743.4

out=glmnet(x,y,alpha=1)

predict(out,type="coefficients",s=bestlam)[1:20,]

#Lasso的系数估计结果是稀疏的