머신러닝 예제

한치록(2022). 계량경제학강의, 제4판, 박영사, 19장 부록 “머신러닝 예제”, 버전 0.2.2.

이 문서에 제시된 R 명령들은 주로 James, Witten, Hastie and Tibshirani (2013, An Introduction to Statistical Learning: with Applications in R, Springer)을 참고하였으며, 그 책으로 충분하지 않은 부분은 패키지 문서에 기초하여 직접 코딩하였거나 인터넷 검색으로부터 얻은 정보를 확인하여 작성하였다.

회귀(regression) 분류(classifiction)
데이터 준비
Subset Selection
Splines
Ridge와 Lasso
PCR과 PLS
Decision Tree
Tree Ensemble
Support Vector Regression
Neural Networks
Super Learner
 데이터 준비
Logit, LDA, QDA
ROC와 Cutoff 값
Class Imbalance
Ridge와 Lasso
Decision Tree
Tree Ensemble
Support Vector Machines
Neural Networks
Super Leaner
실습용 데이터 준비는 여기 참조. 본 소절의 전체 코드는 여기에 있음.

Decision Tree

Decision tree를 실습해 보자. Decision tree는 각 구간(잎)마다 하나의 값만을 주므로 분류(classification)에 적절하다. 앞에서와 마찬가지로 다음과 같이 random over-sampling 기법으로 클래스 균형화를 하고 시작하자.

set.seed(1)
idx1 <- which(TrainSet$deny=='no')
idx2 <- with(TrainSet, sample(which(deny=='yes'), sum(deny=='no'), replace=TRUE))
Over <- TrainSet[c(idx1,idx2), ]
summary(Over$deny)
#   no  yes 
# 1948 1948

큰 나무

나무를 만들자. 우선 큰 나무를 만들어 본다. 아래 명령에서 mindev=.005 옵션(mindev의 디폴트 값은 0.01)은 잔 가지가 많은 나무를 만들도록 해 준다. mindev 값이 작을수록 나무를 크게 만든다.

library(tree) # install.packages("tree") if necessary
tr.full <- tree(deny~., data=Over, control=tree.control(nrow(Over), mindev=.005)) # large tree
plot(tr.full, col='royalblue')
text(tr.full, pretty = 0)
큰 나무

나뭇잎은 많다(27개).

가지치기

만들어진 나무(tr.full)를 그대로 사용하지 않고 CV를 통하여 가지치기를 해 주자. 이를 위해 우선 CV (10-fold)를 한다.

set.seed(1)
cv.tr <- cv.tree(tr.full, FUN = prune.tree) # K = 10
plot(cv.tr)
나무 CV

그림을 보면 잎(말단 노드)이 10개일 때 deviance가 가장 작다. 잎 10개 나무를 최적 나무로 하자.

tr.pruned <- prune.tree(tr.full, best=10)
plot(tr.pruned, col='royalblue')
text(tr.pruned, pretty = 0)
가지치기 완료된 나무

가지치기가 완료된 나무 그림을 보면 다음 변수들이 중요하다.

이 가지치기된 나무(tr.pruned)를 TrainSet 자료(Over 자료가 아니라)에 적용하여 예측하고, 성능을 살펴보자. 여기 써 먹으려고 데이터 준비 단원에서 Performance() 함수가 tree 객체를 처리할 수 있도록 이미 만들어 두었다.

Performance(tr.pruned, TrainSet) # defined in "데이터 준비"
# $ConfusionMatrix
#       pred
# actual   no  yes
#    no  1589  359
#    yes   68  197
# 
# $Summary
# Sensitivity Specificity   Precision    Accuracy 
#   0.7433962   0.8157084   0.3543165   0.8070493

전체 변수를 사용한 로짓(Over 데이터 사용)과 ROC 곡선을 비교하자.

logit.over <- glm(deny~., family = binomial, data = Over)
phat.logit <- predict(logit.over, TrainSet, type='r')
phat.ptree <- predict(tr.pruned, TrainSet)[,'yes']
library(ROCit) # install.packages("ROCit") if necessary
roc.logit <- rocit(phat.logit, TrainSet$deny)
roc.ptree <- rocit(phat.ptree, TrainSet$deny)
plot(roc.logit, YIndex=F)
with(roc.ptree, lines(TPR~FPR, col=2))
legend('right', c('Logit', 'Tree'), lwd=c(2,1), col=1:2, cex=.8, bty='n')
로짓과 나무의 ROC 곡선 비교

몇 가지 점이 눈에 띈다. 첫째, 로짓과 나무의 결과가 약간 다르다. 둘째, sensitivitiy와 specificity의 합의 최댓값은 나무의 경우가 더 크다. 셋째, 나무는 직선 여러 개를 합쳐 놓은 것처럼 보인다. 이는 나무가 ’piecewise constant’이기 때문이다. 위의 가지치기한 나무는 잎이 10개여서 unique한 확률예측값은 10개뿐이다.

table(phat.ptree)
# phat.ptree
# 0.171647509578544 0.323340471092077 0.333333333333333  0.65962441314554 
#              1113               338               206               184 
# 0.724137931034483 0.771428571428571 0.824039653035936 0.851063829787234 
#                22                11               226                11 
# 0.886904761904762 0.976377952755906 
#                82                20

그리고 roc.ptree가 고려하는 cutoff 값도 10개뿐이다.

roc.ptree$Cutoff
#  [1]       Inf 0.9763780 0.8869048 0.8510638 0.8240397 0.7714286 0.7241379
#  [8] 0.6596244 0.3333333 0.3233405 0.1716475

Cutoff가 0인 점은 좌하귀(0, 0)이고 cutoff가 1 (≥일 때 yes이므로 1이 아니라 Inf)인 점은 우상귀(1, 1)이다.

with(roc.ptree, cbind(Cutoff, TPR, FPR))
#          Cutoff        TPR         FPR
#  [1,]       Inf 0.00000000 0.000000000
#  [2,] 0.9763780 0.06415094 0.001540041
#  [3,] 0.8869048 0.23018868 0.021047228
#  [4,] 0.8510638 0.24528302 0.024640657
#  [5,] 0.8240397 0.56226415 0.097535934
#  [6,] 0.7714286 0.57358491 0.101642710
#  [7,] 0.7241379 0.59622642 0.109856263
#  [8,] 0.6596244 0.74339623 0.184291581
#  [9,] 0.3333333 0.79622642 0.282854209
# [10,] 0.3233405 0.87924528 0.445071869
# [11,] 0.1716475 1.00000000 1.000000000

이 가지치기된 나무를 test set에 적용할 때의 performance는 다음과 같다.

Performance(tr.pruned, TestSet) # defined in "데이터 준비"
# $ConfusionMatrix
#       pred
# actual  no yes
#    no  120  27
#    yes   8  12
# 
# $Summary
# Sensitivity Specificity   Precision    Accuracy 
#   0.6000000   0.8163265   0.3076923   0.7904192

Original train set 사용하고 cutoff 조정

Over-sample된 데이터가 아니라 원래 TrainSet 데이터를 사용하여 나무를 만들고 cutoff를 조정하는 방법을 사용해 보자. 원래 데이터에는 클래스 불균형이 현저함을 기억하라.

summary(Over$deny) # random over-sampling
#   no  yes 
# 1948 1948 
summary(TrainSet$deny) # original train set
#   no  yes 
# 1948  265

우선 큰 나무를 만들면 다음과 같다.

library(tree) # install.packages("tree") if necessary
tr.full.orig <- tree(deny~., data=TrainSet, control=tree.control(nrow(TrainSet), mindev=.005)) # large tree
plot(tr.full.orig, col='royalblue')
text(tr.full.orig, pretty = 0)
큰 나무(원래 자료 사용)

10-fold CV로 가지치기할 때의 deviance는 다음과 같다.

set.seed(1)
cv.tr.orig <- cv.tree(tr.full.orig, FUN = prune.tree) # K = 10
plot(cv.tr.orig)
나무 CV(원래 자료 사용)

나무 크기가 6개일 때 최적화된다. 가지치기된 나무는 다음과 같다.

tr.pruned.orig <- prune.tree(tr.full.orig, best = 6)
plot(tr.pruned.orig, col='royalblue')
text(tr.pruned.orig, pretty = 0)
가지치기 완료된 나무(원래 자료 사용)

나무 모양이 Over 자료를 사용할 때와 다르다. 이것으로 확률 0.5 기준으로 예측을 하면 train set에서 성능은 다음과 같다. 아마도 specificity는 높고 sensitivity는 낮을 것이다.

Performance(tr.pruned.orig, TrainSet) # defined in "데이터 준비"
# $ConfusionMatrix
#       pred
# actual   no  yes
#    no  1897   51
#    yes  174   91
# 
# $Summary
# Sensitivity Specificity   Precision    Accuracy 
#   0.3433962   0.9738193   0.6408451   0.8983281

예상한 대로 specificity는 높고 sensitivity는 낮다. 이는 0.5 임계값이 너무 높아서 좀처럼 yes로 예측을 하지 않기 때문이다. (그렇게 된 이유는 train set에 no가 대부분이어서 yes의 올바른 예측보다 no의 올바른 예측이 중요하기 때문이다.)

Youden Index를 최대화하는 cutoff 값을 구하자.

phat.ptree.orig <- predict(tr.pruned.orig, TrainSet)[,'yes']
library(ROCit) # install.packages("ROCit") if necessary
roc.ptree.orig <- rocit(phat.ptree.orig, TrainSet$deny)
(cutoff <- with(roc.ptree.orig, Cutoff[which.max(TPR-FPR)]))
# [1] 0.3291139

생각보다 값이 크다. 이 cutoff 값을 사용하여 예측하면 다음 결과를 얻는다.

Performance(tr.pruned.orig, TrainSet, cutoff = cutoff) # defined in "데이터 준비"
# $ConfusionMatrix
#       pred
# actual   no  yes
#    no  1791  157
#    yes  112  153
# 
# $Summary
# Sensitivity Specificity   Precision    Accuracy 
#   0.5773585   0.9194045   0.4935484   0.8784455

여전히 sensitivity가 낮다. 3개의 ROC를 그림으로 그려서 좀 더 살펴보자.

plot(roc.ptree.orig, YIndex=T, col=3)
with(roc.ptree, lines(TPR~FPR, col=2))
with(roc.logit, lines(TPR~FPR, col=1))
legend('right', c('Logit (oversampling)', 'Tree (oversampling)', 'Tree (original train set)'), lty=1, col=1:3, bty='n', cex=.75)
ROC 곡선의 비교

그림에 “Optimal (Youden Index) point”로 표시된 것은 original train set에 해당한다. 약간 왼쪽으로 치우친 느낌을 주는데, 이는 (우리 생각보다) cutoff가 높다는 것을 의미하고, 그로 인해서 sensitivity가 낮다. 이건 사람 눈으로 볼 때 그런 것이고, 아무튼 이를 test set에 적용하면 결과는 다음과 같다.

Performance(tr.pruned.orig, TestSet, cutoff = cutoff) # defined in "데이터 준비"
# $ConfusionMatrix
#       pred
# actual  no yes
#    no  133  14
#    yes  11   9
# 
# $Summary
# Sensitivity Specificity   Precision    Accuracy 
#   0.4500000   0.9047619   0.3913043   0.8502994
  
## --------------------------------------------------------------------
set.seed(1)
idx1 <- which(TrainSet$deny=='no')
idx2 <- with(TrainSet, sample(which(deny=='yes'), sum(deny=='no'), replace=TRUE))
Over <- TrainSet[c(idx1,idx2), ]
summary(Over$deny)
library(tree)
tr.full <- tree(deny~., data=Over, control=tree.control(nrow(Over), mindev=.005)) # large tree
plot(tr.full, col='royalblue')
text(tr.full, pretty = 0)
set.seed(1)
cv.tr <- cv.tree(tr.full, FUN = prune.tree) # K = 10
plot(cv.tr)
tr.pruned <- prune.tree(tr.full, best=10)
plot(tr.pruned, col='royalblue')
text(tr.pruned, pretty = 0)
Performance(tr.pruned, TrainSet)
logit.over <- glm(deny~., family = binomial, data = Over)
phat.logit <- predict(logit.over, TrainSet, type='r')
phat.ptree <- predict(tr.pruned, TrainSet)[,'yes']
library(ROCit)
roc.logit <- rocit(phat.logit, TrainSet$deny)
roc.ptree <- rocit(phat.ptree, TrainSet$deny)
plot(roc.logit, YIndex=F)
with(roc.ptree, lines(TPR~FPR, col=2))
legend('right', c('Logit', 'Tree'), lwd=c(2,1), col=1:2, cex=.8, bty='n')
table(phat.ptree)
roc.ptree$Cutoff
with(roc.ptree, cbind(Cutoff, TPR, FPR))
Performance(tr.pruned, TestSet)
summary(Over$deny) # random over-sampling
summary(TrainSet$deny) # original train set
library(tree)
tr.full.orig <- tree(deny~., data=TrainSet, control=tree.control(nrow(TrainSet), mindev=.005)) # large tree
plot(tr.full.orig, col='royalblue')
text(tr.full.orig, pretty = 0)
set.seed(1)
cv.tr.orig <- cv.tree(tr.full.orig, FUN = prune.tree) # K = 10
plot(cv.tr.orig)
tr.pruned.orig <- prune.tree(tr.full.orig, best = 6)
plot(tr.pruned.orig, col='royalblue')
text(tr.pruned.orig, pretty = 0)
Performance(tr.pruned.orig, TrainSet)
phat.ptree.orig <- predict(tr.pruned.orig, TrainSet)[,'yes']
library(ROCit)
roc.ptree.orig <- rocit(phat.ptree.orig, TrainSet$deny)
(cutoff <- with(roc.ptree.orig, Cutoff[which.max(TPR-FPR)]))
Performance(tr.pruned.orig, TrainSet, cutoff = cutoff)
plot(roc.ptree.orig, YIndex=T, col=3)
with(roc.ptree, lines(TPR~FPR, col=2))
with(roc.logit, lines(TPR~FPR, col=1))
legend('right', c('Logit (oversampling)', 'Tree (oversampling)', 'Tree (original train set)'), lty=1, col=1:3, bty='n', cex=.75)
Performance(tr.pruned.orig, TestSet, cutoff = cutoff)
## --------------------------------------------------------------------