Tree-based Ensemble

stat17_hb @ korea.ac.kr

이 글은 고려대학교 강필성 교수님의 2018년 2학기 비즈니스 어낼리틱스 강의 내용을 바탕으로 작성되었습니다.

1 Random Forests

1.1 Intro

랜덤 포레스트는 트리들의 상관성을 제거하는(decorrlate) 단순한 방법으로 배깅 트리보다 더 좋은 성능을 보여준다. 트리들의 상관성을 제거하는 이유는 모델의 다양성을 증가시키기 위함이다. 랜덤 포레스트에서는 모델의 다양성(diversity)을 증가시키기 위해 두 가지 방법을 사용한다. 하나는 배깅이고, 다른 하나는 각 분할(split)에서 사용할 예측변수들을 랜덤하게 선택하는 것이다. 배깅은 부트스트랩(bootstrap)을 통해 각각의 트리를 만들때 사용되는 인스턴스의 다양성을 증가시킨다. 각 분할에서 사용할 예측변수들을 랜덤하게 선택하는 것은 모델에 사용되는 변수의 다양성을 증가시키는 효과를 준다. 이 글에서는 모델의 다양성을 증가시키기 위해 랜덤포레스트가 사용하는 방법들을 살펴보고, Out of Bag, Variable importance와 관련된 내용들을 정리해보고자 한다.

1.2 Algorithm

랜덤 포레스트 알고리즘([1] p.588)은 다음과 같다:

1.3 Diversity in Ensemble Model

모델의 다양성을 증가시키기 위한 두 가지 방법에 대해 더 자세하게 살펴보자. 배깅의 핵심 아이디어는 모델의 분산을 줄이기 위해 분산이 크지만 편향(bias)이 작은 모델들을 많이 만들어서 평균을 내는 것이다. 트리 알고리즘은 모델의 분산이 큰 알고리즘 중 하나이기 때문에 평균을 내는 것이 매우 효과적이다. 게다가, 배깅을 통해 생성된 트리 각각은 동일한 분포를 따르기 때문에(identically distributed) B개 트리들의 평균(sample average)의 기댓값(expectation)은 각각의 트리의 기댓값과 같다. 이 말은 배깅된 트리들의 편향은 개별 트리의 편향과 같고, 따라서 모델을 개선시킬 수 있는 유일한 방법은 분산을 줄이는 것이라는 의미이다. 반면, 부스팅에서는 트리들이 적응적인(adaptive) 방식으로 편향을 제거하면서 커지기 때문에 동일한 분포를 따르지 않는다.

만약 B개의 트리들이 독립적이고 동일한 분포를 따르고(i.i.d.; independently and identically distributed), 각 트리들의 분산이 \(\sigma^2\)이면 B개 트리들의 평균의 분산은 \(\frac{1}{B} \sigma^2\)이 될 것이다. 하지만 B개의 트리들이 독립적이지 않다고 하면, B개 트리들의 평균의 분산은 \(\rho\sigma^2+\frac{1-\rho}{B}\sigma^2\)이 된다. 여기서 \(\rho\)는 트리 쌍들의 상관계수이다. B가 커지면 두 번째 항이 매우 작아지게 되겠지만 여전히 첫 번째 항이 남아있게 된다. 즉, 평균을 내는 것을 통해 분산이 얼마나 줄어들지는 트리 쌍들의 상관계수에 따라 달라진다. 따라서 트리 쌍들의 상관계수를 줄일 수 있다면 앙상블 모델 전체의 분산도 줄일 수 있을 것이다.

배깅에서는 각 분할에 p개의 예측변수를 모두 사용한다. 만약 데이터에 매우 강한 설명변수가 하나 있고, 다수의 적당히 강한 설명변수가 있다고 해보자. 그러면 모든 부트스트랩 데이터세트의 맨 위 분할(top split)에 가장 강한 설명 변수가 사용될 것이다. 따라서 배깅된 트리들은 모두 서로 상당히 유사할 것이고, 결과적으로 단일 트리에 비해 모델의 분산을 크게 줄이지 못할 것이다. 랜덤포레스트는 각 분할에서 랜덤하게 선택된 p보다 작은 m개의 예측변수를 사용하는 것을 통해 트리들의 상관성을 줄인다. 이를 통해 앙상블 모델의 분산을 배깅에서보다 크게 줄 일 수 있게 된다.

일반적으로 m의 기본 값은 분류문제에서는 \(\sqrt{p}\), 회귀문제에서는 \(p/3\)이다(랜덤 포레스트를 만든 Leo Breiman과 Adele Cutler가 추천하는 값이다). 랜덤 포레스트가 m에 민감하지 않기 때문에 굳이 m에 대한 미세조정(fine tuning)을 할 필요는 없다고 한다 [2]. 하지만 튜닝 자체를 하지 않아도 되는 것은 아니다. 문제에 따라 m의 기본값을 사용할 때 성능이 좋지 않은 경우가 있을 수 있기 때문이다.

이를 Boston{MASS} 데이터를 통해 살펴보자.

library(MASS)
data(Boston)
str(Boston)
## 'data.frame':    506 obs. of  14 variables:
##  $ crim   : num  0.00632 0.02731 0.02729 0.03237 0.06905 ...
##  $ zn     : num  18 0 0 0 0 0 12.5 12.5 12.5 12.5 ...
##  $ indus  : num  2.31 7.07 7.07 2.18 2.18 2.18 7.87 7.87 7.87 7.87 ...
##  $ chas   : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ nox    : num  0.538 0.469 0.469 0.458 0.458 0.458 0.524 0.524 0.524 0.524 ...
##  $ rm     : num  6.58 6.42 7.18 7 7.15 ...
##  $ age    : num  65.2 78.9 61.1 45.8 54.2 58.7 66.6 96.1 100 85.9 ...
##  $ dis    : num  4.09 4.97 4.97 6.06 6.06 ...
##  $ rad    : int  1 2 2 3 3 3 5 5 5 5 ...
##  $ tax    : num  296 242 242 222 222 222 311 311 311 311 ...
##  $ ptratio: num  15.3 17.8 17.8 18.7 18.7 18.7 15.2 15.2 15.2 15.2 ...
##  $ black  : num  397 397 393 395 397 ...
##  $ lstat  : num  4.98 9.14 4.03 2.94 5.33 ...
##  $ medv   : num  24 21.6 34.7 33.4 36.2 28.7 22.9 27.1 16.5 18.9 ...

각 분할에서 사용될 변수의 개수는 mtry인자에 명시해주면 된다. medv를 타겟변수로 나머지 13개의 변수를 예측변수로 한 후, 다른 옵션은 모두 기본 값으로 두었다. 5-fold Cross-Validation을 사용하여 CV error를 계산하였다.

library(randomForest)
## randomForest 4.6-14
## Type rfNews() to see new features/changes/bug fixes.
set.seed(1)
folds <- sample(1:5, size = nrow(Boston), replace = T) #5 fold CV
grid.mtry <- c(2,4,6,8,10,12)
cv.error <- matrix(0, nrow=length(grid.mtry), ncol=5)

for (j in 1:5){
  print(paste0(j, "th folds"))
  for (i in 1:length(grid.mtry)){
    cv.rf <- randomForest(medv~., mtry=grid.mtry[i], data=Boston[folds != j, ])
    cv.pred <- predict(cv.rf, newdata=Boston[folds == j,])
    cv.error[i,j] <- mean((Boston$medv[folds == j]-cv.pred)^2)
  }
}
## [1] "1th folds"
## [1] "2th folds"
## [1] "3th folds"
## [1] "4th folds"
## [1] "5th folds"
plot(grid.mtry, rowMeans(cv.error), type="l", 
     xlab="Number of predictors for each split (mtry)", ylab="CV error")

mtry의 기본 값은 max(floor(ncol(x)/3), 1)=4인데, 8일 때의 CV error가 가장 낮은 것을 볼 수 있다.

1.4 Out of Bag Samples

배깅된 트리는 평균적으로 관측치들의 약 2/3을 이용한다. Out of Bag(이하 OOB) 샘플은 배깅된 트리를 적합하는데 사용되지 않은 나머지 1/3의 관측치들이다(이에 관한 설명은 여기에 정리해 두었다).

랜덤 포레스트는 배깅을 기반으로 하기 때문에 OOB 샘플을 test set 대신 사용하여 OOB error를 구할 수 있다. OOB error는 test error와 거의 비슷한 값을 준다고 알려져 있다. 이를 Boston 데이터를 통해 살펴보자.

library(randomForest)
set.seed(1)
fit.rf <- randomForest(medv~., data=Boston)
fit.rf
## 
## Call:
##  randomForest(formula = medv ~ ., data = Boston) 
##                Type of random forest: regression
##                      Number of trees: 500
## No. of variables tried at each split: 4
## 
##           Mean of squared residuals: 9.814389
##                     % Var explained: 88.37

여기서 나오는 Mean of squared residuals(MSE; Mean Squared Error) 값이 OOB sample을 통해 계산된 OOB error이다. 이 값은 아래 코드를 통해서도 확인해 볼 수 있다.

fit.rf$mse[500]
## [1] 9.814389

fit.rf$mse에는 트리 개수를 증가시키면서 구한 OOB error 값들이 저장되어 있는데 plot 함수에 fit.rf를 넣으면 이에 대한 그래프를 보여준다.

plot(fit.rf, main="OOB Error")

# plot(1:500, fit.rf$mse, xlab="trees", ylab="Error", main="fit.rf", type="l")

이제 데이터세트를 train과 test로 분할 한 후, test MSE를 계산해서 OOB error와 비교해보자. 데이터세트 분할 비율은 OOB 샘플의 비율로 하였다(전체 샘플의 2/e). 트리 개수는 ntree 인자에 명시해주면 된다.

set.seed(1)
test.idx <- sample(nrow(Boston), nrow(Boston)*(2/exp(1)))
train <- Boston[-test.idx, ]
test <- Boston[test.idx, ]

test.mse <- c()
for (i in 1:500){
  if (i %% 100 == 0) print(paste0("ntree=", i))
  fit.rf2 <- randomForest(medv~., ntree=i, data=train)
  rf.pred <- predict(fit.rf2, newdata = test)
  test.mse[i] <- mean((test$medv-rf.pred)^2)
}
## [1] "ntree=100"
## [1] "ntree=200"
## [1] "ntree=300"
## [1] "ntree=400"
## [1] "ntree=500"
plot(fit.rf, main="")
lines(1:500, test.mse, col="red")
legend("topright", legend=c("OOB error", "Test error"), col=c("black","red"), lwd=2)

Test error가 OOB error보다 큰 경향을 보인다. Cross-Validation error를 구해서 비교해보자. 여기서는 5-fold CV를 사용하였다.

set.seed(1)
folds <- sample(1:5, size = nrow(Boston), replace = T) #5 fold CV
cv.error <- matrix(0, nrow=500, ncol=5)

for (j in 1:5){
  for (i in 1:500){
    if (i %% 100 == 0) print(paste0(j, "th folds ", "ntree=", i))
    cv.rf <- randomForest(medv~., ntree=i, data=Boston[folds != j, ])
    cv.pred <- predict(cv.rf, newdata=Boston[folds == j,])
    cv.error[i,j] <- mean((Boston$medv[folds == j]-cv.pred)^2)
  }
}
## [1] "1th folds ntree=100"
## [1] "1th folds ntree=200"
## [1] "1th folds ntree=300"
## [1] "1th folds ntree=400"
## [1] "1th folds ntree=500"
## [1] "2th folds ntree=100"
## [1] "2th folds ntree=200"
## [1] "2th folds ntree=300"
## [1] "2th folds ntree=400"
## [1] "2th folds ntree=500"
## [1] "3th folds ntree=100"
## [1] "3th folds ntree=200"
## [1] "3th folds ntree=300"
## [1] "3th folds ntree=400"
## [1] "3th folds ntree=500"
## [1] "4th folds ntree=100"
## [1] "4th folds ntree=200"
## [1] "4th folds ntree=300"
## [1] "4th folds ntree=400"
## [1] "4th folds ntree=500"
## [1] "5th folds ntree=100"
## [1] "5th folds ntree=200"
## [1] "5th folds ntree=300"
## [1] "5th folds ntree=400"
## [1] "5th folds ntree=500"
plot(fit.rf, main="")
lines(1:500, test.mse, col="red")
lines(1:500, rowMeans(cv.error), col="blue")
legend("topright", legend=c("OOB error", "Test error", "CV error"), col=c("black","red", "blue"), lwd=2)

CV error와 OOB error와 거의 비슷한 경향을 보이는 것을 알 수 있다.

1.5 Variable Importance

변수 중요도를 계산하려면 importance 인자를 TRUE로 설정해주면 된다.

set.seed(1)
fit.rf <- randomForest(medv~., importance=T, data=Boston)

변수 중요도 값은 importance() 함수를 사용해서 확일 할 수 있다.

importance(fit.rf)
##           %IncMSE IncNodePurity
## crim    16.393291     2549.9030
## zn       5.103123      225.6366
## indus   11.487445     2351.3389
## chas     4.955784      246.3893
## nox     19.353551     2603.5509
## rm      36.375678    12089.9862
## age     12.772651     1126.2815
## dis     17.927151     2478.4243
## rad      6.840204      352.7493
## tax     12.818632     1420.6800
## ptratio 14.758557     2685.9002
## black    9.883777      751.3808
## lstat   31.069552    13264.0767

두 가지의 변수 중요도가 계산되는데 %IncMSE는 OOB 데이터를 permuting 해서 계산되는 값이다. 각각의 트리에 대해, OOB 데이터에 대한 prediction error(분류인 경우 error rate, 회귀인 경우 MSE)를 계산한다. 그런 다음 같은 작업을 각 변수들의 순서를 섞어서(permutation) 진행한다. 이렇게 계산된 두 개의 prediction error의 차이를 모든 트리에 대해 평균을 낸 후 차이의 표준편차로 표준화를 해서 계산된 값이 %IncMSE이다.

IncNodePurity는 해당 변수가 분할에 사용되었을 때 감소된 노드 불순도(node impurity)의 총량을 모든 트리에 대해 평균 낸 것이다. 노드 불순도는 분류인 경우 지니 지수(Gini index)로 측정되고, 회귀인 경우 잔차 제곱합(residual sum of square)로 측정된다.

varImpPlot() 함수는 변수 중요도를 그래프로 표현해준다.

varImpPlot(fit.rf, main="Variable Importance")

2 Rotation Forest

2.1 Intro

로테이션 포레스트는 특징 추출(feature extraction)에 기반한 앙상블 기법이다. base learner가 사용할 훈련 데이터(training data)를 만들기 위해, 예측변수들을 K개의 부분집합(subset)으로 랜덤하게 분할한 후 주성분 분석(PCA; Principal Component Analysis)를 각각의 부분집합에 적용한다. 데이터의 변동성(variablility) 정보를 보존하기 위해 모든 주성분(Principal Components)들이 사용된다. 결과적으로, base learner가 사용할 새로운 예측변수들을 만들기 위해 K개의 축 회전이 발생하는 것이다. 축을 회전시키는 접근은 앙상블 모델의 다양성을 높이기 위한 것이다. 다양성이 높아지는 이유는 축 회전을 통해 새로운 변수들이 추출(feature extraction)되었기 때문이다[3]). 로테이션 포레스트는 사선 방향으로 결정 경계(decision boundary)를 만들 수 있다는 장점이 있지만, 분산 방향이 특정 구조로 나타나는 데이터가 아니라면 성능향상을 기대하기 어렵다는 단점이 있다.

2.2 Algorithm

2.3 Example

로테이션 포레스트는 rotationForest 패키지의 rotationForest() 함수를 이용하여 구현할 수 있는데, 이항 분류(binary classification)만 지원하고 있다.

이를 kernlab 패키지의 spam 데이터를 통해 확인해보자.

library(rotationForest)
## rotationForest 0.1.3
## Change log: rotationForestNews()
data(spam, package = "kernlab")
str(spam)
## 'data.frame':    4601 obs. of  58 variables:
##  $ make             : num  0 0.21 0.06 0 0 0 0 0 0.15 0.06 ...
##  $ address          : num  0.64 0.28 0 0 0 0 0 0 0 0.12 ...
##  $ all              : num  0.64 0.5 0.71 0 0 0 0 0 0.46 0.77 ...
##  $ num3d            : num  0 0 0 0 0 0 0 0 0 0 ...
##  $ our              : num  0.32 0.14 1.23 0.63 0.63 1.85 1.92 1.88 0.61 0.19 ...
##  $ over             : num  0 0.28 0.19 0 0 0 0 0 0 0.32 ...
##  $ remove           : num  0 0.21 0.19 0.31 0.31 0 0 0 0.3 0.38 ...
##  $ internet         : num  0 0.07 0.12 0.63 0.63 1.85 0 1.88 0 0 ...
##  $ order            : num  0 0 0.64 0.31 0.31 0 0 0 0.92 0.06 ...
##  $ mail             : num  0 0.94 0.25 0.63 0.63 0 0.64 0 0.76 0 ...
##  $ receive          : num  0 0.21 0.38 0.31 0.31 0 0.96 0 0.76 0 ...
##  $ will             : num  0.64 0.79 0.45 0.31 0.31 0 1.28 0 0.92 0.64 ...
##  $ people           : num  0 0.65 0.12 0.31 0.31 0 0 0 0 0.25 ...
##  $ report           : num  0 0.21 0 0 0 0 0 0 0 0 ...
##  $ addresses        : num  0 0.14 1.75 0 0 0 0 0 0 0.12 ...
##  $ free             : num  0.32 0.14 0.06 0.31 0.31 0 0.96 0 0 0 ...
##  $ business         : num  0 0.07 0.06 0 0 0 0 0 0 0 ...
##  $ email            : num  1.29 0.28 1.03 0 0 0 0.32 0 0.15 0.12 ...
##  $ you              : num  1.93 3.47 1.36 3.18 3.18 0 3.85 0 1.23 1.67 ...
##  $ credit           : num  0 0 0.32 0 0 0 0 0 3.53 0.06 ...
##  $ your             : num  0.96 1.59 0.51 0.31 0.31 0 0.64 0 2 0.71 ...
##  $ font             : num  0 0 0 0 0 0 0 0 0 0 ...
##  $ num000           : num  0 0.43 1.16 0 0 0 0 0 0 0.19 ...
##  $ money            : num  0 0.43 0.06 0 0 0 0 0 0.15 0 ...
##  $ hp               : num  0 0 0 0 0 0 0 0 0 0 ...
##  $ hpl              : num  0 0 0 0 0 0 0 0 0 0 ...
##  $ george           : num  0 0 0 0 0 0 0 0 0 0 ...
##  $ num650           : num  0 0 0 0 0 0 0 0 0 0 ...
##  $ lab              : num  0 0 0 0 0 0 0 0 0 0 ...
##  $ labs             : num  0 0 0 0 0 0 0 0 0 0 ...
##  $ telnet           : num  0 0 0 0 0 0 0 0 0 0 ...
##  $ num857           : num  0 0 0 0 0 0 0 0 0 0 ...
##  $ data             : num  0 0 0 0 0 0 0 0 0.15 0 ...
##  $ num415           : num  0 0 0 0 0 0 0 0 0 0 ...
##  $ num85            : num  0 0 0 0 0 0 0 0 0 0 ...
##  $ technology       : num  0 0 0 0 0 0 0 0 0 0 ...
##  $ num1999          : num  0 0.07 0 0 0 0 0 0 0 0 ...
##  $ parts            : num  0 0 0 0 0 0 0 0 0 0 ...
##  $ pm               : num  0 0 0 0 0 0 0 0 0 0 ...
##  $ direct           : num  0 0 0.06 0 0 0 0 0 0 0 ...
##  $ cs               : num  0 0 0 0 0 0 0 0 0 0 ...
##  $ meeting          : num  0 0 0 0 0 0 0 0 0 0 ...
##  $ original         : num  0 0 0.12 0 0 0 0 0 0.3 0 ...
##  $ project          : num  0 0 0 0 0 0 0 0 0 0.06 ...
##  $ re               : num  0 0 0.06 0 0 0 0 0 0 0 ...
##  $ edu              : num  0 0 0.06 0 0 0 0 0 0 0 ...
##  $ table            : num  0 0 0 0 0 0 0 0 0 0 ...
##  $ conference       : num  0 0 0 0 0 0 0 0 0 0 ...
##  $ charSemicolon    : num  0 0 0.01 0 0 0 0 0 0 0.04 ...
##  $ charRoundbracket : num  0 0.132 0.143 0.137 0.135 0.223 0.054 0.206 0.271 0.03 ...
##  $ charSquarebracket: num  0 0 0 0 0 0 0 0 0 0 ...
##  $ charExclamation  : num  0.778 0.372 0.276 0.137 0.135 0 0.164 0 0.181 0.244 ...
##  $ charDollar       : num  0 0.18 0.184 0 0 0 0.054 0 0.203 0.081 ...
##  $ charHash         : num  0 0.048 0.01 0 0 0 0 0 0.022 0 ...
##  $ capitalAve       : num  3.76 5.11 9.82 3.54 3.54 ...
##  $ capitalLong      : num  61 101 485 40 40 15 4 11 445 43 ...
##  $ capitalTotal     : num  278 1028 2259 191 191 ...
##  $ type             : Factor w/ 2 levels "nonspam","spam": 2 2 2 2 2 2 2 2 2 2 ...

rotationForest 패키지에서는 target variable을 0 또는 1의 값을 갖는 factor로 제한하고 있기 때문에 이에 맞춰어 데이터를 바꿔줘야 한다. 또한 범주형 예측변수들은 미리 더미코딩을 해서 넣어주어야 한다.

  • y : A factor containing the response vector. Only {0,1} is allowed.
spam$type <- factor(ifelse(spam$type=="spam", 1, 0))
set.seed(1)
train.idx <- sample(1:nrow(spam), nrow(spam)*0.7)
train.X <- spam[train.idx,-58]
train.Y <- spam[train.idx,58]
test.X <- spam[-train.idx,-58]
test.Y <- spam[-train.idx,58]
fit.rot <- rotationForest(x=train.X, y=train.Y, K = round(ncol(spam)/3, 0), L = 10) # default K, L
pred.rot <- predict(fit.rot, newdata = test.X)
head(pred.rot)
## [1] 0.9285774 0.8945287 0.7293693 0.9283171 0.9341454 0.6893023

분류 문제의 성능 평가를 위해 주로 사용하는 AUC를 계산하기 위해 ROCR 패키지를 사용하였다.

먼저 ROC 곡선을 그려보자.

library(ROCR)
## Loading required package: gplots
## 
## Attaching package: 'gplots'
## The following object is masked from 'package:stats':
## 
##     lowess
pred.obj <- prediction(pred.rot, test.Y)
plot(performance(pred.obj, "tpr", "fpr"))
abline(a=0,b=1)

AUC 값은 performance 함수에서 ’auc’를 지정하여 계산된 객체에서 y.values이다.

auc <- performance(pred.obj, "auc")
auc@y.values
## [[1]]
## [1] 0.9614531

Accuracy 값은 ’acc’를 지정하여 구할 있다. cutoff 값에 따라 Accuracy가 어떻게 변하는지 살펴보자.

acc <- performance(pred.obj, "acc")
plot(acc)

  • cutoff 값이 0.4에서 0.7 사이일 때는 Accuracy가 거의 일정한 것을 확인할 수 있다.

로테이션 포레스트에서 지정해줘야 하는 하이퍼 파라미터에는 L과 K가 있는데 이 값을 조정함에 따라 AUC 값이 어떻게 변하는지 알아보자.

L : The number of base classifiers (trees using the rpart package). The default is 10.

K : The number of variable subsets. The default is the value K that results in three features per subset.

먼저 K값은 기본값인 round(ncol(spam)/3, 0)(spam 데이터의 경우 19)으로 고정시키고, L값을 변화시켜보자.

set.seed(1)
folds <- sample(1:5, size = nrow(spam), replace = T) # 5 fold CV
cv.auc <- matrix(0, nrow=6, ncol=5)
L.grid <- c(5,10,15,20,25,30)

for (j in 1:5){
  fold.idx <- folds==j
  for (i in 1:6){
    print(paste0(j, "th folds ", "L=", i))
    cv.rot <- rotationForest(x=spam[-fold.idx,-58], y=spam[-fold.idx,"type"], 
                             K = round(ncol(spam)/3, 0), 
                             L = L.grid[i])
    cv.pred <- predict(cv.rot, newdata=spam[fold.idx,-58])
    cv.auc[i,j] <- unlist(performance(prediction(cv.pred, spam[fold.idx,58]), "auc")@y.values)
  }
}
## [1] "1th folds L=1"
## [1] "1th folds L=2"
## [1] "1th folds L=3"
## [1] "1th folds L=4"
## [1] "1th folds L=5"
## [1] "1th folds L=6"
## [1] "2th folds L=1"
## [1] "2th folds L=2"
## [1] "2th folds L=3"
## [1] "2th folds L=4"
## [1] "2th folds L=5"
## [1] "2th folds L=6"
## [1] "3th folds L=1"
## [1] "3th folds L=2"
## [1] "3th folds L=3"
## [1] "3th folds L=4"
## [1] "3th folds L=5"
## [1] "3th folds L=6"
## [1] "4th folds L=1"
## [1] "4th folds L=2"
## [1] "4th folds L=3"
## [1] "4th folds L=4"
## [1] "4th folds L=5"
## [1] "4th folds L=6"
## [1] "5th folds L=1"
## [1] "5th folds L=2"
## [1] "5th folds L=3"
## [1] "5th folds L=4"
## [1] "5th folds L=5"
## [1] "5th folds L=6"
rowMeans(cv.auc)
## [1] 0.9548381 0.9636836 0.9695326 0.9684165 0.9685175 0.9697847
plot(L.grid, rowMeans(cv.auc), type="l", xlab="L", ylab="AUC")

  • L 값이 너무 작지만 않다면 괜찮아 보인다.

이번에는 L값을 기본값인 10으로 고정시키고, K값을 변화시켜 보자.

set.seed(1)
folds <- sample(1:5, size = nrow(spam), replace = T) # 5 fold CV
cv.auc <- matrix(0, nrow=7, ncol=5)
K.grid <- c(4,8,12,16,20,24,28)

for (j in 1:5){
  fold.idx <- folds==j
  for (i in 1:7){
    print(paste0(j, "th folds ", "K=", i))
    cv.rot <- rotationForest(x=spam[-fold.idx,-58], y=spam[-fold.idx,"type"], 
                             K = K.grid[i], 
                             L = 10)
    cv.pred <- predict(cv.rot, newdata=spam[fold.idx,-58])
    cv.auc[i,j] <- unlist(performance(prediction(cv.pred, spam[fold.idx,58]), "auc")@y.values)
  }
}
## [1] "1th folds K=1"
## [1] "1th folds K=2"
## [1] "1th folds K=3"
## [1] "1th folds K=4"
## [1] "1th folds K=5"
## [1] "1th folds K=6"
## [1] "1th folds K=7"
## [1] "2th folds K=1"
## [1] "2th folds K=2"
## [1] "2th folds K=3"
## [1] "2th folds K=4"
## [1] "2th folds K=5"
## [1] "2th folds K=6"
## [1] "2th folds K=7"
## [1] "3th folds K=1"
## [1] "3th folds K=2"
## [1] "3th folds K=3"
## [1] "3th folds K=4"
## [1] "3th folds K=5"
## [1] "3th folds K=6"
## [1] "3th folds K=7"
## [1] "4th folds K=1"
## [1] "4th folds K=2"
## [1] "4th folds K=3"
## [1] "4th folds K=4"
## [1] "4th folds K=5"
## [1] "4th folds K=6"
## [1] "4th folds K=7"
## [1] "5th folds K=1"
## [1] "5th folds K=2"
## [1] "5th folds K=3"
## [1] "5th folds K=4"
## [1] "5th folds K=5"
## [1] "5th folds K=6"
## [1] "5th folds K=7"
plot(K.grid, rowMeans(cv.auc), type="l", xlab="K", ylab="AUC")

K값이 작을 때 AUC 값이 높은 것을 볼 수 있다. subset 개수인 K를 너무 크게 잡으면 안된다고 볼 수 있다.

3 Decision Jungle

3.1 Intro

앞에서 살펴본 랜덤 포레스트는 모델을 구축하는데 드는 시간이 짧고, 예측 정확도가 높다는 장점이 있다. 하지만 트리를 전개시키는 과정에서 노드 수가 기하급수적으로 증가하기 때문에 메모리가 많이 필요하다는 단점이 있다. 모바일이나 임베디드 프로세서의 경우 메모리 리소스가 제한적이기 때문에 트리의 깊이를 제한 할 수 밖에 없고, 이에 따라 모델의 정확도도 낮아지게 된다. 디시전 정글은 다범주(multi-class) 분류 문제를 풀기 위해 만들어진 DAG(Directed Acyclic Graphs) 모델의 앙상블이다. 마이크로소프트 리서치에서 개발한 알고리즘이고, 마이크로소프트의 머신러닝 플랫폼인 Azure에 적용되어 있다. 전통적인 트리 모델과는 달리 디시전 정글은 parent node들이 모든 child node들과 연결되는 것을 허용한다. 좀 더 정확히 말하면, 두 개 이상의 child node 사이의 병합을 허용하는 것이다. 디시전 정글은 랜덤 포레스트보다 모델 구축에 걸리는 시간은 증가하지만 메모리 사용은 획기적으로 감소하고, 예측 성능도 확보할 수 있다[4].

3.2 Algorithm

디시전 정글의 original paper([4])를 보면 아래의 Figure를 통해 알고리즘을 설명하고 있다.

(b)에서 \(N_p\)는 parent node들의 집합이고, \(N_c\)는 child node들의 집합이다. 여기서 child node의 maximum 값을 지정하여 노드 수가 기하급수적으로 늘어나는 것을 방지한다. \(\theta_i\)는 i번째 parent node가 split에 어떤 변수를 사용했는지, split point가 어디인지를 나타내는 파라미터들이다. \(S_i^L\)은 i번째 parent node의 left child node이고, \(S_i^R\)은 i번째 parent node의 right child node이다.

\(l_i\)와 \(r_i\)는 i번째 parent node에서 어떻게 split이 되었고, i번째 parent node의 left child node와 right child node가 무엇인지가 다 로그로 기록된다는 것을 의미한다. 이 로그 데이터를 바탕으로 최적의 split을 찾는 과정을 반복하는데, 이러한 점 때문에 모델 구축 과정이 랜덤 포레스트보다 더 오래 걸리게 된다.

디시전 정글의 objective function은 다음과 같다:

\[E(\{ \theta \}, \{l_i\}, \{r_i\})=\sum_{j \in N_c}|S_j|H(S_j)\]

여기서 \(|S_j|\)는 j번째 child node의 인스턴스 개수이고, \(H(S_j)\)는 child node의 entropy이다. 이제 ojective function을 minimize 시켜야 하는데 이를 direct로 푸는 것은 쉽지 않아서 근사적인 방법을 사용해야 한다.

3.2.1 Lsearch

Lsearch의 핵심은 feasible solution으로 부터 트리 분할을 시작한다는 것이다. 다시 말해, random한 split point에서 출발한다는 것을 의미한다. 파라미터들을 랜덤하게 할당한 후 split optimization과 branch optimization을 반복적으로 수행하는 것이 Lsearch이다. split optimization 단계에서는 parent node에서 child node로 가는 화살표는 모두 고정한 상태에서 각각의 parent node의 split point를 바꿔가면서 전체 entropy가 가장 낮게 하는 split 조합을 찾는 것을 목적으로 한다. branch optimization 단계에서는 left child node와 right child node를 어떻게 할당하는 것이 최적인지 찾는다. 하나의 parent node의 최적 left/right child node 할당을 찾는 과정에서 나머지 parent 노드들의 화살표는 고정된다. 그런 다음 split optimization과 branch optimization을 더 이상 전체 entropy의 감소가 없을 때까지 반복한다.