이번에 살펴볼 예제는 Best Subset Selection에 이은 Forward Subset Selection이라는 방법입니다. 앞서 Best Subset Selection에서는 총 $p$개의 predictors에 대해서 $2^p$개의 models를 fitting해야하기 때문에 computationally expensive하다고 했었습니다.
그것 외에도 모든 가능한 predictors의 조합들을 고려하다 보니(searching space가 매우 넓겠죠?) 어쩌다가 낮은 error를 갖는 조합이 걸릴 수 있게 되고, 결국 predictors의 조합 자체를 overfitting하게 되는 결과를 초래할 수도 있게 됩니다.
이런 저런 문제들을 보완하기 위해 나온 방법이 바로 Forward Subset Selection입니다.(Backward Subset Selection도 거의 동일한 방법이니 여기선 스킵하겠습니다.) Greedy algorithm의 일종인 Forward Subset Selection은, null model(intercept만 있는 상황. 어떤 predictors도 선택되지 않은 상황)에서 시작해서 full model(모든 predictors를 선택한 모델)까지 각 stage에서 가장 작은 에러를 갖는 predictor들을 추가로 선택해 나가면서 searching하는 방법입니다.
처음 null model에서 하나의 predictor들만 사용해 $p$개의 모델을 만들어보고 하나의 predictor를 선택합니다. 다음엔 선택한 predictor를 제외한 $p-1$개의 predictors를 각각 기존에 선택했던 predictor와 조합하여 2개의 predictors로 모델을 $p-1$개 만들어보고 가장 낮은 MSE를 갖는 model을 선택합니다. 이런식으로 $p$개의 predictors를 선택할 때까지 계속 greedy하게 predictor들을 추가해 나갑니다. 결국 총 만드는 모델의 수는 \[p+(p-1)+(p-2)+…+2+1=\frac{p(p+1)}{2}\approx \frac{p^2}{2}\] 개가 되어 앞서 exponential 개수보다 훨씬 작은 polynomial 개수가 나오게 됩니다 :) 코드를 보면서 어떻게 fitting하는지 알아봅시다.
library(leaps) #regsubsets함수를 쓰기 위해 라이브러리를 포함합니다.
library(ISLR) #Hitters 데이터를 쓰기 위해 포함합니다.
# method를 forward로 함으로써 forward subset selection을 시행합니다.
regfit.fwd=regsubsets(Salary~.,data=Hitters,nvmax=19,method="forward")Summary에서는 각 stage에 어떤 variables가 선택되었는지 *를 통해 확인할 수 있습니다.
summary(regfit.fwd)## Subset selection object
## Call: regsubsets.formula(Salary ~ ., data = Hitters, nvmax = 19, method = "forward")
## 19 Variables (and intercept)
## Forced in Forced out
## AtBat FALSE FALSE
## Hits FALSE FALSE
## HmRun FALSE FALSE
## Runs FALSE FALSE
## RBI FALSE FALSE
## Walks FALSE FALSE
## Years FALSE FALSE
## CAtBat FALSE FALSE
## CHits FALSE FALSE
## CHmRun FALSE FALSE
## CRuns FALSE FALSE
## CRBI FALSE FALSE
## CWalks FALSE FALSE
## LeagueN FALSE FALSE
## DivisionW FALSE FALSE
## PutOuts FALSE FALSE
## Assists FALSE FALSE
## Errors FALSE FALSE
## NewLeagueN FALSE FALSE
## 1 subsets of each size up to 19
## Selection Algorithm: forward
## AtBat Hits HmRun Runs RBI Walks Years CAtBat CHits CHmRun CRuns
## 1 ( 1 ) " " " " " " " " " " " " " " " " " " " " " "
## 2 ( 1 ) " " "*" " " " " " " " " " " " " " " " " " "
## 3 ( 1 ) " " "*" " " " " " " " " " " " " " " " " " "
## 4 ( 1 ) " " "*" " " " " " " " " " " " " " " " " " "
## 5 ( 1 ) "*" "*" " " " " " " " " " " " " " " " " " "
## 6 ( 1 ) "*" "*" " " " " " " "*" " " " " " " " " " "
## 7 ( 1 ) "*" "*" " " " " " " "*" " " " " " " " " " "
## 8 ( 1 ) "*" "*" " " " " " " "*" " " " " " " " " "*"
## 9 ( 1 ) "*" "*" " " " " " " "*" " " "*" " " " " "*"
## 10 ( 1 ) "*" "*" " " " " " " "*" " " "*" " " " " "*"
## 11 ( 1 ) "*" "*" " " " " " " "*" " " "*" " " " " "*"
## 12 ( 1 ) "*" "*" " " "*" " " "*" " " "*" " " " " "*"
## 13 ( 1 ) "*" "*" " " "*" " " "*" " " "*" " " " " "*"
## 14 ( 1 ) "*" "*" "*" "*" " " "*" " " "*" " " " " "*"
## 15 ( 1 ) "*" "*" "*" "*" " " "*" " " "*" "*" " " "*"
## 16 ( 1 ) "*" "*" "*" "*" "*" "*" " " "*" "*" " " "*"
## 17 ( 1 ) "*" "*" "*" "*" "*" "*" " " "*" "*" " " "*"
## 18 ( 1 ) "*" "*" "*" "*" "*" "*" "*" "*" "*" " " "*"
## 19 ( 1 ) "*" "*" "*" "*" "*" "*" "*" "*" "*" "*" "*"
## CRBI CWalks LeagueN DivisionW PutOuts Assists Errors NewLeagueN
## 1 ( 1 ) "*" " " " " " " " " " " " " " "
## 2 ( 1 ) "*" " " " " " " " " " " " " " "
## 3 ( 1 ) "*" " " " " " " "*" " " " " " "
## 4 ( 1 ) "*" " " " " "*" "*" " " " " " "
## 5 ( 1 ) "*" " " " " "*" "*" " " " " " "
## 6 ( 1 ) "*" " " " " "*" "*" " " " " " "
## 7 ( 1 ) "*" "*" " " "*" "*" " " " " " "
## 8 ( 1 ) "*" "*" " " "*" "*" " " " " " "
## 9 ( 1 ) "*" "*" " " "*" "*" " " " " " "
## 10 ( 1 ) "*" "*" " " "*" "*" "*" " " " "
## 11 ( 1 ) "*" "*" "*" "*" "*" "*" " " " "
## 12 ( 1 ) "*" "*" "*" "*" "*" "*" " " " "
## 13 ( 1 ) "*" "*" "*" "*" "*" "*" "*" " "
## 14 ( 1 ) "*" "*" "*" "*" "*" "*" "*" " "
## 15 ( 1 ) "*" "*" "*" "*" "*" "*" "*" " "
## 16 ( 1 ) "*" "*" "*" "*" "*" "*" "*" " "
## 17 ( 1 ) "*" "*" "*" "*" "*" "*" "*" "*"
## 18 ( 1 ) "*" "*" "*" "*" "*" "*" "*" "*"
## 19 ( 1 ) "*" "*" "*" "*" "*" "*" "*" "*"Cp를 기준으로 어떤 모델이 가장 적합한지 plot해볼 수 있습니다.
plot(regfit.fwd,scale="Cp")