만약 우리가 가진 데이터가 너무 많은 predictor들을 가지고 있으면 어떻게 할까요? 이런 경우 overfitting에 시달릴 가능성도 커지고, 분석 또한 어려워지겠죠. 그렇기 때문에 다양한 selection method가 등장하게 되었습니다. 그 중 오늘은 가장 직관적인 Best Subset Selection방법에 대해 알아보겠습니다.
우리가 가진 데이터에 총 $p$개의 predictor가 존재한다고 가정해 봅시다. 그렇다면 아무 predictor도 선택하지 않는 null model부터 시작해서, 모든 $p$개의 predictor를 다 사용하는 full model까지 총 $2^p$가지의 model을 fitting할 수 있습니다. 결국 brute-force하게 모든 가능한 model들을 만들어 보고 성능을 비교하여 최고의 model을 사용하겠다는 방법입니다. Complexity가 exponentially 증가하기 때문에 적당한 수($p\le 20$)의 predictors를 가지고 있을 때만 사용이 가능합니다.
일단 우리가 이번에 예제로 사용할 dataset은 Hitters입니다.
library(ISLR)
summary(Hitters)## AtBat Hits HmRun Runs
## Min. : 16.0 Min. : 1 Min. : 0.00 Min. : 0.00
## 1st Qu.:255.2 1st Qu.: 64 1st Qu.: 4.00 1st Qu.: 30.25
## Median :379.5 Median : 96 Median : 8.00 Median : 48.00
## Mean :380.9 Mean :101 Mean :10.77 Mean : 50.91
## 3rd Qu.:512.0 3rd Qu.:137 3rd Qu.:16.00 3rd Qu.: 69.00
## Max. :687.0 Max. :238 Max. :40.00 Max. :130.00
##
## RBI Walks Years CAtBat
## Min. : 0.00 Min. : 0.00 Min. : 1.000 Min. : 19.0
## 1st Qu.: 28.00 1st Qu.: 22.00 1st Qu.: 4.000 1st Qu.: 816.8
## Median : 44.00 Median : 35.00 Median : 6.000 Median : 1928.0
## Mean : 48.03 Mean : 38.74 Mean : 7.444 Mean : 2648.7
## 3rd Qu.: 64.75 3rd Qu.: 53.00 3rd Qu.:11.000 3rd Qu.: 3924.2
## Max. :121.00 Max. :105.00 Max. :24.000 Max. :14053.0
##
## CHits CHmRun CRuns CRBI
## Min. : 4.0 Min. : 0.00 Min. : 1.0 Min. : 0.00
## 1st Qu.: 209.0 1st Qu.: 14.00 1st Qu.: 100.2 1st Qu.: 88.75
## Median : 508.0 Median : 37.50 Median : 247.0 Median : 220.50
## Mean : 717.6 Mean : 69.49 Mean : 358.8 Mean : 330.12
## 3rd Qu.:1059.2 3rd Qu.: 90.00 3rd Qu.: 526.2 3rd Qu.: 426.25
## Max. :4256.0 Max. :548.00 Max. :2165.0 Max. :1659.00
##
## CWalks League Division PutOuts Assists
## Min. : 0.00 A:175 E:157 Min. : 0.0 Min. : 0.0
## 1st Qu.: 67.25 N:147 W:165 1st Qu.: 109.2 1st Qu.: 7.0
## Median : 170.50 Median : 212.0 Median : 39.5
## Mean : 260.24 Mean : 288.9 Mean :106.9
## 3rd Qu.: 339.25 3rd Qu.: 325.0 3rd Qu.:166.0
## Max. :1566.00 Max. :1378.0 Max. :492.0
##
## Errors Salary NewLeague
## Min. : 0.00 Min. : 67.5 A:176
## 1st Qu.: 3.00 1st Qu.: 190.0 N:146
## Median : 6.00 Median : 425.0
## Mean : 8.04 Mean : 535.9
## 3rd Qu.:11.00 3rd Qu.: 750.0
## Max. :32.00 Max. :2460.0
## NA's :59요약 결과를 보시면 아시겠지만, missing value들이 존재하지요? 따라서 더 진행하기 전에 이것들을 일단 제거해줍시다.(왜냐하면 linear model들은 NA를 다룰 수 없기 때문입니다.) 이것은 na.omit함수로 쉽게 할 수 있습니다.
Hitters=na.omit(Hitters)
with(Hitters,sum(is.na(Salary))) #NA인 것이 있다면 0보다 큰 수가 나올 것입니다.## [1] 0자 이제 본격적으로 best subset selection을 위해 leaps 라이브러리를 불러옵시다. 그리고 regsubsets함수에 우리가 친숙한 Y~X...형태로 공식을 넣어주면 알아서 계산을 해 줍니다 :)
library(leaps)
regfit.full=regsubsets(Salary~.,data=Hitters)
summary(regfit.full)## Subset selection object
## Call: regsubsets.formula(Salary ~ ., data = Hitters)
## 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 8
## Selection Algorithm: exhaustive
## 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 ) "*" "*" " " " " " " "*" " " " " " " "*" "*"
## CRBI CWalks LeagueN DivisionW PutOuts Assists Errors NewLeagueN
## 1 ( 1 ) "*" " " " " " " " " " " " " " "
## 2 ( 1 ) "*" " " " " " " " " " " " " " "
## 3 ( 1 ) "*" " " " " " " "*" " " " " " "
## 4 ( 1 ) "*" " " " " "*" "*" " " " " " "
## 5 ( 1 ) "*" " " " " "*" "*" " " " " " "
## 6 ( 1 ) "*" " " " " "*" "*" " " " " " "
## 7 ( 1 ) " " " " " " "*" "*" " " " " " "
## 8 ( 1 ) " " "*" " " "*" "*" " " " " " "디폴트값으로 총 8개까지의 predictor의 best model만 구하게 됩니다. 동일한 개수의 predictor로 구성된 모델들 중에서는 training RSS에 의해 best model이 결정되게 됩니다. 요약된 것 중에 *로 표시된 것들이 선택된 predictor들입니다.
우리가 가진 데이터셋은 총 19개의 predictor들이 있으니, 이 개수에 대해서 best model들을 구해봅시다. nvmax값을 변경함으로써 가능합니다 :)
regfit.full=regsubsets(Salary~.,data=Hitters, nvmax=19)
reg.summary=summary(regfit.full)
names(reg.summary)## [1] "which" "rsq" "rss" "adjr2" "cp" "bic" "outmat" "obj"Summary의 다양한 항목들이 있음을 확인할 수 있죠? 앞서 동일한 수의 predictor를 가지는 model에 대해서는 training RSS로 predictor들을 선택한다고 했는데요, 그렇다면 몇 개의 predictor를 사용해야 할까요? 총 1~19개까지 19개의 model들이 있는데요, 이것들도 마찬가지로 training RSS를 사용해서 비교가 가능할까요?
정답은 결코! 아닙니다!! 그 이유는 당연히 full model(모든 predictor들을 포함한 model)이 가장 낮은 training RSS를 기록할 것이기 때문입니다! (predictor를 추가하는 것은 최소한 동일한 training RSS를 얻습니다. 하지만 이는 overfitting으로 연결되서 test RSS는 악화될 수 있죠.) 따라서 여기서는 다양한 방법이 존재하지만, 가장 많이 사용되는 방법이 Mallows’s_Cp, Bayesian_information_criterion 등 변형된 지표를 사용해야 합니다. 이번에는 Cp값을 이용해서 몇 개의 predictor를 사용한 model이 가장 좋을지 판단해 보죠.
#일단 전체적인 cp값을 plot하고,
plot(reg.summary$cp,xlab="Number of Variables",ylab="Cp")
#min값을 찾아서 동그라미 쳐 줍니다.
which.min(reg.summary$cp)## [1] 10points(10,reg.summary$cp[10],pch=20,col="red")
10개의 predictor를 사용했을 때가 가장 좋다는 결론이 Cp값을 사용했을 때 나왔네요 :) 그럼 마지막으로 조금 혁신적인 plot(?)을 끝으로 마무리하겠습니다. 각 Cp값마다 어떤 predictor들을 사용했는지 블록으로 표시해주는 좋은 plot을 regsubsets model이 제공합니다.
plot(regfit.full,scale="Cp")
또한 predictor의 수에 따른 선택된 predictor들의 계수들도 coef함수를 통해 쉽게 얻을 수 있습니다.
coef(regfit.full,10)## (Intercept) AtBat Hits Walks CAtBat
## 162.5354420 -2.1686501 6.9180175 5.7732246 -0.1300798
## CRuns CRBI CWalks DivisionW PutOuts
## 1.4082490 0.7743122 -0.8308264 -112.3800575 0.2973726
## Assists
## 0.2831680