본문으로 바로가기



능선회귀분석 with R - 3. Coefficient paths 그리기

category Statistics/Linear models 2017.06.16 10:58
능선회귀분석 with R - 3. Coefficient paths 그리기

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



이전 시간(http://www.issactoast.com/80)에 우리는 실제 데이터를 가지고 주어진 \(\lambda\)에 대한 능선 회귀의 베타에 대한 추정값을 구하는 방법에 대하여 알아보았다. 이번 시간에는 능선회귀 계수들을 \(\lambda\)를 변화시켜가면서 어떻게 달라지는 지에 대하여 알아보도록 하자.

능선 회귀의 Effective degrees of freedom

먼저 알아야 할 개념은 바로 effective degrees of freedom이라는 개념인데 능선회귀의 계수의 자유도에 관한 개념이다. 이전 포스팅에서 능선회귀의 해는 패널티 조절 모수인 \(\lambda\)에 의하여 영향을 받는다는 것을 보인 적이 있다.

\(\lambda\)의 값이 점점 커짐에 따라서 베타의 값들은 점점 0으로 수렴한다. 이런 측면에서 보았을 때 베타가 가질 수 있는 값들의 자유도가 점점 줄어든다고 말 할 수 있을 것이다. 반대로 \(\lambda\)의 값이 0일 경우에는 계수는 p개의 값을 갖게 된다. 이러한 현상을 잘 반영하는 측정 도구로서 제안된 함수가 바로 effective degrees of freedom이다. 정의는 다음과 같다.

\[ df\left(\lambda\right) =tr\left[X\left(X^{T}X+\lambda I\right)^{-1}X^{T}\right] \]

위와 같이 정의된 함수를 행렬의 SVD분해(singular value decomposition)를 이용하여 좀 더 간단하게 만들 수 있다. 이 이론에 따르면 데이터 매트릭스 \(X\)를 다음과 같이 세 개 행렬의 곱으로 나타낼 수 있다.

\[ X=U\Sigma V^{T} \]

위의 식에서 \(U\)\(V\) 행렬은 직교 행렬이며, \(\Sigma\) 경우는 대각행렬 \(\Sigma=diag\left(\sigma_{1},...,\sigma_{p}\right)\)를 의미한다. 이것을 가지고 effective degrees of freedom 함수를 다음과 같이 간단히 표현 할 수 있다.

(이제부터의 전개되는 수식이 이해가 안 되시는분들은 이전 능선회귀(Ridge regression)에 대하여 2 - λ 와 β 의 관계 포스팅에 SVD의 개략적 설명와 함께 수식에 대해 자세하게 설명이 되어 있으니 같이 보시면 이해가 될 것이라 생각합니다.)

\[ \begin{align*} df\left(\lambda\right) & =tr\left[X\left(X^{T}X+\lambda I\right)^{-1}X^{T}\right]\\ & =tr\left[\left(X^{T}X+\lambda I\right)^{-1}X^{T}X\right]\\ & =tr\left[\left(V\Sigma U^{T}U\Sigma V^{T}+\lambda I\right)^{-1}V\Sigma U^{T}U\Sigma V^{T}\right]\\ & =tr\left[\left(V\Sigma^{2}V^{T}+\lambda I\right)^{-1}V\Sigma^{2}V^{T}\right]\\ & =tr\left[\left(V\left(\Sigma^{2}+\lambda I\right)V^{T}\right)^{-1}V\Sigma^{2}V^{T}\right]\\ & =tr\left[V\left(\Sigma^{2}+\lambda I\right)^{-1}V^{T}V\Sigma^{2}V^{T}\right]\\ & =tr\left[\left(\Sigma^{2}+\lambda I\right)^{-1}\Sigma^{2}\right]\\ & =tr\left[diag\left(\frac{\sigma_{1}^{2}}{\sigma_{1}^{2}+\lambda},...,\frac{\sigma_{p}^{2}}{\sigma_{p}^{2}+\lambda}\right)\right]\\ & =\sum_{i=1}^{p}\left(\frac{\sigma_{i}^{2}}{\sigma_{i}^{2}+\lambda}\right) \end{align*} \]

위와 같이 얻게 된 자유도 함수는 앞에서 언급했던 능선회귀 계수와 \(\lambda\)와의 관계를 잘 반영하고 있다는 것을 알 수 있다.

    1. 먼저 람다가 무한대로 갈 경우 \(df(\lambda)\)의 값으로 수렴한다.
    1. 람다의 값이 0일 경우(이 경우 능선회귀는 일반 회귀분석의 경우가 된다.), 위의 자유도 함수의 값은 \(p\)가 된다. 이는 일반 회귀분석에서의 자유도와 같아진다.
    1. 위 자유도 함수는 monotone decreasing하며, 람다와 1 대 1 대응 함수이다.

3의 내용은 다음과 같이 R 코드를 사용하여 쉽게 볼 수 있다.

# generate arbitrary data points
X <- matrix(rnorm(100), ncol = 5)

# apply SVD to get d's
d_vec <- svd(X)$d

# df function def.
df_lambda <- function(lambda, dvec){
    sum(dvec^2 / (dvec^2 + lambda))    
}

# Draw plot
lambda <- seq(0, 100, by = 0.1)
df <- sapply(lambda, df_lambda, dvec = d_vec)
plot(lambda, df, type = "l")

능선회귀 계수 그래프 그리기

이제 능선회귀 coefficient paths를 그릴 준비가 다 되었다. 이전 포스팅에서 사용했던 코드를 불러오도록 하자. (코드 앞에 붙은 주석의 URL은 코드과 관련된 포스팅 주소입니다.)

library(magrittr)

# 데이터 준비 하기 (http://www.issactoast.com/79)
url_Data <- readLines("https://raw.githubusercontent.com/issactoast/Machine-Learning-with-R/master/Data/house.txt") %>% paste(collapse = "\n")

houseData <- read.table(text = url_Data, header = FALSE, sep = ":")

for (i in 1:dim(houseData)[2]){
    houseData[,i] <- gsub(paste("",i), "", houseData[,i])
}

houseData <- houseData %>% unlist() %>% as.numeric() %>% matrix(ncol = 14)

colnames(houseData) <- c("MEDV", "CRIM", "ZN", "INDUS",
                         "CHAS", "NOX", "RM", "AGE",
                         "DIS", "RAD", "TAX", "PTRATIO",
                         "B", "LSTAT")

# 데이터 분리 하기 (http://www.issactoast.com/80)
set.seed(1234)
index <- sample(506, 355)
train_Data <- houseData[index,]
test_Data <- houseData[-index,]

y <- train_Data[,1]
y_bar <- mean(y)
y <- y - y_bar

X <- train_Data[,-1]
x_bar <- colMeans(X)
X <- scale(X)
#X <- sweep(X, 2, x_bar)

# 능선회귀 계수추정 함수 정의(http://www.issactoast.com/80)
get_ridgeSol <- function(lambda, x_data, y_data){
  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))
  names(b_hat) <- colnames(x_data)
  b_hat
}

위의 코드에서 데이터 매트릭스는 X변수에 저장되어 있다. 따라서 앞에서 했던 것과 똑같이 SVD를 다음과 같이 적용하여 d값들을 구할 수 있다.

# apply SVD to get d's
d_vec <- svd(X)$d

# df function def.
df_lambda <- function(lambda, dvec){
    sum(dvec^2 / (dvec^2 + lambda))    
}

이제 그래프를 그리기 위하여, 그래프의 X축에 대응하는 df(X)벡터를 발생시키자. 자유도 함수가 람다와 1대1 대응을 한다는 사실을 이용하여 발생시킨 df값에 대응하는 람다값을 uniroot()를 사용하여 구한다. 아래의 f 함수와 g 함수는 람다값을 한꺼번에 구하기 위하여 사용한 sapply()를 이용하기 위하여 정의하였다.

위의 방법으로 람다를 얻게 된 후, 이와 대응하는 능선회귀의 계수들은 위의 코드 마지막 부분에 정의된 get_ridgeSol() 함수를 이용하여 구한다. 다음의 코드를 살펴보면서 말한 내용을 이해해보도록 하자.

# generate df values for x axis of coefficient path plot
df <- seq(0.5, 13, by= 0.5)

# obtain lambda values corresponding df 
f <- function(x, y){df_lambda(x, d_vec) - y}
g <- function(df){ uniroot(f, c(0, 10^8), tol = 0.0001, y = df)$root}
lambda <- sapply(df, g)

# obtain coefficients w.r.t. lambda
beta_hat <- sapply(lambda, get_ridgeSol, x_data = X, y_data = y)
datafor_path <- cbind(df, t(beta_hat)) %>% as.data.frame()
datafor_path <- rbind(0, datafor_path)

마지막 0을 추가해준 코드는 이론상 람다가 무한대일 경우 df값은 0에 대응하고, 모든 계수는 0으로 수렴하는 것을 반영해주는 코드이다. 이것들을 이용하여 이번 포스팅의 마지막 단계인 능선회귀의 Coefficient paths 그래프를 그려보자.

library(reshape2)
library(ggplot2)

# melt data for applying ggplot
melted <- melt(datafor_path, id.var = "df")
p <- ggplot(data = melted, aes(x = df,
                               y = value,
                               group = variable,
                               color = variable)) + geom_line()

# text postion & label setting
t_position <- rep(14, 13)
t_position[(1:6)*2] <- 16
varNames <- sort(datafor_path[27, -1]) %>% colnames()

# Draw coeffiecnts paths of ridge regression
p + annotate("text", x = t_position, y = as.numeric(sort(datafor_path[27, -1])), label = varNames)

자유도, 결국, 람다의 값에 따른 각 계수들 변화를 한눈에 알아볼 수 있게 되었다. 이번 포스팅을 마치면서 몇가지 짚고 넘어가야 할 것이 있는데, 예전 포스팅에서 필자는 데이터를 평균만 뺀 centered 데이터를 가지고 계수를 구하였다. 하지만 이번 능선회귀 계수 그래프를 그리면서 scale()함수를 사용하여 정규화된 데이터를 사용하였다. 그 이유는 특정 변수(NOX)의 계수가 너무 높게 나와 변수들 간의 변화가 한눈에 들어오지 않게 되었기 때문이다.

사실, 필자 자신도 어떤 경우 센터링만을, 어떤 경우에는 스케일링까지 해야 하는가에 대한 질문에 답을 하기가 쉽지 않다. 나중에 기회가 된다면, 데이터 스케일링에 대한 포스팅을 한번 쓰고 싶은 마음이다. 이 토픽에 대한 좋은 토론은 이곳을 참고하자.

Reference

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

[2] PennState STAT 897D(https://onlinecourses.science.psu.edu/stat857/node/155 )


SHARE TO



티스토리 툴바