子集選擇方法：最優子集選擇
#Hitters (棒球)數據集實踐最優於集選擇方法 library(ISLR) fix(Hitters) names(Hitters)
dim(Hitters)
sum(is.na(Hitters$Salary))
Hitters<-na.omit(Hitters) #刪除缺失值 dim(Hitters)
sum(is.na(Hitters)) #檢驗是否含有缺失值
library(leaps)
#R語言中的 regsubset ()函數(於leaps庫)通過建立一系列包含給定數目預測變量的最優模型，來實現最優預測變量子集的篩選，其中“最優”這一概念使用RSS來量化。
regfit.full=regsubsets(Salary~.,Hitters) summary(regfit.full)
結果分析：星號*表示列對應的變量包含於行對應的模型當中。例如，以上輸出結果表示最優的兩變量模型僅包含 Hits (前一年究成的安打次數)和 CRBI (職業生涯總跑壘次數)兩個變量。
regfit.full<-regsubsets(Salary~.,data=Hitters,nvmax=19) #默認果只給出截至最優八變量模型的八種變量篩選結果，nvmax 選項可以設置用戶所需的預測變量個數，這里設置為19。 reg.summary<-summary(regfit.full) names(reg.summary)
reg.summary$rsq
#par(mfrow=c(2,2)) #畫布分成4個區域，第一行有兩個區域，第二行也是二個區域
#畫出所有模型的 RSS 、調整 R2 ，及Cp 、BIC 的圖像可以輔助確定最終選擇哪一個模型
plot(reg.summary$rss,xlab="Number of Variables",ylab="RSS",type="l") #RSS
plot(reg.summary$adjr2,xlab="Number of Variables",ylab="Adjusted RSq",type="l") #R2
which.max(reg.summary$adjr2) #which. max ()函數可以用於識別一個向量中最大值所對應點的位置
points(11,reg.summary$adjr2[11], col="red",cex=2,pch=20)#points ()用於將點加在已有圖像上
plot(reg.summary$cp,xlab="Number of Variables",ylab="Cp",type='l') #Cp
which.min(reg.summary$cp) #which min() 函數標示出統計指標最小的模型。
points(10,reg.summary$cp[10],col="red",cex=2,pch=20) #points ()用於將點加在已有圖像上
which.min(reg.summary$bic) #which min() 函數標示出統計指標BIC最小的模型。
plot(reg.summary$bic,xlab="Number of Variables",ylab="BIC",type='l') #BIC
points(6,reg.summary$bic[6],col="red",cex=2,pch=20) #points ()用於將點加在已有圖像上
plot(regfit.full,scale="r2")
plot(regfit.full,scale="adjr2")
plot(regfit.full,scale="Cp")
plot(regfit.full,scale="bic")
結果分析：每個圖像第一行的黑色方塊表示根據相應統計指標選擇的最優模型所包含的變量。例如，多個模型的BIC都接近與-150。而BIC 指標最小的模型為包含變量AtBat (前一年的打數)、 Hits (前一年完成的安打次數)、 Walks (前一年的保送次數)、 CRBI (職業生涯單跑線次數)、 Division (年來翻聯賽級別)和 PutOuts (前一年的接殺次數)的六變量模型。
coef(regfit.full,6) #coef ()函數可以提取該模型的參數估計值。
子集選擇：向前逐步選擇
#向前逐步選擇和向后逐步選擇可以分別通過設定 regsubsets 函數中的參數 method =" forward" 和method =" backward" 來實現。
regfit.fwd=regsubsets(Salary~.,data=Hitters,nvmax=19,method="forward") summary(regfit.fwd)
結果分析：使用向前逐步選擇篩選變量時，最優的單變量模型只包含預測變量 CRBI ，最優的兩變量模型在此基礎上引人了預測變盤Hits。
子集選擇：向后逐步選擇
regfit.bwd=regsubsets(Salary~.,data=Hitters,nvmax=19,method="backward") summary(regfit.bwd)
coef(regfit.full,7) #最優子集選擇
coef(regfit.fwd,7) #向前逐步選擇
coef(regfit.bwd,7) #向后逐步選擇
結果分析：對於這個數據集，使用最優子集選擇方法和向前逐步選擇所得的單變量模型到六變盤模型的結果是完全一致的。但是，使用向前逐步選擇、向后逐步選擇和最優子集選擇所得的最優七變量模型是不同的。
模型選擇：驗證集方法
# 驗證集方法選擇模型，需要划分訓練集和測試集
set.seed(1)
Train<-sample(c(TRUE,FALSE), nrow(Hitters),rep=TRUE)
Test<-(!train)
regfit.best<-regsubsets(Salary~.,data=Hitters[train,],nvmax=19)
test.mat<-model.matrix(Salary~.,data=Hitters[test,])
val.errors<-rep(NA,19)
for(i in 1:19){
Coefi<-coef(regfit.best,id=i)
Pred<-test.mat[,names(coefi)]%*%coefi
val.errors[i]<-mean((Hitters$Salary[test]-pred)^2)
}
val.errors
which.min(val.errors) #通過驗證集方法，可以得出最優的模型含有7個預測變量
coef(regfit.best,7)
#編寫預測函數
predict.regsubsets=function(object,newdata,id,...){
form=as.formula(object$call[[2]])
mat=model.matrix(form,newdata)
coefi=coef(object,id=id)
xvars=names(coefi)
mat[,xvars]%*%coefi
}
regfit.best=regsubsets(Salary~.,data=Hitters,nvmax=19) #對整個數據集使用最優子集選擇，選出最優的七變量模型
coef(regfit.best,7)
模型選擇：交叉驗證法
k=10 set.seed(1) folds=sample(1:k,nrow(Hitters),replace=TRUE) ##定義一個向量將數據集中每個觀測歸為k= 10 折中的某一折，從1到k中隨機籌集有放回的抽取nrow(Hitters)個樣本 length(folds)
folds
cv.errors=matrix(NA,k,19, dimnames=list(NULL, paste(1:19))) #定義一個存儲計算結果的矩陣，k行19列。 cv.errors
#通過一個循環語句實現了交叉驗證。在第j 折中，數據集中對應 folds向量中等於j 的元素歸於測試集，其他元素歸於訓練集中
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) #使用apply ()函數求該矩陣的列平均（1表示對行，2表示1對列），從而得到一個向量，該列向量的第j個元素表示j變量模型的交叉驗證誤差。
mean.cv.errors
par(mfrow=c(1,1)) plot(mean.cv.errors,type='b')
結果分析：從上圖可知，交叉驗證選擇了十變最模型。
reg.best=regsubsets(Salary~.,data=Hitters, nvmax=19) #對整個數據集使用最優子集選擇，以獲得該十變最模型的參數估計結果。 coef(reg.best,10)
