ROC 곡선을 이용해 최적의 컷오프 찾는 법
개요
ROC 곡선을 그리면 트레이닝 데이터로 얻은 모형이 테스트 데이터를 어느정도로 잘 설명하는지 한 눈에 보여서 유용하다. 하지만 이 곡선은 모든 컷오프에 대한 분류율을 계산하고 이은 것이므로 결국 ‘어떤 컷오프로 0과 1을 분류할 것인가’는 알 수 없다. 이것을 알아내기 위해 교차검증의 방법론을 응용해보자.
검증 데이터
트레이닝 데이터에서 0과 1을 가장 잘 분류해주는 최적 컷오프를 찾아내봤자 그것은 트레이닝 데이터를 잘 설명할 컷오프에 지나지 않는다. 당연히 테스트 데이터에서 분류율을 계산해봐도 그 최적 컷오프는 테스트 데이터만을 잘 설명하는 컷오프일 뿐이다. 따라서 밸리데이션 데이터validation Data라는 것을 별도로 만들어 어떤 데이터에도 편향되지 않는 컷오프를 찾는다.
- Step 1.
트레이닝 데이터를 학습시켜 모형을 만든다. - Step 2.
밸리데이션 데이터에 적용시켜 최적의 컷오프를 찾는다. - Step 3.
트레이닝 데이터로 얻은 모형을 밸리데이션 데이터로 얻은 컷오프로 분류해 테스트 데이터에서 얼마나 퍼포먼스가 높은지 확인한다. 여러가지 모형 후보가 있다면 이 퍼포먼스가 가장 높은 모형을 최종 모형으로 선택한다.
이 때 Step 1에서 얻은 모형은 Step 2에서 얻은 컷오프를 포함하고서야 비로소 제대로 된 모형이 된다.
실습
(ROC 곡선 그리는 법에 이어서)
위 그림은 전 포스트에서 데이터를 트레이닝 데이터와 테스트 데이터로 나누었을 때의 ROC 곡선이다. 이제는 최적의 컷오프를 구하기 위해 데이터를 삼분할 필요가 있다. 이때의 비율은 데이터마다 다르지만 별다른 문제가 없다면 3:1:1 정도가 무난하다.
DATANUM<-nrow(DATA)
numeric(DATANUM)
DATANUM*c(0.6,0.2,0.2)
slicing<-sample(1:DATANUM)
slicing[slicing>(DATANUM*0.8)]<-(-1)
slicing[slicing>(DATANUM*0.6)]<-(-2)
slicing[slicing>0]<-(-3)
slicing<-slicing+4
train<-DATA[slicing==1,]; head(train)
valid<-DATA[slicing==2,]; head(valid)
test<-DATA[slicing==3,]; head(test)
데이터의 전처리를 끝내고 위의 코드를 실행하면 다음과 같이 트레이닝, 밸리데이션, 테스트 세 가지 데이터으로 쪼개진다.
Step 1.
트레이닝 데이터를 학습시켜 모형을 만든다.
out0<-glm(vote~.,family=binomial(),data=train); summary(out0)
vif(out0)
out1<-step(out0, direction="both"); summary(out1)
qchisq(0.95,df=1454)
vif(out1)
위의 코드를 실행하면 다음과 같이 변수선택절차를 거쳐 모형을 만들고 적합도 검정, 다중공선성을 체크해준다.
모형을 살펴보면 특별히 문제될 곳이 없다.
Step 2.
밸리데이션 데이터에 적용시켜 최적의 컷오프를 찾는다.
p <- predict(out1, newdata=valid, type="response")
pr <- prediction(p, valid$vote)
prf <- performance(pr, measure = "tpr", x.measure = "fpr")
win.graph(); plot(prf, main='ROC of Validation Data')
위의 코드를 실행하면 마치 밸리데이터 데이터를 테스트 데이터처럼 취급해서 다음과 같이 ROC 커브를 그려준다.
최적의 컷오프는 데이터와 목적에 따라 다르게 결정될 수 있지만, 별도의 주안점이 없다면 왼쪽 위의 $(0,1)$ 에서 가장 가까운 점을 찾아 그 점의 컷오프를 최적의 컷오프로 삼는다. 거리를 계산해야하기 때문에 코드는 다소 복잡하다.
optid<-(1:length(prf@y.values[[1]][-1]))[((prf@x.values[[1]][-1])^2 + (1-prf@y.values[[1]][-11])^2)
==min((prf@x.values[[1]][-1])^2 + (1-prf@y.values[[1]][-1])^2)]
points(prf@x.values[[1]][-1][optid],prf@y.values[[1]][-1][optid], col='red', pch=15)
optcut<-prf@alpha.values[[1]][-1][optid]; optcut
위의 코드를 실행하면 위의 설명대로 $(0,1)$ 에서 가장 가까운 점을 빨간색으로 표시해주고 그 지점의 컷오프를 출력해준다.
코드가 많이 복잡하지만 이해하려고 노력할 필요는 없다. 복잡한 것과 어려운 것은 다른 일이다. 위의 코드는 길기만하지 개념적으로는 전혀 어렵지 않다. 그냥 곡선 위의 모든 점과 $(0,1)$ 사이의 거리를 잰 후 그 거리가 가장 짧은 점을 선택한 것 뿐이다. 그 점에서의 $alpha.values를 참조하면 컷오프를 알 수 있다. 이렇게 구한 컷오프는 데이터를 가장 적절하게 잘 분류하는 컷오프로써 받아들일 수 있다. (다시 한 번 강조하지만, 이것은 절대적인 척도가 아니다. 사용자의 목적에 따라서 ‘최적의 컷오프’ 자체가 완전히 새롭게 정의될 수 있다.)
이 예제에서 얻은 최적의 컷오프는 $0.4564142$ 로써, 이보다 높으면 1로 판정하고 낮으면 0으로 판정하는게 가장 좋다고 받아들여도 무방하다.(벌써 세번째 강조하는건데, 가장 좋다고 받아들일 수 있는 것이지 가장 좋은 것이 아니다. 이에 대해 타당한 해석을 내놓는 것은 전적으로 분석자에게 달려있다.)
Step 3.
테스트 데이터에서 얼마나 잘 맞는지 확인해본다.
p <- predict(out1, newdata=test, type="response"); head(p,48)
table(test$vote, p>optcut)
위의 코드를 실행하면 테스트 데이터에서 확률을 계산해주고 최적 컷오프에 따른 오류행렬을 출력해준다.
위 오류행렬의 정분류율은 약 $81 \%$ 로써 꽤 쓸만하고, 분석자가 만족할만하다면 최종모형으로 받아들여봄직하다.눈치챘겠지만 엄밀히 말해 최적 컷오프를 구하는데 있어서 꼭 ROC 곡선을 그릴 필요는 없다. 어차피 계산을 위한 데이터는 데이터 프레임으로써 다 구해놨기 때문에 코드만 잘 돌려서 값만 얻어내도 전혀 상관 없다.
코드
아래는 예제 코드다.
install.packages("car")
install.packages("ResourceSelection")
install.packages("ROCR")
library(car)
library(ResourceSelection)
library(ROCR)
set.seed(150421)
?Chile
str(Chile)
nrow(Chile)
DATA<-na.omit(Chile)
DATA$vote[DATA$vote!='Y']<-'N'
DATA$vote<-factor(DATA$vote)
DATANUM<-nrow(DATA)
numeric(DATANUM)
DATANUM*c(0.6,0.2,0.2)
slicing<-sample(1:DATANUM)
slicing[slicing>(DATANUM*0.8)]<-(-1)
slicing[slicing>(DATANUM*0.6)]<-(-2)
slicing[slicing>0]<-(-3)
slicing<-slicing+4
train<-DATA[slicing==1,]; head(train)
valid<-DATA[slicing==2,]; head(valid)
test<-DATA[slicing==3,]; head(test)
out0<-glm(vote~.,family=binomial(),data=train); summary(out0)
vif(out0)
out1<-step(out0, direction="both"); summary(out1)
qchisq(0.95,df=1454)
vif(out1)
p <- predict(out1, newdata=valid, type="response")
pr <- prediction(p, valid$vote)
prf <- performance(pr, measure = "tpr", x.measure = "fpr")
win.graph(); plot(prf, main='ROC of Validation Data')
optid<-(1:length(prf@y.values[[1]][-1]))[((prf@x.values[[1]][-1])^2 + (1-prf@y.values[[1]][-11])^2)
==min((prf@x.values[[1]][-1])^2 + (1-prf@y.values[[1]][-1])^2)]
points(prf@x.values[[1]][-1][optid],prf@y.values[[1]][-1][optid], col='red', pch=15)
optcut<-prf@alpha.values[[1]][-1][optid]; optcut
p <- predict(out1, newdata=test, type="response"); head(p,48)
table(test$vote, p>optcut)