[ISLR] Linear Model Selection and Regularization (2) (Validation set approach, k-fold CV)

이 글은 Introduction to Statistical Learning with R (ISLR)의 내용을 바탕으로 작성되었습니다.

 

저자들은 웹사이트를 통해 pdf 파일과 예제데이터, 예제코드를 제공하고 있습니다.

http://www-bcf.usc.edu/~gareth/ISL/

---

 Linear Model Selection and Regularization (1)에서는 크기가 다른 모델들 중에서 가장 좋은 모델을 고르기 위해 test MSE를 간접적으로 추정하는 방법들에 대해 알아봤습니다.

이번에는 Validation Set 기법과 k-fold Cross-Validation 기법을 이용해서 직접적으로 test MSE를 계산하는 방법에 대해 알아보겠습니다.

먼저 Validation Set 기법을 사용하기 위해서는 관측치들을 train set과 test set으로 분할해야 합니다.

set.seed(1)
train=sample(c(TRUE,FALSE), nrow(Hitters),rep=TRUE)
test=(!train)

Sample 함수를 사용해서 TRUE와 FALSE를 랜덤하게 Hitters의 관측치 수만큼 뿌려서 train이라는 벡터를 만듭니다.

이 벡터는 인덱스의 역할을 하게 됩니다.

> sum(train)
[1] 134

> length(train)
[1] 263

확인해보니 전체 263개의 데이터 중에 134개가 TRUE입니다. train을 인덱스로 사용하게 되면 134개의 관측치가 train set으로 할당되는 것입니다.

test set에는 train set에 포함되지 않은 나머지 관측치들이 들어가게 됩니다.

이제, regsubsets() 함수를 train set에 적용해서 best subset selection을 해줍니다.

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)
coef(regfit.best,10)

> val.errors
 [1] 220968.0 169157.1 178518.2 163426.1 168418.1 171270.6 162377.1 157909.3 154055.7
[10] 148162.1 151156.4 151742.5 152214.5 157358.7 158541.4 158743.3 159972.7 159859.8
[19] 160105.6

> which.min(val.errors)
[1] 10

> coef(regfit.best,10)
(Intercept)       AtBat        Hits       Walks      CAtBat       CHits      CHmRun 
-80.2751499  -1.4683816   7.1625314   3.6430345  -0.1855698   1.1053238   1.3844863 
     CWalks     LeagueN   DivisionW     PutOuts 
 -0.7483170  84.5576103 -53.0289658   0.2381662

크기가 다른 모델들의 test MSE를 계산해서 test MSE가 가장 작은 모델을 선택합니다. 여기서는 predictor 10개가 포함된 모델의 test MSE가 가장 작다는 것을 알 수 있습니다.

가장 좋은 모델이 무엇인지 알았다면 모든 데이터를 사용해서 다시 모델을 만들어줘야 합니다. 

regfit.best=regsubsets(Salary~.,data=Hitters,nvmax=19)
coef(regfit.best,10)

이제 k-fold Cross-Validation 기법을 사용해보겠습니다.

10-fold Cross-Validation을 사용하기 위해 10개의 fold를 만들어 줍니다.

k=10
set.seed(1)
folds=sample(1:k,nrow(Hitters),replace=TRUE)

각각의 fold는 다음과 같은 개수로 나뉘어 졌습니다.

> sum(folds==1)
[1] 13

> sum(folds==2)
[1] 25

> sum(folds==3)
[1] 31

> sum(folds==4)
[1] 32

> sum(folds==5)
[1] 33

> sum(folds==6)
[1] 27

> sum(folds==7)
[1] 26

> sum(folds==8)
[1] 30

> sum(folds==9)
[1] 22

> sum(folds==10)
[1] 24

cross-validation error를 계산해줍니다.

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

> mean.cv.errors
       1        2        3        4        5        6        7        8        9 
160093.5 140196.8 153117.0 151159.3 146841.3 138302.6 144346.2 130207.7 129459.6 

      10       11       12       13       14       15       16       17       18 
125334.7 125153.8 128273.5 133461.0 133974.6 131825.7 131882.8 132750.9 133096.2 

      19 
132804.7

plot을 그려봅니다.

plot(mean.cv.errors,type='b')

 

 

predictor 개수가 11개인 모델이 가장 좋다는 것을 알 수 있습니다.

전체 데이터를 사용해서 모델을 다시 만들어줍니다.

 

reg.best=regsubsets(Salary~.,data=Hitters, nvmax=19)
coef(reg.best,11)

 

> coef(reg.best,11)

 (Intercept)        AtBat         Hits        Walks       CAtBat        CRuns 
 135.7512195   -2.1277482    6.9236994    5.6202755   -0.1389914    1.4553310 

        CRBI       CWalks      LeagueN    DivisionW      PutOuts      Assists 
   0.7852528   -0.8228559   43.1116152 -111.1460252    0.2894087    0.2688277