#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的系数估计结果是稀疏的