머신러닝 예제

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

이하에서는 회귀(regression) 문제에서 ridge와 lasso 회귀를 실습한다. glmnet 패키지를 사용하는 방법과 h2o 패키지(설치)를 사용하는 방법을 모두 설명한다. 두 패키지 중 glmnet이 더 간편하고 안정적인 것으로 보임에도 h2o 패키지를 소개하는 것은, h2o가 deep learning을 제대로 구현하고 있어, 독자들이 일단 여기에 노출되기를 바라기 때문이다.

Ridge regression

그림이 특이하다. λ가 증가함에 따라 보통은 변수들의 계수가 수축되는데 지금은 어떤 변수의 계수가 수축되는 것이 아니라 확장되었다가 수축된다. 매우 특이하며, 필자는 이런 경우를 별로(거의) 보지 못했다.

다음으로 10-fold CV 방법으로 최적 λ를 정한다. 먼저 cv.glmnet 명령을 이용하여 CV를 하고 결과를 그림으로 그려 보자.

CV를 최소화시키는 λ 값은 아래에서 보듯이 33.27139이다. Test set (z15)에 대하여 CV를 이용한 최적 ridge를 적용한 예측의 성과는 앞의 변수선택 방법보다 나쁘며 단순한 random walk와 거의 비슷하다.

보통은 ridge의 결과가 좋은데 여기서 이렇게 안 좋은 것은 특이하다. 이는 lambda가 제대로 선택되지 않았기 때문으로 보인다. 특히 위 그림에서 CV를 최소화시키는 λ (왼쪽 점선)가 후보군 내에 거의 최솟값인 것 같다. 이는 최적 λ는 그보다 더 작을 수도 있음을 의미한다. 이하에서 이 점을 개선해 보자.

Ridge에서 λ 선택 개선

위 ridge의 성능이 좋지 않다. 그 이유가 λ 선택이 좋지 않기 때문이라는 의심이 든다. glmnet 명령이 자동으로 선택한 lambda 값 후보들 중 best.lambda가 어디쯤인지 살펴보기 위해 후보값들 중 가장 작은 6개를 보았더니 다음과 같았다.

CV에 의한 최적 λ값(cv.ridge$lambda.min) 33.27139는 이들 λ 후보 중 가장 작다. 더 작은 λ 값에서는 예측 성능이 더 좋을 수도 있다. 그러니 glmnet으로 하여금 더 작은 λ 값들도 후보군에 포함시키도록 해야 할 것 같다. 이를 위해 lambda.min.ratio를 기본값(0.0001)보다 작은 값(예를 들어 1e-6)으로 설정하고, 이와 더불어 λ 값들을 더 촘촘히 하기 위해 nlambda 옵션(기본값은 100)을 증가시키자. 아래에서는 nlambda를 200으로 설정하였다.

선택된 λ 값은 전체 200개 중 3번째로 작은 값이어서 여전히 불만스러우나, 어쨌든 가장 작은 값은 아니니 됐다. 이 lambda 값을 이용한 test set 예측의 성능은 다음과 같다.

이렇게 조정해 주었더니 ridge의 결과가 훨씬 나아졌다. 계수들과 표준화 계수들은 다음과 같다.

bhat.ridge <- coef(cv.ridge2, s = 'lambda.min')
options(scipen=1)  # fixed digits (not exponential notation)
cbind(bhat.ridge, bhat.ridge * c(1, apply(X, 2, sd)))
# 20 x 2 sparse Matrix of class "dgCMatrix"
#                         1           1
# (Intercept) -282.50349431 -282.503494
# grdp          -1.17756612   -9.311269
# regpop        86.84145413   18.055504
# popgrowth     -2.68510803   -5.301252
# eq5d          -2.05815174   -2.097936
# deaths        -0.01410179  -10.733191
# drink         -1.39124153   -6.211952
# hdrink         2.22434327    7.894561
# smoke          2.34114522    6.159724
# aged          10.99928869   85.067698
# divorce      -14.02562562   -5.159043
# medrate        0.80357585    7.268658
# gcomp         -4.54592170  -23.546566
# vehipc        76.94518991    6.547744
# accpv         -1.34752466   -3.625625
# dumppc        21.42175764    8.453333
# stratio       -0.62597700   -2.995986
# deathrate      0.73503960  237.090334
# pctmale       18.81388329   23.451831
# accpc         -1.20022043   -1.563887

위에서 첫 번째 칼럼은 계수, 두 번째 칼럼은 표준화 계수이다(상수항은 표준화하지 않음). 단, glmnet은 표준화를 위해 표본분산을 구할 때 n − 1로 나누는 것이 아니라 n으로 나누는 반면, R의 sd 명령은 n − 1로 나눈 값에 제곱근을 취하므로, 위의 두 번째 열이 정확한 glmnet의 표준화 계수인 것은 아니나, 크기 비교만 할 것이므로 문제 없다. 계수 자체(변수의 측정단위에 따라 규모가 달라짐)가 아니라 표준화 계수값의 크기가 클수록 중요한 변수이다. deathrate와 aged가 가장 중요한 두 변수이다.

표준화 계수들의 상대적 크기(변수 중요도, 가장 중요한 변수의 중요도를 1.0으로 표준화)를 다음과 같이 그림으로 표현할 수 있다.

참고: glmnet의 ridge 계산 방법에 대한 주석¹

수학에 의하면 ridge 기울기 추정값은 β̃ = (n⁻¹X̃′X̃+λI)⁻¹n⁻¹X̃′y 수식으로 구할 수 있다(X̃는 x_j에서 그 표본평균을 차감한 값들의 행렬). 다음으로 절편은 α̃ = ȳ − β̃₁x̄₁ − ⋯ − β̃_px̄_p 수식에 따라 구한다. 그런데 특별히 디자인한 경우가 아니면, 같은 자료와 같은 λ를 사용한다 할지라도 glmnet이 구해 주는 ridge 추정값들은 위 수식으로 구한 ridge 추정값들과 다르다. 이는 glmnet이 변수 표준화와 관련하여 특별한 일을 하기 때문이다.

우선, glmnet은 ridge 계산 전에 목표변수(y)를 표준화하고 계산 종료 후 되돌린다. 그러므로 glmnet에서 사용하는 ’λ’는 수식에 나오는 벌점이 아니다. glmnet이 사용한 λ에 대하여 위 행렬 연산을 통하여 직접 ridge 추정값을 구하려면, σ̂_y을 목표변수의 표본표준편차(표본분산 계산 시 n − 1이 아니라 n으로 나눔)라 할 때, 벌점을 λ^* = λ/σ̂_y로 바꾸어야 한다. 이렇게 벌점을 바꾸고 위 행렬 연산을 하면 glmnet에서 x_j를 표준화하지 않은(standardize = FALSE 옵션 사용) 계수 추정값을 구할 수 있다.

그런데 glmnet은 디폴트로 특성변수들(x_j)도 표준화한다. 즉, σ̂_j를 x_j의 표본표준편차(표본분산 계산 시 n − 1이 아니라 n으로 나눔)라 하면 표준화는 모형의 β_jx_j를 (β_jσ̂_j)(x_j/σ̂_j)로 변환한 것과 마찬가지여서, 변수가 x_j/σ̂_j가 되고 계수는 β̂_jσ̂_j가 되는 효과를 갖는다. 표준화의 영향을 되돌리리면 계수 추정값을 다시 σ̂_j로 나누어야 한다.

이상을 정리하면, glmnet에 standardize = TRUE 옵션(디폴트)을 주고 계산하는 ridge 추정값은 다음과 같이 구할 수 있다. ① x̃_ij^* = (x_ij−x̄_j)/σ̂_j로 표준화한다. 여기서 x̄_j는 x_ij의 표본평균, σ̂_j는 x_ij의 표본평균편차(표본분산 계산 시 n으로 나눔). 표본평균을 빼는 것은 절편을 처리하기 위함이며, σ̂_j로 나누는 것은 변수 표준화. ② β̂^* = (n⁻¹X̃^*′X̃^*+λ^*I)⁻¹n⁻¹X̃^*′y를 구한다(X̃^*는 x̃_ij^*의 행렬, λ^* = λ/σ̂_y). 이것은 ‘베타계수’ β_jσ̂_j들의 추정값 벡터이다. ③ 표준화를 되돌리기 위해 β̂_j = β̂_j^*/σ̂_j를 계산하면 이것이 최종적인 기울기 계수의 ridge 추정값. ④ 절편 추정값은 앞과 동일하게 α̂ = ȳ − β̂₁x̄₁ − ⋯ − β̂_px̄_p로 계산.

참고로, 위 ①∼③의 β̂은 D = diag(σ̂₁²,…,σ̂_p²)라 할 때 (n⁻¹X̃′X̃+λ^*D)⁻¹n⁻¹X̃′y으로 계산한 것과 동일하다.

이상의 내용은 R을 이용하여 수치적으로 확인하였다. 확인한 바는 여기에 있다.

Lasso

Lasso를 위해서는 glmnet과 cv.glmnet에서 alpha를 0에서 1로 바꾸면 된다.

이로부터 최적 λ를 구하고 이에 대응하는 lasso 추정 결과를 test set (z15)에 적용하여 사망률을 예측한 성과는 다음과 같다.

테스트셋에서 lasso의 RMSE는 random walk보다 더 낫고 Best Subset Selection보다도 더 낫다.

예측오차 제곱합이 최저인 λ(lambda.min)보다는 “one standard error” (1se)를 사용하는 경우도 많다. 이를 위해서는 'lambda.min' 대신에 'lambda.1se'를 사용하면 된다.

Ridge와 달리 lasso는 변수선택도 해 준다. Lasso에 의하여 선택된 변수들과 그 계수는 coef 명령으로써 확인할 수 있다. 또한, 변수들의 중요도를 표준화 계수(beta coefficient)의 크기, 즉 X_j를 표준화한 상태에서의 계수 추정값(β_j에 X_j의 표본표준편차를 곱함)의 크기로써 표현할 수 있다.

표준화 계수 추정값의 크기에 의하면, 직전연도 사망률(deathrate)이 가장 중요하고 그 다음으로 고령인구비율(aged)이 중요한 예측변수이다. 여기에 흡연인구 비율(smoke), 고위험 음주인구 비율(hdrink), 인구증가율(popgrowth), 1인당 자동차 수(vehipc)가 사망률 예측에 미세한 도움을 준다.

이 예에서, 2015년 시험용 자료에서 그 다음 해 deathrate (즉, ynext)의 예측에는 lasso가 ridge보다 더 나은 예측성과를 보였지만, 다른 목표변수나 다른 자료의 경우에는 결과가 달라질 수 있다. Lasso가 변수선택을 해 준다는 사실과 예측을 얼마나 잘하느냐는 별개의 문제이다. 참 모형이 sparse (많은 변수들의 계수가 0)할 경우에는 특정 lasso가 마치 어느 변수의 계수가 0인지 정확히 아는 것과 같다는 성질(’오라클’이라 함)이 있다는 결과들도 있으나, 참 모형이 sparse한지 알지 못하므로 실질적으로는 큰 의미가 없다. 그냥, ridge가 더 나은 경우도 있고 lasso가 더 나은 경우도 있다고 이해하는 편이 간편하다. 특정 문제에서 목표변수가 관측되지 않은 상태에서 ridge가 더 나은지 lasso가 더 나은지 알 수는 없다. 하나의 문제에서 lasso가 더 나았다고 하여 다른 문제에서도 그러라는 법은 없다. 그러지 말라는 법도 없다.

Elastic Net

α = 0은 ridge, α = 1은 lasso이다. Elastic net은 일반적인 α를 사용하는 방법을 통칭한다. α = 0.5를 사용하면 결과는 다음과 같다.

예측 오차가 아주 작은데(RMSE = 46.06922), 이는 우연히 그런 것이다. α = 0.5는 체계적인 방법으로 최적화하여 정한 것이 아니라, 한번 시도해 본 것뿐이다. 그 대신 α를 CV 방법으로 정할 수도 있다.

Grid search를 이용한 α 매개변수의 선택

많이 사용되지는 않지만, grid search를 통하여 최적의 α 매개변수를 선택하는 것도 가능하다. 각각의 α 값에서 CV를 통하여 λ를 찾고 이 λ에 해당하는 CV error를 기록하여 비교하면 된다. 이하에서는 α = 0, 0.1, …, 0.9, 1에 대하여 비교하고자 한다.

pset <- data.frame(alpha = seq(0,1,by=.1), lambda.min = NA, cv.err = NA)

set.seed(1)
for (i in 1:nrow(pset)) {
  set.seed(1)  # important for same fold ID for all alpha
  cv1 <- cv.glmnet(X,Y,alpha = pset$alpha[i], lambda.min.ratio = 1e-6, nlambda = 200)
  pset$lambda.min[i] <- cv1$lambda.min
  err <- cv1$cvm[cv1$index['min',1]]
  pset$cv.err[i] <- err
}
pset
#    alpha lambda.min   cv.err
# 1    0.0  0.8204098 2714.922
# 2    0.1  0.8016422 2708.978
# 3    0.2  1.7223300 2689.722
# 4    0.3  1.4140911 2682.232
# 5    0.4  5.6125048 2662.508
# 6    0.5  5.5296669 2642.683
# 7    0.6  4.6080558 2630.376
# 8    0.7  4.2337157 2622.856
# 9    0.8  4.5622807 2617.291
# 10   0.9  4.3469058 2613.710
# 11   1.0  3.9122153 2615.097
which.min(pset$cv.err)
# [1] 10

위 결과에 의하면 α = 0.9일 때 CV error (MSE)가 가장 작고, 이때 최적 λ는 4.3469058이다. (앞에서 최적 α에 해당하는 cv.glmnet 결과를 저장해 놓을 수도 있었으나 시간이 얼마 걸리지 않으므로) α = 0.9로 다시 한번 elastic net을 실행하면 결과는 다음과 같다.

H2O 패키지를 이용한 lasso

H2O.ai (설치는 여기)는 glmnet 패키지보다 훨씬 무겁지만, 방대한 데이터셋을 처리할 수 있고, 더 다양한 옵션을 제공하며, 특히 나중에 살펴볼 인공신경망 등을 포함하므로 한번 연습해 보기로 하자. Lasso (alpha = 1)를 고려한다. 다른 elastic net을 이용하려면 alpha를 바꾸면 된다. Grid search로써 alpha를 정할 수도 있으나, 여기서는 1로 설정한다.

H2O는 데이터셋도 나름의 방식으로 처리하여 사용한다. 이를 위해 as.h2o 명령(약간 시간이 걸림)으로써 데이터셋을 먼저 처리하여 저장해 두는 것이 간편하다. 아래에서는 z14h <- as.h2o(z14) 등처럼 해 두고 z14h 등을 사용할 것이다.

## H2O package
library(h2o) # install.packages("h2o") if necessary
h2o.init()
yvar <- 'ynext'     # name of target variable
xvar <- setdiff(names(z14), yvar) # character vector
z14h <- as.h2o(z14) # h2o's data format
z15h <- as.h2o(z15) # h2o's data format
## lasso with cross validation
glm <- h2o.glm(xvar, yvar, z14h, alpha=1, nfolds=10, lambda_search = TRUE, early_stopping = FALSE, seed=1)
glm@model$lambda_best
# [1] 5.285855
glm@model$coefficients_table
# Coefficients: glm coefficients
#        names coefficients standardized_coefficients
# 1  Intercept    -8.857370                769.412108
# 2       grdp    -0.310505                 -2.455230
# 3     regpop     0.000000                  0.000000
# 4  popgrowth     0.000000                  0.000000
# 5       eq5d     0.000000                  0.000000
# 6     deaths     0.000000                  0.000000
# 7      drink     0.000000                  0.000000
# 8     hdrink     0.279261                  0.991143
# 9      smoke     0.716354                  1.884778
# 10      aged     9.214701                 71.265826
# 11   divorce     0.000000                  0.000000
# 12   medrate     0.000000                  0.000000
# 13     gcomp     0.000000                  0.000000
# 14    vehipc     1.009409                  0.085897
# 15     accpv     0.000000                  0.000000
# 16    dumppc     0.000000                  0.000000
# 17   stratio     0.000000                  0.000000
# 18 deathrate     0.794411                256.240748
# 19   pctmale     0.000000                  0.000000
# 20     accpc     0.000000                  0.000000
## h2o.std_coef_plot(glm)  # xlim bug?
h2o.varimp_plot(glm)
RMSE(z15$ynext, as.vector(h2o.predict(glm, z15h))) # RMSE() defined in index11.php
# [1] 47.86051
h2o.shutdown(prompt = FALSE)

참고로, glmnet과 h2o의 변수 표준화 방식에 약간 차이가 있다. glmnet은 표본표준편차를 구할 때 분모를 n − 1이 아니라 n으로 나누는 반면 h2o에서는 n − 1로 나눈다. 또한, glmnet의 λ를 y의 표본표준편차로 나누어야 h2o의 λ와 유사한 값이 되는 것으로 보인다. lambda_min_ratio (디폴트는 n > p이면 0.0001), nlambdas, solver 등을 잘 지정할 필요가 있어 보이지만, 특이하게도 일반적인 elastic net에서 glmnet처럼 안정적인 결과를 얻기는 매우 어려웠다. 또, glmnet 패키지는 grid search에 사용한 λ 값들을 리턴해 주지만 h2o.glm은 최적 λ만 리턴하는 것 같다(필자가 잘 모르는 것일 수도 있다). 데이터셋이 웬만한 크기여서 glmnet에서도 잘 돌아가기만 한다면 Elastic net (ridge와 lasso 포함)의 경우 h2o보다는 glmnet 패키지를 사용하는 것을 추천한다.

이 내용은 고려대 경제학과 박상수 교수와 함께 알아낸 것임↩︎