기존 변수 선택 방법의 한계
나는 이 글을 통해 다중회귀모델의 변수 선택의 최종적인 Guideline을 제시하고자 한다. 고백하자면 이 전글 회귀모델의 변수 선택 방법을 통해 표준 Guideline을 제시하고자 했었다. 또한 추가 글 UCLA 대학원 입학을 위한 Spec.을 통해 최적 모델이란 전체 데이터가 아닌 검정(test) 데이터에서의 오차가 최소되어야 함을 설명했다. 그러나 이런 방법들은 여전히 문제점을 가지고 있다.
따라서 이 글을 통해 1. 앞서 소개한 방법의 문제점을 보완하고 2.변수 선택시 우리만의 경쟁 우위를 가지는 Guideline을 제시하고자 한다. 특히, 설명 변수가 많을 경우에 이런 경쟁 우위는 더욱 두드러지게 될 것이다. 또한 앞서 사용된 함수, ‘step’ 대신 3.하나의 함수로 광범위한 분석이 가능한 범용성 강한 함수, ‘regsubsets’을 적용하고자 한다.
비록 앞서 소개한 방법이 “일반적인 변수 선택법”으로 널리 사용됨에도 불구하고 다음의 문제점을 가진다.
전진, 후진, 단계를 통한 변수 선택은 전체 부분집합(subset) 중 일부만 비교되므로 실제 최적의 설명변수를 간과할 수 있다.
일반적으로 변수 갯수가 적을(5개 이하) 경우에는 크게 문제되지 않으나 많을(6개 이상) 경우에는 이 방법이 전체를 검토하는 부분집합 결과와 달라질 수 있다.(이 글 아래에서 증명한다)임의 표본(sampling)한 훈련셋 중 테스트셋을 가장 잘 설명한 경우를 최적 모델로 선택할 경우, 이 모델이 전체 데이터를 잘 설명한다는 보장이 없다. “UCLA대학원 입학”글에서는 검증셋 기법(validation Set)을 사용하였다. 이는 테스트셋의 데이터 수를 실제 반응변수(outcome) 빈도를 기준으로 정하는 점에서 훈련셋과 테스트셋을 정확히 두 부분으로 나누는 일반적인 검증셋 기법(validation Set)과 차이가 있으나 이 역시 일종의 검증셋 기법이다. 검증셋 기법은 어느 관측치가 훈련셋과 검증셋에 포함되느냐에 따라 변동이 상당히 크기 때문에 이 전글처럼 가장 우수한 값을 가지는 표본을 선택하여 이를 최적 모델이라 결론 짓는 것은 답을 정해 놓고 우연성이 그 답을 내놓길 기다리는 것과 같은 오류이다.
Cp, AIC, BIC, Adjusted R2통한 변수 선택은 이상값 제거시 결과가 달라 질 수 있으므로 이상값을 먼저 제거할 지, 변수를 먼저 선정할지에 대한 “닭이 먼저냐 달걀이 먼저냐”의 문제에 직면하게 된다.
진화된 변수 선택 방법(Mean Squared Error, MSE 기준)
이전 방법의 문제점을 통해 진화된 변수 선택법은 1.일부가 아닌 전체 부분 집합을 고려해야하며 2.훈련 및 테스트셋의 임의성이 일반화될 수 있어야 하며 3.이상값 제거 후에도 변수 선택의 변동성이 작아야 함을 알 수 있다. “가급적”이란 표현은 이 글의 대한 배려이다.
평균제곱오차 MSE는 각 데이터의 실제값에서 추정값을 뺀 값 제곱의 평균이다.
MSE가 최소라는 것은 실제 추정되는 함수가 실데이터를 가장 잘 설명한다고 할 수 있어 좋은 모델 선정의 하나의 기준이 된다. 그러나 만약 이 이유만으로 최소 MSE를 변수 선택의 기준으로 삼는다면 변수가 많을수록 MSE는 작아지므로 최대한 많은 변수를 선택하는 오류에 직면하게 된다. Criterion 접근법(AIC, BIC, Cp, Adjusted R2)은 이를 보완하고자 변수 갯수와 MSE간의 균형점을 제시하는 방법임을 이전 글에서 소개했다. 만약 Criterion을 전진, 후진, 단계가 아닌 Cpplot등의 전체 부분집합을 통해 확인한다면 앞서 제기된 문제점 1은 해결될 것이다. 그러나 이 방법은 검증을 하지 못하므로 전체셋이 실제 모델을 잘 설명한다는 대전제에서만 유효하다 할 수 있다.
전체 부분집합을 고려하면서 테스트셋의 임의성 또한 일반화(문제점 2)하려면 어떻게 해야할까? 교차 검증(cross validation) 방법 중 5-fold 또는 10-fold CV법이 그 해결책이다. 테크니컬한 요점은 전체 데이터를 5개 또는 10개 그룹으로 나눠 각 그룹이 한번씩 검정셋(test)이 되어 총 5번 또는 10번의 MSE의 평균을 구한다. (참고로 LOOCV(Leave-One-Out)방법은 10개 그룹이 아닌 전체 데이터 수 만큼 그룹화하여 하나하나 검정하는 전수 조사법이다. LOOCV는 거의 동일한 관측치들로 구성된 훈련셋을 사용하기 때문에 적합된 모델의 결과들이 서로 높은 상관성이 있다. 따라서 LOOCV는 k-folds에 비해 더 큰 분산을 가진다. 일반적으로 5-folds나 10folds 가 편향-분상 절충하는 것으로 알려져있다.
교차 검증의 가장 큰 장점은 검증셋 기법과 같이 훈련과 테스트셋으로 구분하나 검증셋 기법이 너무 많은 데이터를 테스트셋에 할당(예: 2분의 1)해 정작 충분한 훈련 데이터 수를 확보하지 못하는데 비해 검증 후에도 전체 데이터를 통해 변수 선택과 계수를 구한다는점이다.
이것은 1차적으로 검증셋 기법이 표본이 달라지면 변동이 큰데 반해 CV방법은 여러번 표본 수행에 대한 변수 갯수에 대한 검정 MSE 변동이 작기 때문에 가능하다.(위 그림 참조) 교차 검증은 변수 갯수에 대한 훈련셋과 테스트셋을 통해 변수의 갯수만을 선택한 후 최종적인 변수 선택은 전체셋을 통해 구한다.
결론을 내기 전에 다음의 의문 3가지를 살펴보자.
- 왜 전수검증법인 LOOCV는 사용하지 않는가?
- 선형회귀의 가정은 어느 정도의 신뢰성이 있을까?
- MSE만 고려한 교차 검증에 Criterion법과 같이 변수 갯수를 추가하면 어떨까?
위 그림은 우리가 찾고자 하는 실제 함수가 각각 유연한 스플라인(좌), 선형(중앙), 강한 스플라인(우)일 때의 추정 함수의 Flexibility(2=선형 회귀, 5=약한 평활스플라인, 20 강한 평활스플라이)에 따른 검정 MSE이다. 이 그래프는 앞서 제기한 의문 1번, 2번에 답을 준다.(진짜 검정MSE(파란색), LOOCV(검정 파선), 10-fold(오렌지))
- 10-fold CV는 LOOCV와 오차가 크게 다르지 않는 좋은 방법이다.
- 선형으로 추정하는 것이(Flexibility=2) 어느 정도의 신뢰성이 있을까?라는 의문에 실제 함수가 강한 스플라인만 아니라면 선형 회귀는 적절한 신뢰성을 가진다(강한 스플라인일 경우 과소추정)
따라서 10-fold 교차 검증은 (실제가 강한 스플라인만 아니라면) 실제 검정 MSE를 잘 설명할 뿐 아니라 표본의 임의성 또한 일반화한다 할 수 있다.
상기 그림은 변수 갯수에 따른 BIC의 제곱, 검증 기법, 교차 검증의 관계를 보여준다. 각 모델이 제시하는 최적 변수 갯수는 모델 별로 차이가 있으나 상기 그림은 교차 검증이 BIC와 유사한 기준임을 보여 준다. 즉, 변수 갯수를 포함하지 않아도 변수 갯수를 잘 설명함을 알 수 있다.(변수 갯수가 증가해도 검증 기법을 통한 MSE는 더이상 작아지지 않는다.)
결론
- 10 fold CV는 변수 선택의 최우선 방법으로 사용한다. 특히, 변수가 많을 경우 이 방법은 막대한 경쟁우위를 가져다 줄 것이다.
- BIC, AIC, CP를 통한 전체부분집합 변수 선택 방법은 여전히 훌륭한 대안이다.
- 디테일하게는 10 fold CV를 30번 수행하여 표준오차를 구해 1 표준오차 범위에서 MSE가 가장 작은 모델 갯수를 확인하면 좋다. (Option)
- 값(예: 작은 MSE나 AIC등)보다도 해석이 잘되는 모델이 더 중요하며 상황에 따라 선택 기준은 달라질 수 있다. (예: MSE최소 대신 적은 변수 갯수 등)
실전 예제
야구 선수 기록(설명변수)과 연봉(반응변수)이 포함된 데이터를 ISLR 패키지에서 불러온다. 그리고 regsubsets 함수 사용을 위해 leaps 패키지를 부른다.
library(ISLR)
library(leaps)
Hitters=na.omit(Hitters)
regfit.full=regsubsets(Salary~.,Hitters,nvmax=19)1개 변수에서 nvmax 변수 갯수까지 최소 RSS가 되는 설명 변수를 선택하여 보여주며 최대 19까지 할 수 있다. 참고로 MSE는 RSS의 평균이다.데이터가 20개이므로 19를 한다. summary함수가 어떤 변수가 선택되었는지 보여 준다.
summary(regfit.full)## 1 subsets of each size up to 19
## 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 ) "*" "*" " " " " " " "*" " " " " " " "*" "*"
## 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 ) "*" "*" "*" "*" "*" "*" "*" "*"reg.summary=summary(regfit.full)결과 값은 변수가 2개일 경우 Hits와 CRBI의 설명변수 포함 모델이 최소의 RSS를 가짐을 말해준다. summary는 또한 아래와 같이 Adjusted R2, BIC, Cp를 포함한다.
names(reg.summary)## [1] "which" "rsq" "rss" "adjr2" "cp" "bic" "outmat" "obj"name을 통해 확인한 변수명을 통해 adjusted R2, Cp, BIC를 plot해보자. plot 전에 각 값이 최적이 되는 위치를 확인해 보자.
which.max(reg.summary$adjr2)## [1] 11which.min(reg.summary$cp)## [1] 10which.min(reg.summary$bic)## [1] 6각각 11, 10, 6개 변수 모델이 각각 최적임을 알 수 있다. plot과 최적점을 표시해 보자.
par(mfrow=c(2,2))
plot(reg.summary$adjr2,type="l",xlab="No. of variables", ylab="Adjusted R2")
points(11,reg.summary$adjr2[11],col="red",cex=2,pch=20)
plot(reg.summary$cp,ylab="Cp",type="l")
points(10,reg.summary$cp[10],col="red",cex=2,pch=20)
plot(reg.summary$bic, ylab="BIC",type="l")
points(6,reg.summary$bic[6],col="red",cex=2,pch=20)
교차검증을 통해 변수 선택을 해보자. 우선 10 folds CV를 통해 최적 모델의 변수 갯수는 몇 개인지를 구한다. 샘플링은 랜덤 시드를 설정하고 복원추출한다.
k<-10
set.seed(1)
folds<-sample(1:k,nrow(Hitters),replace=TRUE)
cv.errors<-matrix(NA,k,19,dimnames=list(NULL,paste(1:19)))#10x19 NA 생성
for(j in 1:k){
best.fit=regsubsets(Salary~.,data=Hitters[folds!=j,],nvmax=19)
for(i in 1:19){
test.mat.fold<-model.matrix(Salary~.,data=Hitters[folds==j,])
coefi.fold<-coef(best.fit,id=i)
pred.fold<-test.mat.fold[,names(coefi.fold)]%*%coefi.fold
cv.errors[j,i]<-mean((Hitters$Salary[folds==j]-pred.fold)^2)
}
}
mean.cv.errors<-apply(cv.errors,2,mean) #각 변수 갯수별 평균값
mean.cv.errors## 1 2 3 4 5 6 7 8
## 160093.5 140196.8 153117.0 151159.3 146841.3 138302.6 144346.2 130207.7
## 9 10 11 12 13 14 15 16
## 129459.6 125334.7 125153.8 128273.5 133461.0 133974.6 131825.7 131882.8
## 17 18 19
## 132750.9 133096.2 132804.7par(mfrow=c(1,1))
plot(mean.cv.errors,type='b')
변수 갯수가 11개 일때 검정 MSE가 최소임을 알 수 있다.(10와는 큰 차이가 없다.) 이는 앞서 검증없는 BIC 결과와(변수 갯수 6개) 크게 다름을 알 수 있다. 이제 전체 데이터를 통해 변수를 선택하고 각 변수의 회귀계수를 구하면 완료가 된다.
coef(regfit.full,11)## (Intercept) AtBat Hits Walks CAtBat
## 135.7512195 -2.1277482 6.9236994 5.6202755 -0.1389914
## CRuns CRBI CWalks LeagueN DivisionW
## 1.4553310 0.7852528 -0.8228559 43.1116152 -111.1460252
## PutOuts Assists
## 0.2894087 0.2688277추가 확인 필요 사항.
- 이상값 제거시 최적 모델의 변수 갯수 변동이 어느 정도인지 확인이 필요하다.
- 알수 없는 실제 값이 강한 스플라인일 경우엔 선형회귀 검정은 과소 추정하므로 (ISLR 5.6참조) 선형이 아닌 다른 모델 적용이 필요하다.(ISLR 그림 2.9, 2.10, 2.11 참조)
- 상호작용항이나 비선형항에 대한 추가가 필요하다.
부록 - 유용한 팁들
1. 전진, 후진, 선택법
regfit.fwd=regsubsets(Salary~.,Hitters,nvmax=19, method="forward")
regfit.back=regsubsets(Salary~.,Hitters,nvmax=19, method="backward")6개 또는 7개가 포함된 변수 모델을 부분집합, 전진, 후진간에 비교해보면 다음과 같다.
coef(regfit.full,6)## (Intercept) AtBat Hits Walks CRBI
## 91.5117981 -1.8685892 7.6043976 3.6976468 0.6430169
## DivisionW PutOuts
## -122.9515338 0.2643076coef(regfit.fwd,6)## (Intercept) AtBat Hits Walks CRBI
## 91.5117981 -1.8685892 7.6043976 3.6976468 0.6430169
## DivisionW PutOuts
## -122.9515338 0.2643076coef(regfit.back,6)## (Intercept) AtBat Hits Walks CRuns
## 78.2664070 -1.8158931 7.3597644 3.5123248 0.6187876
## DivisionW PutOuts
## -113.7958600 0.2995788coef(regfit.full,7)## (Intercept) Hits Walks CAtBat CHits
## 79.4509472 1.2833513 3.2274264 -0.3752350 1.4957073
## CHmRun DivisionW PutOuts
## 1.4420538 -129.9866432 0.2366813coef(regfit.fwd,7)## (Intercept) AtBat Hits Walks CRBI
## 109.7873062 -1.9588851 7.4498772 4.9131401 0.8537622
## CWalks DivisionW PutOuts
## -0.3053070 -127.1223928 0.2533404coef(regfit.back,7)## (Intercept) AtBat Hits Walks CRuns
## 105.6487488 -1.9762838 6.7574914 6.0558691 1.1293095
## CWalks DivisionW PutOuts
## -0.7163346 -116.1692169 0.30288476개 변수 모델은 3개 모델이 모두 같으나 7개 모델은 모두 달라 부분집합에서의 최적 모델이 전진, 후진에서는 간과됨을 볼 수 있다. regsubsets에서는 단계별(stepwise)은 제공하지 않는다.
2. 검증셋 기법
추출된 훈련셋을 기준으로 각 변수 갯수별 최고의 모델의 검정 MSE값을 비교하자. !는 반대되는 벡터 선택을 의미하고 model.matrix는 X 행렬 %*%는 행렬 곱을 의미한다.
set.seed(2)
train=sample(c(TRUE,FALSE),nrow(Hitters),replace=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## [1] 133131.11 122794.01 121494.06 105461.51 111393.31 105881.94 107300.36
## [8] 107352.96 101524.62 99166.38 99856.40 99781.67 100000.87 100787.31
## [15] 101479.85 101508.55 101362.82 101299.74 101347.87which.min(val.errors)## [1] 10coef(regfit.best,10)## (Intercept) AtBat Hits Walks CAtBat
## 193.3398390 -2.0847371 5.9703838 6.3148144 -0.1893914
## CRuns CRBI CWalks DivisionW PutOuts
## 1.5580778 1.1626014 -0.8644379 -148.9749672 0.2496536
## Assists
## 0.504555710-folds와 달리 10개 변수가 최적 모델로 나온다. 물론 이 결과는 검증셋 기법의 특징상 표본을 다시 했을 경우 변동이 클 수 있다. 참고로 10-folds CV에서도 10개 변수 갯수는 거의 최적에 가깝다.
3. Adjusted R2, Cp, BIC Plot 법
하기 각 값에 대한 전체부분집합 plot은 사용할 일이 없을 것 같으나 많은 블로거들이 사용므로 구현하는 방법을 아는 차원에서 정리해본다. 검정색은 포함된 설명 변수를 나타낸다.
par(mfrow=c(1,1))
plot(regfit.full,scale="adjr2")
plot(regfit.full,scale="Cp")
plot(regfit.full,scale="bic")
마지막으로 Cpplot에 대한 명령어는 아래와 같다. car package를 불러온다.
library(car)
subsets(regfit.full, statistic="cp"")
abline(1,1,lty=2,col="red")'R-PAGO 노트' 카테고리의 다른 글
| 검정오차(test error)가 중요한 이유 (0) | 2017.03.30 |
|---|---|
| 회귀트리(regression tree)와 분류트리(classification tree) 예측 (0) | 2017.03.27 |
| Z검정, t검정, 카이제곱 검정 (0) | 2017.03.02 |
| UCLA 대학원 합격을 위한 필요 Spec.은? (0) | 2017.02.18 |
| 회귀모델의 변수 선택(Variable Selection) 방법 (0) | 2017.02.15 |
