逐步选择（向前、向后）

只需拟合p(p+1)/2个模型，运算效率极大提高，但得到的模型可能非最优模型。通常将向前向后逐步选择结合使用，在加入新变量的同时剔除不能提升模型拟合效果的变量。

基于测试误差选择最优模型：

1）  Cp值、AIC、BIC、调整R方，间接估计测试误差

2）  验证集或交叉验证，直接估计测试误差

library(ISLR)    #利用Hitters（棒球）数据集

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"

dim(Hitters)
sum(is.na(Hitters$Salary))
Hitters=na.omit(Hitters)    #先剔除缺失值
dim(Hitters)

1、利用 Cp、AIC、BIC、调整R方 确定最优模型

###最优子集选择#####
library(leaps)
regfit.full=regsubsets(Salary~.,Hitters)
summary(regfit.full)
#图1 默认给出最优8变量模型的8种变量筛选结果

regfit.full=regsubsets(Salary~.,Hitters,nvmax=19) #最优19变量模型的19种变量筛选结果
reg.summary=summary(regfit.full)
names(reg.summary)

#通过比较各变量模型的RSS、调整R方、Cp值和BIC值确定最优变量模型
par(mfrow=c(2,2))
plot(reg.summary$rss,xlab = 'Number of Variables',ylab ='RSS',type='l')
plot(reg.summary$adjr2,xlab = 'Number of Variables',ylab = 'AdjustedRSq',type='l')
which.max(reg.summary$adjr2)
points(11,reg.summary$adjr2[11],col='red',cex=2,pch=20)

plot(reg.summary$cp,xlab='number of variables',ylab='Cp',type='l')
which.min(reg.summary$cp)
points(10,reg.summary$cp[10],col='red',cex=2,pch=20)

plot(reg.summary$bic,xlab='number of variables',ylab='Bic',type='l')
which.min(reg.summary$bic)
points(6,reg.summary$bic[6],col='red',cex=2,pch=20)
#图a   RSS随着解释变量数目增大而单调递减，无法判断
#图b   调整R2——最优模型为11变量模型
#图c   Cp值——最优模型为10变量模型
#图d   BIC——最优模型为6变量模型

par(mfrow=c(2,2))
plot(regfit.full,scale = "r2")
plot(regfit.full,scale = "adjr2")
plot(regfit.full,scale = "Cp")
plot(regfit.full,scale = "bic")


#图2黑色方块表示根据相应的统计指标选择最优模型所含的变量

###向前逐步选择 向后逐步选择####

regfit.fwd=regsubsets(Salary~.,data=Hitters,nvmax=19,method='forward')
summary(regfit.fwd)

regfit.bwd=regsubsets(Salary~.,data=Hitters,nvmax=19,method='backward')
summary(regfit.bwd)

2、 利用 验证集和交叉验证方法 确定最优模型

###使用验证集和交叉验证 ####
set.seed(1)
train=sample(c(TRUE,FALSE),nrow(Hitters),rep=TRUE)
test=(!train)
regfit.best=regsubsets(Salary~.,data=Hitters[train,],nvmax=19) #训练数据模型拟合
summary(regfit.best)

test.mat=model.matrix(Salary~.,data=Hitters[test,])
val.errors=rep(NA,19)
for(i in 1:19){
  coefi=coef(regfit.best,id=i)  # i变量模型的系数
pred=test.mat[,names(coefi)]%*%coefi
val.errors=mean((Hitters$Salary[test]-pred)^2)
}
val.errors
which.min(val.errors)   #图3 测试误差最小的为10模型
coef(regfit.best,10)   #10模型回归系数


#10折交叉验证
k=10
set.seed(1)
folds=sample(1:k,nrow(Hitters),replace = TRUE)
cv.errors=matrix(NA,k,19,dimnames = list(NULL,paste(1:19)))

for(j in 1:k){
best.fit=regsubsets(Salary~.,data=Hitters[folds!=j,],nvmax=19)
  for(i in 1:19){
    pred=predict(best.fit,Hitters[folds==j,],id=i)
   cv.errors[j,i]=mean((Hitters$Salary[folds==j]-pred)^2)
  }
}
mean.cv.errors=apply(cv.errors,2,mean)
mean.cv.errors
which.min(mean.cv.errors)  #图4 11变量模型的10次交叉验证误差均值最小
reg.best=regsubsets(Salary~.,data=Hitters,nvmax = 19)
coef(reg.best,11)