머신러닝 예제

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

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

Tree Bagging

회귀(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

randomForest 패키지를 사용하여 나무 bagging을 할 수 있다. 최대로 사용할 변수 개수(mtry)를 전체 예측변수 개수로 설정하면 된다. 아래에서 학습용 데이터(z14)를 이용하여 500회 부트스트랩을 한 나무 bagging을 실행하고 그 결과를 test set (z15)에 적용하여 예측성능을 구한다.

## Tree bag
library(randomForest) # install.packages("randomForest") if necessary
dim(z14)  # p = ncol(dataset) - 1
# [1] 223  20
set.seed(1)
tr.bag <- randomForest(ynext~., data=z14, mtry=19, importance=TRUE) # mtry=19 for treebag; 500 trees
RMSE(z15$ynext, predict(tr.bag, z15)) # bagging # RMSE() defined in index11.php
# [1] 52.20014
rmspe.rw # random walk, defined in index11.php
# [1] 53.24273

Test set에서 RMSE는 52.20014로, 앞의 가지치기된 나무 한 그루(RMSE = 71.22396)에 비하여 큰 개선이 이루어진 것을 볼 수 있다. Lasso (RMSE = 47.75)보다는 못하지만 임의보행 결과(53.24273)보다는 낫다. 변수 중요도는 다음과 같다.

importance(tr.bag)
#              %IncMSE IncNodePurity
# grdp       4.14643938      17881.77
# regpop     3.64682270      23244.42
# popgrowth -0.51396422      28830.06
# eq5d      -2.07628674      19738.18
# deaths    -0.73762568      21655.24
# drink      2.13062639      19275.66
# hdrink    -1.40168066      26605.05
# smoke     -0.46163688      20722.22
# aged      16.94874268    6834906.70
# divorce    4.00092212      30791.70
# medrate   -0.28777961      18528.55
# gcomp     -0.05347892      10686.13
# vehipc    -1.50337961      17421.78
# accpv     -1.40684012      12708.01
# dumppc    -0.36283830      19869.18
# stratio   -2.14823035      29466.98
# deathrate 44.31904874   17878644.80
# pctmale   -0.41336720      10890.20
# accpc      4.16383701      18820.55

위에 2가지 ‘변수 중요도’ 척도가 제시되어 있다. %IncMSE 열은 각 부트스트랩 표본에서 OOB (out of bag) 표본(부트스트랩 시 제외된 관측치들의 표본)의 다른 예측변수들은 그대로 두고 한 예측변수값만을 무작위로 뒤섞을 때 MSE가 증가하는 정도(모든 나무들에서 평균냄)를 증가분들의 표준편차로 나눈 값(표준편차가 0이면 나누지 않는데 이 경우에는 평균도 거의 늘 0임)이며, IncNodePurity는 나무에서 각 변수로써 분할이 이루어질 때 SSE가 감소하는 정도의 합계를 모든 나무들에 대하여 평균낸 값이다. 이상의 내용은 R 도움말과 Jake Hoarse의 블로그를 참조하였다.

위에 보고된 결과에 의하면 변수들 중 deathrate와 aged 변수의 중요도가 가장 높게 나온다. 다음과 같이 하여 그림을 그릴 수도 있다(그림은 별로 멋이 없다).

varImpPlot(tr.bag)

이상의 bagging에서는 튜닝 매개변수를 선택하지 않았다. 부트스트랩 나무 개수(B)는 많다고 하여 overfit이 되지 않으므로 튜닝 대상이 아니다. 가지치기도 하지 않았다. 매번 가지치기(pruning)를 한다면 가지치기 정도가 튜닝 대상이 되고, CV나 OOB를 이용하여 튜닝 매개변수를 결정할 수 있을 것이다.

참고로, randomForest 명령에서는 학습용 자료에 대한 설명력 지표를 제공한다. 부트스트랩 표본마다 SSR을 표본크기로 나눈 값(MSE)과 이로부터 구한 ‘R제곱 비슷한 지표’(pseudo R-squared)를 제공한다. B회 부트스트랩으로부터 구한 이 ’R제곱 비슷한 지표’는 다음과 같다. 학습자료 설명력은 매우 높음을 확인할 수 있다.

mean(tr.bag$rsq)
# [1] 0.9723799

Bagging의 효과

나무 한 그루와 달리 bagging을 하면 연속변수도 꽤 잘 맞추는 경향이 있다. 다음 예가 bagging의 효과를 극단적으로 보여 준다. 데이터를 다음과 같이 생성해 보자.

## Toy (quadratic function)
smpl <- data.frame(x=seq(-1,1,by=.01))
smpl$y <- 1-smpl$x^2  # y = 1-x^2
plot(y~x, data=smpl)

위 그림에 정확한 y = 1 − x² 관계가 보인다. 특성변수는 1개(x)이다.

나무는 piecewise-constant이므로 나무를 학습하는 것은 토막난 수평선을 그리는 것과 같다. 위 그림에서 토막난 수평선들을 그려 y를 맞추기는 쉽지 않아 보인다. 실제로 위 데이터에 대하여 하나의 tree로 학습하고 학습 데이터 예측치를 그려 보면 다음과 같다.

library(tree) # install.packages("tree") if necessary
tr <- tree(y~x, data=smpl)  # one tree
points(smpl$x, predict(tr, smpl), col='red')

앞서 말한 것처럼 tree는 piecewise constant이므로 여기저기 끊긴 수평선 모양의 점들이 보인다(빨간 색). 수평선들은 실제 y값의 형태를 대략적으로 따라가기는 하지만 제대로 맞추었다고 보기는 힘들다(조각 개수를 아주 많게 하면, 즉 나무 크기를 아주 크게 하면 더 잘 맞출 수는 있을 것이다).

Bagging을 하는 것은 토막난 수평선들을 여러 차례 그린 다음 각 x값마다 여러 수평선들 높이의 평균을 구하는 것과 같다. 많은 수평선들을 평균내면 곡선도 만들 수 있다. 실제 bagging으로 위 데이터를 학습하고 예측하면 다음 파란색과 같이 곡선을 잘 맞춘다(아래에서 예측변수는 1개이므로 mtry는 설정할 필요가 없다).

set.seed(1)
trbag <- randomForest(y~x, data=smpl, ntree = 1000)  # tree bag
points(smpl$x, predict(trbag, smpl), col='blue')

파란색 점들을 보면 tree bag의 효과를 확인할 수 있다. 단일 tree는 piecewise constant여서 제대로 맞추지 못하나 1000개 tree를 평균냄으로써 부드러운 곡선을 잘 맞추었다.

임의의 숲(Random Forest)

Bagging과 Random Forest의 차이는 나무를 만들 때 변수들을 모두 사용하느냐 아니면 고려할 변수들의 개수(mtry)를 제한하느냐이다.

우선 회귀나무 기본 m값인 p/3의 정수값을 사용한 임의의 숲 결과는 다음과 같다.

## Random Forest
set.seed(1)
tr.rf <- randomForest(ynext~., data=z14)
tr.rf
# 
# Call:
#  randomForest(formula = ynext ~ ., data = z14, importance = TRUE) 
#                Type of random forest: regression
#                      Number of trees: 500
# No. of variables tried at each split: 6
#
#           Mean of squared residuals: 3708.214
#                     % Var explained: 96.73

부트스트랩 표본 개수는 500개, 나뭇가지 분기 시 고려한 변수 개수(mtry)는 6개(19/3=6.333)이며, 이는 위 결과에도 표시되어 있다. 이 결과를 이용한 test set 예측 성능은 다음과 같다.

RMSE(z15$ynext, predict(tr.rf, z15, type="response")) # RMSE() defined in index11.php
# [1] 51.84617

OOB 관측치들을 이용하여 mtry를 CV할 수 있다. 이하에 나오는 X와 Y는 ‘데이터 준비’ 부분에서 만들었으므로 혹시 아래 코드가 작동하지 않으면 해당 페이지를 방문하여 ‘클립보드로 복사’ 버튼을 클릭한 후 R 콘솔에 붙여넣어 실행하고 나서 이 페이지로 돌아오기 바란다.

set.seed(1)
tuneRF(X,Y, ntreeTry = 500, stepFactor = 2, improve=.005) # X,Y: 데이터 준비 page
# mtry = 6  OOB error = 3708.214 
# Searching left ...
# mtry = 3  OOB error = 6063.54 
# -0.6351645 0.005 
# Searching right ...
# mtry = 12     OOB error = 3020.347 
# 0.1854981 0.005 
# mtry = 19     OOB error = 3174.927 
# -0.0511793 0.005 
#    mtry OOBError
# 3     3 6063.540
# 6     6 3708.214
# 12   12 3020.347
# 19   19 3174.927

OOB 관측치들을 이용하여 예측능력을 평가했을 때 mtry=12가 좋은 것으로 나타났으므로, 이를 이용하여 다시 임의의 숲 학습을 하자.

set.seed(1)
tr.rf2 <- randomForest(ynext~., data=z14, mtry=12, importance=TRUE)  # mtry=12

varImpPlot(tr.rf2)

위 학습 결과를 test set에 적용하면 예측 성과는 다음과 같다.

RMSE(z15$ynext, predict(tr.rf2, z15, type="response")) # RMSE() defined in index11.php
# [1] 49.29249

mtry를 아예 CV할 수도 있다. 이하 코드에서는 10-fold CV로 mtry를 튜닝한다. 시간이 걸리겠지만 mtry를 3~19에 대하여 CV로 탐색하고자 한다. CPU 코어가 여럿인 경우 병렬(parallel) 처리하면 더 빨리 할 수 있으므로 foreach 패키지와 doParallel 패키지를 이용해서 CV를 병렬로 진행하고자 한다. 우선 10-fold CV를 위한 그룹 구분을 해 준다.

set.seed(1)
group <- sample(1:10, nrow(z14), replace = TRUE)
table(group)
# group
#  1  2  3  4  5  6  7  8  9 10 
# 21 17 22 21 23 24 22 21 28 24

다음으로 병렬처리를 위한 함수를 만들고자 한다. 이 함수는 train set과 validation set을 입력 받아, 3~19까지의 mtry 값을 이용하여 train set에 대하여 임의의 숲 훈련을 하고 각 결과를 validation set에 적용한 예측치를 구하여 리턴해 준다. 행은 validation set에 대응하고 열은 mtry 값에 해당한다.

RFpredFun <- function(DF.train, DF.valid) {
  ans <- matrix(NA, nrow(DF.valid), 2) # do nothing for mtry=1,2
  rownames(ans) <- rownames(DF.valid)
  for (i in 3:19) { # mtry = 3, ..., 19
    set.seed(1)
    rf <- randomForest(ynext~., data=DF.train, mtry=i)
    ans <- cbind(ans, predict(rf, DF.valid))
  }
  ans  # Nx19
}

위 RFpredFun 함수를 RFpredFun(DF1,DF2)와 같이 사용하면 DF1 데이터셋을 이용하여 mtry를 3~19로 설정하여 임의의 숲 훈련을 하고 그 결과를 DF2 데이터셋에 적용하여 구한 예측치들의 행렬(행은 DF2의 행, 열은 1~19, 단 1~2는 NA로 설정)을 리턴한다.

library(foreach) # install.packages("foreach") if necessary
library(doParallel) # install.packages("doParallel") if necessary

## Prepare for parallel processing
(cores <- detectCores())
# [1] 8
cl <- makeCluster(cores[1]-1, outfile="")
registerDoParallel(cl)

## Parallel training and validation
cvrf <- foreach(fold = 1:10, .combine = rbind, .packages = 'randomForest') %dopar% {
  cat('Work for fold', fold, '\n')
  RFpredFun(z14[group!=fold,], z14[group==fold, ])
}

## Clean up
stopCluster(cl) # Don't forget this!

병렬(parallel) 처리를 해도 시간이 좀 걸릴 것이다. 순차(sequential) 처리 시 (CPU 코어가 많은 컴퓨터에서는) 이보다 훨씬 오랜 시간이 걸릴 것이다. 코딩의 명료함을 위해 시간을 희생할 각오가 되어 있다면 cvrf를 만드는 코드 전체를 다음과 같이 하여 순차적으로 처리해도 된다. 단, 시간이 꽤 오래 걸려 상당히 답답할 것이다.

cvrf <- NULL
for (fold in 1:10) {
  cat('Work for fold', fold, '\n')
  cvrf <- rbind(cvrf, RFpredFun(z14[group!=fold,], z14[group==fold,]))
}

위 코딩 결과 그룹 1~10 순서로 행들이 정렬되어(즉, 행들이 뒤섞여서) 원래 z14 데이터셋의 ynext 변수와 그대로 비교할 수 없으므로 행을 z14와 동일하도록 재정리해 주자. 이런 부분을 빠뜨리면 이해하기 어려운 결과가 나올 것이다.

cvrf <- cvrf[rownames(z14), ] # Important!

마지막으로 각 열(mtry 값에 해당)별로 RMSE를 구하고 비교하자.

rmse <- apply(cvrf, 2, function(x) RMSE(x, z14$ynext))
rmse
#  [1]       NA       NA 77.15646 69.51262 64.14057 61.30463 59.44525 57.50600
#  [9] 55.98502 55.50474 55.04596 54.67698 54.77896 54.73767 54.98873 55.41734
# [17] 55.40551 55.89370 55.98228
plot(rmse, type='o', ylab = 'RMSE from 10-fold CV')
which.min(rmse)
# [1] 12

운좋게도, 이렇게 탐색을 해도 mtry=12가 최적으로 판명되었다.

Gradient Boosting

gbm 패키지를 이용하여 gradient boosting을 하고자 한다. James, Witten, Hastie and Tibshirani (2013)을 참고하였다. 3가지 튜닝 매개변수가 있다. 하나는 나무를 업데이트하는 횟수 즉 반복 횟수(n.trees, 기본값은 100), 다른 하나는 나무 개선 시 반영하는 정도(shrinkage, 기본값은 0.1), 나머지 하나는 나무당 최대 나뭇가지의 개수(interaction.depth, 기본값은 1)이다. 일단 이 매개변수들을 (내키는 대로) 특정 값으로 설정한 gradient boosting을 해 보자.

## Gradient boosting
library(gbm) # install.packages("gbm") if necessary
set.seed(1)
boost <- gbm(ynext~., data=z14, distribution='gaussian', n.trees=1000, interaction.depth=4, shrinkage=0.05)
summary(boost)
#                 var    rel.inf
# deathrate deathrate 58.5745386
# aged           aged 34.3998587
# regpop       regpop  1.2023047
# grdp           grdp  0.7583550
# stratio     stratio  0.7421761
# medrate     medrate  0.6866381
# drink         drink  0.4681822
# deaths       deaths  0.4588440
# popgrowth popgrowth  0.4013245
# hdrink       hdrink  0.3827729
# smoke         smoke  0.3567991
# gcomp         gcomp  0.3235111
# accpc         accpc  0.2783915
# dumppc       dumppc  0.2754177
# eq5d           eq5d  0.1798527
# accpv         accpv  0.1773467
# vehipc       vehipc  0.1720970
# divorce     divorce  0.1615892
# pctmale     pctmale  0.0000000
RMSE(z15$ynext, predict(boost, z15)) # RMSE() defined in index11.php
# Using 1000 trees...
#
# [1] 57.50946

Random Forest에서와는 달리 Boosting에서는 반복할수록 training set의 설명력이 좋아진다. 이는 결국 overfit 문제가 발생함을 의미하므로 좋아만 할 일은 아니다. 나무 개선 횟수에 따라 training set의 error를 (로그 스케일로) 그려 보면 다음과 같다.

plot(boost$train.error, type='l', log='y')

predict에서 몇 번 boosting을 한 것을 이용하여 예측할지 옵션(n.trees)을 줄 수 있다. 기본값은 gbm에서 설정한 n.trees값과 동일하다.

이제 CV로써 최적 n.tree 값을 결정하자. 이를 위해서는 cv.fold 옵션을 주면 된다.

set.seed(1)
cv1 <- gbm(ynext~., data=z14, distribution='gaussian', n.trees=1000, interaction.depth=4, shrinkage=0.05, cv.folds=10)
(k <- gbm.perf(cv1))
# [1] 96
cv1$cv.error[k]  # square error loss
# [1] 3736.726

위에서는 10-fold CV로써 n.fold를 정한다. 실행 결과에 의하면 96회 boost한 결과에서 CV 오차가 최소화되고, 그 CV 오차는 3736.726이다. Training set에서의 잔차와 CV 예측오차를 그림으로 표현해 보면 CV의 의미를 잘 이해할 수 있다. (참고로, gbm.perf(cv1)라고 해도 그림을 얻지만, 로그 스케일로 그리기 위해서 별도 작업을 하였다.)

options(scipen = 100)
plot(cv1$train.error, type='l', log='y')  # CV error
lines(cv1$cv.error, lty=2)  # train error
abline(v=k, lty=3)
legend('topright', c('Training error', 'CV error'), lty=1:2)

Train error는 boosting을 반복할수록 점점 줄어들지만 CV error는 일정 횟수 이상이 되면 오히려 증가한다. Overfitting의 전형적인 모습이며, 나중에 살펴볼 인공신경망에서도 이와 유사한 일이 발생한다.

96회(세로 점선) boost한 결과를 이용한 test set에서의 예측 성과는 개선된 바가 없다.

RMSE(z15$ynext, predict(cv1, z15, n.trees=k)) # RMSE() defined in index11.php
# [1] 57.50219

위에서는 shrinkage를 0.05로 설정하였다. 이 값이 작은 것은 개선의 단계들을 더 ‘촘촘하게’ 훑고 지나간다는 뜻이다. 이 shrinkage를 0.1로 증가시키면 training error가 더 빨리 감소하고 최적 횟수는 더 작을 것이다.

set.seed(1)
cv2 <- gbm(ynext~., data=z14, distribution='gaussian', n.trees=1000, interaction.depth=4, shrinkage=0.1, cv.folds=10)
(k <- gbm.perf(cv2, plot=FALSE))
# [1] 37
cv2$cv.error[k]
# [1] 4060.88

최적 횟수(37번) boost한 CV 오차는 4060.88로서 shrink가 0.05일 때(cv1 참조)보다 더 나쁘다.

shrinkage를 0.001로 감소시키면 결과는 다음과 같다. 0.001이 아주 작은 값이어서 업데이트가 천천히 일어나므로 n.trees를 10000으로 설정하였다. 시간이 더 오래 걸릴 것이다.

set.seed(1)
cv3 <- gbm(ynext~., data=z14, distribution='gaussian', n.trees=10000, interaction.depth=4, shrinkage=0.001, cv.folds=10)
(k <- gbm.perf(cv3, plot=FALSE))
# [1] 4421
cv3$cv.error[k]
# [1] 3855.458

이 결과의 최적 CV 오차 정도는 3855.458로서 shrinkage = 0.05인 경우(cv1)보다 더 크다. Shrinkage를 0.01로도 해 보았으나 이 경우에도 최소 CV 오차가 cv1의 경우보다 더 큰 것으로 나타났다. (set.seed의 값을 다르게 하면 다른 결과를 얻을 것으로 예상된다.)

shrinkage를 0.05로 설정하고 interaction.depth를 2로 낮추면 결과는 다음과 같다.

set.seed(1)
cv4 <- gbm(ynext~., data=z14, distribution='gaussian', n.trees=1000, interaction.depth=2, shrinkage=0.05, cv.folds=10)
(k <- gbm.perf(cv4, plot=FALSE))
# [1] 118
cv4$cv.error[k]
# [1] 3731.242

interaction.depth가 4인 경우보다 CV 오차가 근소하게 더 낮다. 위의 cv4를 이용하여 k번(즉, 118번) 업데이트한 결과를 가지고 test set에 대하여 예측한 후 RMSE를 구하면 결과는 다음과 같다.

RMSE(z15$ynext, predict(cv4, z15, n.trees=k)) # RMSE() defined in index11.php
# [1] 56.00893

근소하게 개선되었으나 그리 인상적이지는 않다. 여기서도 gbm.perf(cv4)라고 하면 training error와 CV error 그림을 얻을 수 있다. 검정색은 training error, 녹색은 CV error, 파란색 세로 파선은 최적 반복횟수이다.

gbm.perf(cv4)
# [1] 118

참고로, gbm 패키지는 나무를 개선하면서 모든 관측치들을 사용하는 것이 아니라 임의로 추출된 100p% 관측치만을 사용한다. 여기서 p는 bag.fraction 옵션으로써 설정할 수 있고, 기본값은 0.5이다. gbm 명령에서 bag.fraction=1 옵션을 주면 나무 개선 단계에서 모든 관측치들을 사용한다. 이 bag.fraction도 튜닝 대상으로 삼을 수 있다(여기서는 해 보지 않겠다).

이하에서는 interaction.depth를 1~4로 하고 shrinkage를 0.01, 0.05, 0.1, 0.2 중 하나로 하는 모든 가능한 셋팅에 대하여 “grid search”를 해 본다. 단, bag.fraction은 1로 설정한다.

param <- expand.grid(depth=1:4, shrinkage=c(.01, .05, .1, .2))
param
#    depth shrinkage
# 1      1      0.01
# 2      2      0.01
# 3      3      0.01
# 4      4      0.01
# 5      1      0.05
# 6      2      0.05
# 7      3      0.05
# 8      4      0.05
# 9      1      0.10
# 10     2      0.10
# 11     3      0.10
# 12     4      0.10
# 13     1      0.20
# 14     2      0.20
# 15     3      0.20
# 16     4      0.20

위 16개 셋팅 각각에 대하여 최적 boosting 횟수를 CV로 결정하고 이 최적 횟수에서 얻어지는 CV 예측오차 정도를 비교하자. CPU 코어가 많은 경우 parallel processing을 하면 상당한 시간 절약이 가능하다.

## Prepare parallel processing
cores <- detectCores()
cl <- makeCluster(cores[1]-1, outfile="")
registerDoParallel(cl)

## Do parallel processing
gbmcv <- foreach(i = 1:nrow(param), .combine = rbind, .packages = 'gbm') %dopar% {
  set.seed(1)
  cv <- gbm(ynext~., data=z14, distribution='gaussian', n.trees = 1000, interaction.depth = param$depth[i], shrinkage = param$shrinkage[i], cv.folds=10, bag.fraction = 1)
  k <- gbm.perf(cv, plot = FALSE)
  cat('Setting', i, 'done\n')
  c(unlist(param[i,]), best.iter = k, cv.error = cv$cv.error[k])
}

## Clean up
stopCluster(cl) # Don't forget this

gbmcv <- as.data.frame(gbmcv)
gbmcv
#           depth shrinkage best.iter cv.error
# result.1  1     0.01      859       3658.911
# result.2  2     0.01      501       3723.234
# result.3  3     0.01      526       3829.116
# result.4  4     0.01      442       3812.664
# result.5  1     0.05      169       3648.469
# result.6  2     0.05      99        3754.657
# result.7  3     0.05      90        3812.563
# result.8  4     0.05      96        3763.043
# result.9  1     0.1       84        3754.808
# result.10 2     0.1       49        3809.053
# result.11 3     0.1       45        3754.504
# result.12 4     0.1       44        3865.514
# result.13 1     0.2       37        3785.591
# result.14 2     0.2       23        3775.115
# result.15 3     0.2       21        3671.347
# result.16 4     0.2       21        3818.903

(opt <- gbmcv[which.min(gbmcv$cv.error),])
#          depth shrinkage best.iter cv.error
# result.5     1      0.05       169 3648.469

위 결과에 의하면 interaction.depth는 1, shrinkage는 0.05인 경우이며, 이때 최적 개선(boosting) 횟수는 169이다. 이를 이용하여 test set에 대한 예측을 하면 결과는 다음과 같다.

#set.seed(1)  # unnecessary with bag.fraction = 1
gb <- gbm(ynext~., data=z14, distribution = 'gaussian', n.trees = opt$best.iter, interaction.depth = opt$depth, shrinkage = opt$shrinkage, bag.fraction = 1)
RMSE(z15$ynext, predict(gb, z15, type = 'response')) # RMSE() defined in index11.php
# Using 169 trees...
#
# [1] 52.96056

약간 개선이 되기는 했는데, 여기서 서술하지 않은 여러 실험 결과들에 의하면 이는 bag.fraction을 1로 설정했기 때문인 것 같다.

요약

필자의 경험에 의하면 gradient boosting은 train set의 학습에 엄청난 성과를 보이고, 그 결과 overfitting 문제가 쉽게 발생한다. 그래서 boosting 횟수(gbm 명령의 n.trees 매개변수) 튜닝이 핵심적이다(아래의 xgboost 패키지는 nrounds, h2o 패키지는 ntrees). 반면 random forest는 튜닝 매개변수에 덜 민감하여 사용자 입장에서 더 편안한 마음으로 사용할 수 있다. 각자 장단점이 있으며, 사용자는 유행에 휩쓸리기보다는 이것저것 다 해 보는 것이 좋다.

Extreme Gradient Boosting (XGBoost)

Gradient Boosting을 여러 면에서 개선하고, overfitting을 막기 위한 regularization을 도입한 것으로 Extreme Gradient Boosting (XGBoost)이 있다. XGBoost는 R에 xgboost 패키지로 구현되어 있다. Ian Johnson의 2020/11/29 블로그를 참고하였다. 실습할 데이터는 앞에서와 같다(필요하면 여기 방문).

xgboost의 튜닝 매개변수로 부스팅 횟수(nrounds), 나무의 maximum depth (max_depth, 기본값은 6), 개선속도(eta, 기본값은 0.3)가 있다. 그 외에도 gamma, min_child_weight, subsample 등 매개변수가 있으나 처음 3개 매개변수로 CV 대상을 제한한다.

우선 주어진 max_depth=4와 eta=.05에 대하여 최적 반복횟수를 찾아 보자. 최적 반복횟수를 찾을 때, 앞에서 gbm 패키지를 사용할 때에는 n.trees를 1000으로 설정하고, 가장 작은 CV error를 주는 반복횟수를 선택하였다. 이보다 약간 간편한 방법으로 CV error가 일정 단계 동안 개선되지 않으면 boosting을 중단하는 방법(early stopping)도 있다. 그렇게 하면 시간이 훨씬 단축된다. 다음 명령에서는 최대 1,000회 boosting을 반복하지만 5회(early_stopping_rounds) 이상 CV error의 개선이 이루어지지 않으면 중단(early stopping)하도록 설정한다. 참고로, early_stopping_rounds 값이 크면 stopping이 일어나기 더 어려우므로 최적 반복횟수가 더 커진다. 단, 최적 반복횟수는 난수의 시드만 바꾸어도 상당히 달라질 수 있으므로 결과에 예민하게 반응할 것은 없다.

## Extreme Gradient Boosting
library(xgboost) # install.packages("xgboost") if necessary
set.seed(1)
xgbcv <- xgb.cv(data=X, label=Y, nfold=10, nrounds=1000, early_stopping_rounds = 5, max_depth = 6, eta = .05, verbose=F)
xgbcv$best_iteration
# [1] 156
xgbcv$evaluation_log$test_rmse_mean[xgbcv$best_iteration]
# [1] 61.01081
with(xgbcv$evaluation_log, plot(iter, train_rmse_mean, type='l'))
with(xgbcv$evaluation_log, lines(iter, test_rmse_mean, lty=2))
abline(v=xgbcv$best_iteration, lty=3)

위 그림에서 실선은 training error, 파선은 CV error, 세로 점선은 ’best iteration’인 156회에 해당한다.

최적 반복횟수(156회)를 구하였으므로 다시 Gradient Boosting을 하고, test set에서 예측 성과를 확인하자.

xgb <- xgboost(data=X, label=Y, nrounds=xgbcv$best_iteration, max_depth=6, eta=.05, verbose=F)
RMSE(z15$ynext, predict(xgb, X15)) # RMSE() defined in index11.php
# [1] 59.51182

결과가 그리 좋지는 않다(xgboost 패키지가 안 좋다는 뜻이 아니다).

이제 여러 max_depth와 nrounds 값들을 CV로 비교해 보자. 각각의 주어진 매개변수 셋팅에서 nrounds는 early stopping으로 찾는다. 고려할 매개변수 값들은 다음과 같다.

param <- expand.grid(max_depth = 2:6, eta = c(.01, .05, .1, .2))
param
#    max_depth  eta
# 1          2 0.01
# 2          3 0.01
# 3          4 0.01
# 4          5 0.01
# 5          6 0.01
# 6          2 0.05
# 7          3 0.05
# 8          4 0.05
# 9          5 0.05
# 10         6 0.05
# 11         2 0.10
# 12         3 0.10
# 13         4 0.10
# 14         5 0.10
# 15         6 0.10
# 16         2 0.20
# 17         3 0.20
# 18         4 0.20
# 19         5 0.20
# 20         6 0.20

위 20가지 매개변수 셋팅에 대하여 최적 CV 오차 RMSE를 구하자(’최적’은 nrounds를 최적으로 선택한다는 뜻).

## Prepare parallel processing
cores <- detectCores()
cl <- makeCluster(cores[1]-1, outfile="")
registerDoParallel(cl)

## Do parallel processing
xgbcv.full <- foreach(i = 1:nrow(param), .combine = rbind, .packages = 'xgboost') %dopar% {
  set.seed(1)
  cv <- xgb.cv(data=X, label=Y, nfold=10, nrounds=5000, early_stopping_rounds=5, max_depth=param$max_depth[i], eta=param$eta[i], verbose = F)
  best.iter <- cv$best_iteration
  cv.error <- cv$evaluation_log$test_rmse_mean[cv$best_iteration]
  cat('Setting', i, 'done\n')
  c(unlist(param[i,]), best.iter = best.iter, cv.error = cv.error)
}

## Clean up
stopCluster(cl) # Don't forget this

xgbcv.full <- as.data.frame(xgbcv.full) # not necessary but convenient
xgbcv.full
#           max_depth  eta best.iter cv.error
# result.1          2 0.01       562 55.43059
# result.2          3 0.01       539 58.76363
# result.3          4 0.01       532 61.01573
# result.4          5 0.01       540 61.78269
# result.5          6 0.01       603 61.92822
# result.6          2 0.05       109 55.47342
# result.7          3 0.05       101 59.00048
# result.8          4 0.05        98 61.85777
# result.9          5 0.05       113 61.91448
# result.10         6 0.05       156 61.01081
# result.11         2 0.10        54 55.16352
# result.12         3 0.10        52 58.89465
# result.13         4 0.10        51 60.35394
# result.14         5 0.10        62 61.28542
# result.15         6 0.10        60 60.96928
# result.16         2 0.20        27 56.40606
# result.17         3 0.20        25 59.15859
# result.18         4 0.20        30 62.68687
# result.19         5 0.20        30 63.62557
# result.20         6 0.20        31 64.89133

(opt <- xgbcv.full[which.min(xgbcv.full$cv.error),])
#           max_depth eta best.iter cv.error
# result.11         2 0.1        54 55.16352

이 셋팅을 이용하여 Gradient Boosting 하고 예측하면 다음과 같다.

xgb <- xgboost(data = X, label = Y, nrounds = opt$best.iter, max_depth = opt$max_depth, eta = opt$eta, verbose = F)
RMSE(z15$ynext, predict(xgb, X15)) # RMSE() defined in index11.php
# [1] 51.74302

’체면’은 차렸다. 변수 중요도는 다음과 같이 시각화해 볼 수 있다.

xgb.plot.importance(xgb.importance(model = xgb))

참고로, xgboost에서도 나무를 개선할 때 일부만 사용하도록 할 수 있다. subsample 옵션을 사용하는데, 기본값은 1이다. 만약 subsample을 0.5 등으로 바꾸면 xgboost를 실행할 때마다 결과가 달라지므로 결과를 재연할 수 있도록 하려면 set.seed(1) 등으로 통제해야 할 것이다.

H2O 패키지를 이용한 Tree 앙상블

앞에서 소개한 h2o 패키지(설치)를 이용하여 tree 앙상블을 구현해 보자. 이 패키지는 데이터셋 크기가 클 때 특히 유용할 것으로 추측한다. h2o 패키지는 통일된 인터페이스를 제공하므로 한 번 익숙해지면 그 다음부터는 편리하다. 앞에서 lasso를 H2O로 구현해 본 바 있다.

## H2O package
library(h2o) # install.packages("h2o") if necessary
h2o.init()  # does lots of things
z14h <- as.h2o(z14)  # training set
z15h <- as.h2o(z15)  # test set
yvar <- 'ynext'
xvar <- setdiff(names(z14), yvar)  # character vector

Random Forest 학습을 해 보자. Boosting과 달리 부트스트랩은 반복횟수(ntrees)가 많아도 overfit 문제를 야기하지 않으므로 ntrees 매개변수와 관련된 CV를 하지 않는다. 그냥 큰 값을 주면 된다. CV는 max_depth나 mtries 매개변수를 대상으로 grid search 방식으로 할 수 있다. 전부 기본값을 준 상태에서 ntrees = 500으로 설정하여 학습을 하고 변수 중요도를 그림으로 그리면 다음과 같다.

rf.h2o <- h2o.randomForest(xvar, yvar, z14h, ntrees = 500, seed = 1)
h2o.varimp_plot(rf.h2o)  # only first 10 (give option to view more)

위에서 변수 중요도는 10개 변수까지만 그렸다. 변수 개수를 바꾸고 싶으면 h2o.varimp_plot(rf.h2o, 19)처럼 옵션을 주면 된다.

OOB training sample들의 RMSE는 rf.h2o라고 하면 바로 화면에 표시할 수 있고, h2o.rmse(rf)라고 하면 숫자를 구할 수 있다.

rf.h2o
# Model Details:
# ==============
#
# H2ORegressionModel: drf
# Model ID:  DRF_model_R_1632619166365_5 
# Model Summary: 
#   number_of_trees number_of_internal_trees model_size_in_bytes min_depth
# 1             500                      500              912692        11
#   max_depth mean_depth min_leaves max_leaves mean_leaves
# 1        18   13.26200        123        161   140.72200
#
#
# H2ORegressionMetrics: drf
# ** Reported on training data. **
# ** Metrics reported on Out-Of-Bag training samples **
#
# MSE:  3528.036
# RMSE:  59.39727
# MAE:  42.41623
# RMSLE:  0.07142407
# Mean Residual Deviance :  3528.036
#
h2o.rmse(rf.h2o)  # OOB
# [1] 59.39727

Test set에 학습 결과를 적용하여 예측 성능을 구하면 다음과 같다.

h2o.performance(rf.h2o, newdata = z15h)
# H2ORegressionMetrics: drf
# 
# MSE:  2596.366
# RMSE:  50.95455
# MAE:  36.51445
# RMSLE:  0.06642564
# Mean Residual Deviance :  2596.366

RMSE는 50.95455로서 앞에서 randomForest 패키지를 활용한 경우와 유사하다. 참고로, mtries를 12로 설정하면 결과는 다음과 같다.

rf2.h2o <- h2o.randomForest(xvar, yvar, z14h, ntrees = 500, mtries = 12, seed = 1)
h2o.performance(rf2.h2o, newdata = z15h)
# H2ORegressionMetrics: drf
# 
# MSE:  2520.726
# RMSE:  50.20683
# MAE:  35.97571
# RMSLE:  0.06004474
# Mean Residual Deviance :  2520.726

Random number가 발생되는 방식이 달라서 randomForest 패키지의 경우(RMSE는 49.29249)와 근소하게 다르다.

Gradient Boosting을 하려면 h2o.gbm 명령을 사용한다. Boosting은 boost할수록 학습성과가 좋아지므로 ntrees가 너무 크면 overfitting 문제가 발생한다. ntrees가 중요한 튜닝 매개변수이므로 잘 선택해 주어야 한다. 한 가지 방법은 다음과 같이 CV로써 찾는 것이다.

b0.h2o <- h2o.gbm(xvar, yvar, z14h, nfolds = 10, ntrees = 250, score_each_iteration = TRUE, max_depth = 2, learn_rate = 0.1, seed = 1)
h2o.learning_curve_plot(b0.h2o)

이렇게 하기는 했으나 CV error를 가장 작게 해 주는 반복횟수는 별도로 구해야 한다. 다음이 필자가 찾은 방법인데 더 쉬운 방법이 있는지는 모르겠다.

cv.mse <- sapply(b0.h2o@model$cv_scoring_history, function(x) x$validation_rmse^2)
dim(cv.mse)  # (iter = 0, 1, ..., 250) x (nfolds = 10)
# [1] 251   10
(n <- which.min(rowMeans(cv.mse))-1)
# [1] 48

위에서 48개라고 하였으므로 48개 tree를 만드는 Gradient Boosting을 학습하고 test set에 적용하면 다음 결과를 얻는다.

b.h2o <- h2o.gbm(xvar, yvar, z14h, ntrees = n, max_depth = 2, learn_rate = 0.1, seed = 1)
h2o.rmse(b.h2o, train = TRUE)  # train set
# [1] 38.87085
h2o.performance(b.h2o, newdata = z15h)  # test set
# H2ORegressionMetrics: gbm
# 
# MSE:  2792.558
# RMSE:  52.84466
# MAE:  37.81112
# RMSLE:  0.06441821
# Mean Residual Deviance :  2792.558

Early stopping을 사용할 수도 있다. 아래 명령에서는 stopping_rounds 옵션을 ’3’으로 주어, boost 시 최근 ’3’회 CV 에러들의 단순 이동평균 감소 정도가 ’3’회 이상 stopping_tolerance (기본값은 0.001)에 미달하면 boosting을 멈추도록 한다(h2o.gbm 도움말 참조). 이 방법으로 구한 boosting 횟수(나무 개수)는 40회이다(아래 결과 참조). 참고로, stopping_rounds가 클수록, 그리고 stopping_tolerance가 클수록 early stopping이 일어나기 어려우므로 결과적으로 선택되는 boosting 횟수(나무 개수)가 많아진다(overfitting 가능성 ↑).

b2.h2o <- h2o.gbm(xvar, yvar, z14h, nfolds = 10, ntrees = 1000, score_each_iteration = TRUE, max_depth = 2, learn_rate = 0.1, stopping_rounds = 3, seed = 1)
b2.h2o@model$model_summary
# Model Summary: 
#   number_of_trees number_of_internal_trees model_size_in_bytes min_depth
# 1              40                       40                4325         2
#   max_depth mean_depth min_leaves max_leaves mean_leaves
# 1         2    2.00000          3          4     3.95000
h2o.learning_curve_plot(b2.h2o)

h2o.performance(b2.h2o, newdata = z15h)
# H2ORegressionMetrics: gbm
# 
# MSE:  2754.549
# RMSE:  52.4838
# MAE:  37.70973
# RMSLE:  0.06577292
# Mean Residual Deviance :  2754.549

위에서는 max_depth와 learn_rate를 2와 0.1로 설정한 상태에서 10-fold CV에 의하여 최적의 ntrees (boosting 횟수, 혹은 나무 개수)를 선택하였다. max_depth와 learn_rate도 튜닝 매개변수인데, 이 값들을 바꾸어 가면서 각 셋팅에서 최소 CV MSE를 구하고(early stopping 없이 수동으로 CV MSE를 구했다면 최소 CV MSE 사용, early stopping을 사용했다면 h2o.mse(..., xval = TRUE) 또는 h2o.rmse(..., xval = TRUE) 사용), 이 최소 CV MSE가 가장 작은 셋팅을 선택하는 grid search도 가능하다.