在之前的內容泛化誤差和留出法，交叉驗證中，我們介紹了泛化誤差的概念以及用以估計泛化誤差的留出法、留一法和K折交叉驗證法，本節我們討論如何在R語言中進行實現。
本文以An Introduction to Statistical Learning with Applications in R中的Wage數據為例進行說明。我們用Age作為自變量對響應變量Wage進行建模，考慮多項式回歸模型：
當
首先考慮留出法，程序如下：
############hold-out####################library(ISLR)attach(Wage)set.seed(1)n=nrow(Wage)k=2subn=floor(n/k)train=sample(1: n, subn)cv.errorsho<-rep(0,5)
for (i in 1:5) { fit=lm(wage~poly(age,i), data=Wage, train) pred=predict(fit,newdata=Wage[-train,]) cv.errorsho[i]<-mean((pred-wage[-train])^2)}cv.errorshoplot(1:5,cv.errorsho,type="l")cv.errorsho
[1] 1712.319 1625.294 1622.104 1620.998 1629.320
可見，
接下來，我們考慮留一法，即Leave-one-out-cross-validatioin。對此，我們有兩種做法，一種是基於公式：
這裡
############LOOCV####################set.seed(1)cv.error<-rep(0,5)for (i in 1:5) { fit=lm(wage~poly(age,i),data=Wage) pred=predict(fit) X<-rep(1,n); for (l in 1: i){ X<-matrix(c(X,age^l),n) } H=X%*%ginv(t(X)%*%X)%*%t(X); h=diag(H); cv.error[i]<-mean((pred-wage)^2/(1-h)^2)}cv.errorplot(1:5,cv.error,type="l")cv.error
[1] 1676.235 1599.913 1594.727 1592.641 1592.153
儘管
另一種方法則直接使用R中自帶的cv.glm()函數。但需注意的是，cv.glm()並沒有使用上述簡潔的計算公式，所以在樣本量很大時，計算量會非常大。另外在調用cv.glm()時，需使用glm()來擬合模型。glm()命令中有family這個參數，若不設置則此時glm()的擬合結果和lm()的結果一致。程序如下：
##########cv.glm################set.seed(1)library(boot)cv.error=rep(0,5)for(i in 1:5){ glm.fit=glm(wage~poly(age,i)) cv.error[i]=cv.glm(Wage,glm.fit)$delta[1]}cv.errorplot(1:5,cv.error,type="l",xlab="Polynomial order", ylab="Mean Squared Error")cv.error
[1] 1676.235 1600.529 1595.960 1594.596 1594.879
可見，
最後，我們考慮K折交叉驗證。對此，我們也採用了兩種方法。第一個程序命令是按照K折交叉驗證的思路來編寫的，容易知道這種思路也適用於其他問題的模型評價和模型選擇，而第二個程序則採用cv.glm，此命令專門針對廣義線性回歸模型。具體程序如下：
##############10-fold CV###################set.seed(1)k=10n=nrow(Wage)subn=floor(n/k)location=sample(1: n, n)cv.errors1<-matrix(0,k,5)
for(j in 1:k){ for (i in 1:5) { test<-location[(subn*(j-1)+1): (subn*j)] train<-location[-(subn*(j-1)+1: subn*j)] fit=lm(wage~poly(age,i), data = Wage, train) pred = predict(fit,newdata=Wage[test,]) cv.errors1[j,i]<-mean((pred-wage[test])^2) }}cv.error1=apply(cv.errors1,2, mean)cv.error1plot(1:5,cv.error1,type="l",xlab="Polynomial order", ylab="Mean Squared Error")cv.error1
[1] 1673.669 1597.541 1592.690 1590.358 1590.418
###################################set.seed(1)cv.error2=rep(0,5)for(i in 1:5){ glm.fit=glm(wage~poly(age,i)) cv.error2[i]=cv.glm(Wage,glm.fit,K=10)$delta[1]}cv.error2plot(1:5,cv.error2,type="l",xlab="Polynomial order", ylab="Mean Squared Error")cv.error2
[1] 1676.826 1600.763 1598.399 1595.651 1594.977
以上討論均針對回歸問題，當響應變量是離散的時，cv.glm也是適用的，此時需增加family參數，並修改損失函數。具體如何調用，可通過?cv.glm或者cv.glm進行查看。
如果覺得本文不錯，請點讚關注！