본문으로 바로가기



능선회귀분석 with R - 5. Cross Validation을 통한 람다값 정하기

category Statistics/Linear models 2017.06.22 22:52
능선회귀분석 with R - 5. Cross Validation을 통한 람다값 정하기

(이 포스팅은 The Elements of Satistical Learning의 흐름을 따라가고 있음을 밝힙니다. pdf파일은 공개되어 있으니, 다운받아 참고하시면 좋을것 같습니다.)

이전 포스팅에서 우리는 실제 데이터를 가지고, 람다(\(\lambda\))의 변화에 따른 능선회귀의 계수를 나타내주는 coefficient path를 그려보았다. 그렇다면 과연 어떤 람다값을 사용해야할까? 오늘은 람다를 선택하는 방법인 Cross-validation(교차타당성) 기법에 대하여 알아보도록 하자.

Cross validation

최적의 람다를 찾는 방법에 대한 이론은 필자가 아는 한 아직 정립된 방법이 없는 것으로 알고 있다. 단, 휴리스틱(Heuristic이라 쓰고, 야매라고 읽는다)한 방법으로 많은 사람들이 cross-validation이라 부르는 방법을 사용하고 있다.

Cross-validation은 데이터를 분할시켜서 주어진 모델의 성능을 검사하는 방법이라 생각 할 수 있다. 우리는 앞선 포스팅에서 데이터셋을 testtrain 두 부분으로 나누었다. 이는 train 데이터셋을 이용하여 모수를 추정하고, test 데이서셋을 이용하여 모델을 검증하기 위함이었다. 여기서 주의해야 할 점이 있는데, cross-validation 기법에서 데이터를 분할한다는 것은 train 데이터셋을 분할한다는 것이지, test 셋을 건드린다는 것은 아니다.

즉, 데이터를 n등분하여 교차타당성 검증(n-fold cross-validation)을 수행하였다는 의미는 아래의 그림에서와 같이 train data set을 n개의 하위 데이터셋으로 나눈다는 의미이다.

Cross-validation의 핵심 아이디어는 이렇게 나눈 하위 데이터 셋 하나하나를 각각의 독립된 모의 test set으로 가정을 한다는 것이다(n개의 모의고사를 만든거나 마찬가지이다). 그리고 각 데이터 셋들을 사용하여 우리가 후보로 생각하는 모델들의 성능을 검증한 후, 평균적으로 가장 뛰어난 모델을 고르자는 것이다(평균 모의고사 성적이 가장 뛰어난 아이를 뽑아 대회(test set)에 보내는 것과 같은 이치).

4-fold Cross-validation with Boston housing data

이전 포스팅에서 계속해서 사용하고 있는 보스턴 부동산 데이터를 가지고 4-fold cross-validation을 사용하여 람다를 결정해보기로 하자.

1. 데이터셋 분할하기

일단 사용할 데이터를 로딩하자. 데이터 분리 하기 포스팅에서 했던 데이터 클리닝이 적용된 데이터를 필자의 github에 올려놓았다.

library(readr)
library(magrittr)

houseData <- read_csv("https://raw.githubusercontent.com/issactoast/Machine-Learning-with-R/master/Data/houseData.csv", col_names = TRUE) %>% as.matrix.data.frame()

이론상으로 하위 데이터셋을 나눌때는 램던하게 섞어서 나눠야한다. 하지만 실무에서는 그냥 안섞고 한다고 들었다. 이유는 사실 우리가 traintest를 나눌 때 이미 랜덤하게 선택해서 나눠줬기 때문에, 보통 교차타당성 검증을 할 데이터셋이 이미 랜덤하게 섞여있기 때문이다.

4-fold cross-validation을 해야하므로 편의상 traintest를 400개와 106개로 나누는게 편할 것이다. 따라서 나누어진 train 데이터셋의 1에서 100까지의 인텍스는 그룹 1, 101에서 200까지는 그룹 2, 이런 식으로 4개의 그룹이 있다고 생각하자.

set.seed(1234)
index <- sample(506, 400)
train_Data <- houseData[ index,]
test_Data  <- houseData[-index,]

# data split
train_y <- train_Data[, 1]
train_X <- train_Data[,-1]

test_y <- test_Data[, 1]
test_X <- test_Data[,-1]

2. get_ridgeSol() 함수 정의하기

지난 시간에 정의했던 get_ridgeSol() 함수를 살짝 변형시켜서 train_ytrain_X 직접 넣을 수 있도록 변형시키자. 데이터 센터링은 함수 안에서 실행된다.

# 능선회귀 계수추정 함수 정의(http://www.issactoast.com/80)
get_ridgeSol <- function(x_data, y_data, lambda = 0){
    # beta_0 estimation
    b_0 <- mean(y_data)
    y_data <- y_data - b_0
    
    # x_data centering
    x_data <- scale(x_data, center = T, scale = F)
    scale_info <- attributes(x_data)
    
    # solve the eqaution
    a <- t(x_data) %*% x_data + lambda * diag(1, ncol(x_data))
    b <- t(x_data) %*% y_data
    b_hat <- as.numeric(solve(a,b))
    
    # back to the model with intercept
    b_0 <- b_0 - (b_hat %*% scale_info$`scaled:center`)
    b_hat <- c(b_0, b_hat)
    
    # give them colnames
    names(b_hat) <- c("intercept", colnames(x_data))
    
    b_hat    
}

3. 함수 이론값 확인

위에서 정의한 함수의 이론값을 확인해보자. 이론적으로 람다값이 충분히 클 경우 b_0를 제외한 모든 계수들의 값은 0에 근사해야하며, 람다값이 0일 경우는 회귀분석의 계수와 일치해야한다. 다음 코드의 결과에서 우리가 정의한 함수가 잘 정의되었음을 알 수 있다.

# everything goes to zero except intercept
get_ridgeSol(train_X, train_y, lambda = 10^7)
##     intercept          CRIM            ZN         INDUS          CHAS 
##  2.568567e+01 -6.234720e-04  2.445124e-03 -7.346641e-04  1.570695e-05 
##           NOX            RM           AGE           DIS           RAD 
## -1.050903e-05  1.508362e-04 -2.480162e-03  8.124107e-05 -4.378959e-04 
##           TAX       PTRATIO             B         LSTAT 
## -1.217770e-02 -3.061984e-04  5.380606e-03 -1.522511e-03
# lambda = 0, OLS coefficients
my_beta <- get_ridgeSol(train_X, train_y) 

# OLS vis lm() function
beta_lm <- lm(as.formula(paste("MEDV ~ ", paste(colnames(train_X), collapse = "+"))), data = as.data.frame(train_Data))
beta_lm$coefficients
##   (Intercept)          CRIM            ZN         INDUS          CHAS 
##  45.471476946  -0.083686751   0.047491174   0.039260011   2.747514613 
##           NOX            RM           AGE           DIS           RAD 
## -19.703295137   2.946320840  -0.007291839  -1.577448521   0.324745357 
##           TAX       PTRATIO             B         LSTAT 
##  -0.013242670  -1.001927150   0.008531269  -0.564716374

4. fitting evaluation

(image from: http://sebastianraschka.com/Articles/2014_intro_supervised_learning.html)

데이터를 4개로 분할하였으므로, 총 4번의 테스트를 실행할 것이다. 편의상 이렇게 나눈 데이터를 서브데이터 1~4 라고 부르도록 하자. (즉, train_Data[1:100,]은 서브데이터 1이 되는 것이다.) 이들 각각의 서브데이터들은 한번씩 test 데이터 셋의 역할을 하게 된다.

다시말해, 1번째 테스트에서는 train_Data[101:400,]의 데이터(서브데이터 2~4)를 사용하여 계수를 추정하고, train_Data[1:100,]데이터(서브데이터 1)를 사용하여 적합도를 측정한다. 비슷하게 2번째 테스트의 테스트 데이터셋은 서브데이터 2가 되고, 나머지 서브데이터 1,3,4를 사용하여 계수를 추정하게된다.

이때의 적합도를 측정하는 기준은 보통 RMSE(Root-mean-square error)를 사용하는데, 실제값과 예측값의 차이(즉, 에러)값들의 제곱의 평균의 제곱근 값이다.

\[ \operatorname{RMSE}=\sqrt{\frac{\sum_{t=1}^n (\hat y_t - y_t)^2}{n}}. \]

이것을 바탕으로 i = 1 ~ 4 까지의 테스트를 실행 할 예정이다. 이에 앞서 특정 케이스에 대한 코드를 짜보는 것이 좋은데, i = 2에 대한 케이스를 기준으로, \(\lambda\)c(0:19)^2에 대하여 계수를 추정하고, 적합도인 RMSE를 구하는 것을 먼저 해보도록 하자.

# lambda define
lambda <- c(0:19)^2

# special case when i = 2
i <- 2

# train set for cross-validation
cross_trainX <- train_X[-((i - 1) * (100) + (1:100)), ]
cross_trainy <- train_y[-((i - 1) * (100) + (1:100))]

# test set for cross-validation
cross_testX  <- train_X[(i - 1) * (100) + (1:100), ]
cross_testy  <- train_y[((i - 1) * (100) + (1:100))]

# fitting for each lambda
f <- function(x){
     get_ridgeSol(cross_trainX, cross_trainy, lambda = x)
}
coef_ridge <- sapply(lambda, f)
colnames(coef_ridge) <-  paste0("l=", lambda)

# calculating RMSE for each lambda
RMSE_ridge <- (cross_testy - cbind(1, cross_testX) %*% coef_ridge)^2 %>% colMeans() %>% sqrt()

RMSE_ridge
##      l=0      l=1      l=4      l=9     l=16     l=25     l=36     l=49 
## 5.074906 5.020741 4.987244 4.931293 4.868646 4.811913 4.766620 4.734109 
##     l=64     l=81    l=100    l=121    l=144    l=169    l=196    l=225 
## 4.713561 4.703263 4.701334 4.706049 4.715966 4.729922 4.746998 4.766471 
##    l=256    l=289    l=324    l=361 
## 4.787770 4.810440 4.834116 4.858505

이렇게 해서 i = 2 경우, \(\lambda\)의 값이 100근처일 때, RMSE값이 가장 작다는 것을 알 수 있다. 이제 나머지 케이스를 다 고려한 후 평균적으로 가장 작은 RMSE 값을 만들어내는 람다값을 찾아보도록 하자.

cross_valid <- function(lambda){
    # prepare result stoage
    RMSE_ridge <- matrix(0, nrow = 4, ncol = length(lambda))
    colnames(RMSE_ridge) <-  paste0("l=", lambda)

    for (i in 1:4){
        # train set for cross-validation
        cross_trainX <- train_X[-((i - 1) * (100) + (1:100)), ]
        cross_trainy <- train_y[-((i - 1) * (100) + (1:100))]
    
        # test set for cross-validation
        cross_testX  <- train_X[(i - 1) * (100) + (1:100), ]
        cross_testy  <- train_y[((i - 1) * (100) + (1:100))]
    
        # fitting for each lambda
        f <- function(x){
            get_ridgeSol(cross_trainX, cross_trainy, lambda = x)
        }
        coef_ridge <- sapply(lambda, f)
        colnames(coef_ridge) <-  paste0("l=", lambda)
    
        # calculating RMSE for each lambda
        RMSE_ridge[i, ] <- (cross_testy - cbind(1, cross_testX) %*%
                            coef_ridge)^2 %>% colMeans() %>% sqrt()
    }
    RMSE_ridge %>% colMeans()
}

앞에서 정의했던 람다 벡터에 대하여 평균 RMSE값을 그려보도록 하자.

l <- c(0:19)^2
cross_valid(l) %>% plot()

위의 그래프를 보면 초기 다섯개의 람다값에서 값이 발생하는 것으로 보인다. 따라서 좀 더 자세한 람다를 찾기위해 0에서부터 10까지의 값에 대한 평균 RMSE값을 그려보면 다음과 같다.

l <- seq(0, 10, by = 0.05)
mean_RMSE <- cross_valid(l)
plot(l, mean_RMSE)

mean_RMSE %>% which.min()
## l=0.4 
##     9

우리가 그렇게 찾던 람다의 값은 0.4. 따라서 최종 능선회귀의 계수와 RMSE는 다음과 같다.

beta_ridge <- get_ridgeSol(train_X, train_y, lambda = 0.4)
RMSE_ridge <- (test_y - cbind(1, test_X) %*% beta_ridge)^2 %>%
               colMeans() %>% sqrt()

beta_ridge
##     intercept          CRIM            ZN         INDUS          CHAS 
##  41.776845084  -0.078906943   0.048604977   0.016824518   2.637410244 
##           NOX            RM           AGE           DIS           RAD 
## -14.776447588   3.004605444  -0.011482983  -1.507247267   0.313439484 
##           TAX       PTRATIO             B         LSTAT 
##  -0.013714527  -0.946006018   0.008920617  -0.567528627
RMSE_ridge
## [1] 4.307237

이 결과값을 앞에서 구했던 회귀분석의 RMSE와 비교해보자.

beta_lm$coefficients
##   (Intercept)          CRIM            ZN         INDUS          CHAS 
##  45.471476946  -0.083686751   0.047491174   0.039260011   2.747514613 
##           NOX            RM           AGE           DIS           RAD 
## -19.703295137   2.946320840  -0.007291839  -1.577448521   0.324745357 
##           TAX       PTRATIO             B         LSTAT 
##  -0.013242670  -1.001927150   0.008531269  -0.564716374
RMSE_OLS <- (test_y - cbind(1, test_X) %*% beta_lm$coefficients)^2 %>%
               colMeans() %>% sqrt()
RMSE_OLS
## [1] 4.330871

위의 결과에서 우리는 능선회귀의 RMSE의 값이 회귀분석의 RMSE의 값보다 더 작은 것을 확인할 수 있다. 즉, 능선회귀식을 이용한 예측값이 발생시키는 에러값이 회귀분석을 사용하여 예측했을때의 에러보다 평균적으로 더 작다는 것이다.

이제까지 Cross-validation기법을 통한 능선회귀의 \(\lambda\)값을 구하는 방법에 대하여 알아보았다.

Reference

[1] Friedman, Jerome, Trevor Hastie, and Robert Tibshirani. The elements of statistical learning. Vol. 1. Springer, Berlin: Springer series in statistics, 2001.

[2] Bishop, C. “Pattern Recognition and Machine Learning (Information Science and Statistics), 1st edn. 2006. corr. 2nd printing edn.” Springer, New York (2007).


SHARE TO



티스토리 툴바