逐步選擇（向前、向后）

只需擬合p(p+1)/2個模型，運(yùn)算效率極大提高，但得到的模型可能非最優(yōu)模型。通常將向前向后逐步選擇結(jié)合使用，在加入新變量的同時剔除不能提升模型擬合效果的變量。

基于測試誤差選擇最優(yōu)模型：

1）??Cp值、AIC、BIC、調(diào)整R方，間接估計測試誤差

2）??驗(yàn)證集或交叉驗(yàn)證，直接估計測試誤差

library(ISLR)? ? #利用Hitters（棒球）數(shù)據(jù)集

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、調(diào)整R方 確定最優(yōu)模型

###最優(yōu)子集選擇#####
library(leaps)
regfit.full=regsubsets(Salary~.,Hitters)
summary(regfit.full)
#圖1 默認(rèn)給出最優(yōu)8變量模型的8種變量篩選結(jié)果

regfit.full=regsubsets(Salary~.,Hitters,nvmax=19) #最優(yōu)19變量模型的19種變量篩選結(jié)果
reg.summary=summary(regfit.full)
names(reg.summary)

#通過比較各變量模型的RSS、調(diào)整R方、Cp值和BIC值確定最優(yōu)變量模型
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隨著解釋變量數(shù)目增大而單調(diào)遞減，無法判斷
#圖b? ?調(diào)整R2——最優(yōu)模型為11變量模型
#圖c? ?Cp值——最優(yōu)模型為10變量模型
#圖d? ?BIC——最優(yōu)模型為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黑色方塊表示根據(jù)相應(yīng)的統(tǒng)計指標(biāo)選擇最優(yōu)模型所含的變量

###向前逐步選擇 向后逐步選擇####

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、 利用 驗(yàn)證集和交叉驗(yàn)證方法 確定最優(yōu)模型

###使用驗(yàn)證集和交叉驗(yàn)證 ####
set.seed(1)
train=sample(c(TRUE,FALSE),nrow(Hitters),rep=TRUE)
test=(!train)
regfit.best=regsubsets(Salary~.,data=Hitters[train,],nvmax=19) #訓(xùn)練數(shù)據(jù)模型擬合
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變量模型的系數(shù)
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模型回歸系數(shù)


#10折交叉驗(yàn)證
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次交叉驗(yàn)證誤差均值最小
reg.best=regsubsets(Salary~.,data=Hitters,nvmax = 19)
coef(reg.best,11)