데이터 셋 분할하기

능선회귀 모형을 앞에서 얻는 houseData 데이터 셋에 적용하기 위하여 데이터 셋을 분할해야 하는데, 보통 training과 test 용도로 나눈다. trainging 셋을 이용하여 능선회귀의 패러미터를 추정하고, test 셋은 추정한 모델의 성능을 판단하는 용도로 쓰인다. 보통 주어진 데이터의 70%를 추정에 사용하고, 나머지 30%를 평가하는데에 사용하게 된다.

houseData의 경우 총 506개의 관측치가 있으므로 대략 70% 정도인 355개를 무작위로 뽑아 train_Data에 나머지 관측값들은 test_Data로 설정하도록 하자.

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

## [1] 355  14

## [1] 151  14

능선회귀의 패러미터 추정하기

예전 포스팅에서 학습했던 능선회귀분석에 대하여 잠깐 언급해보면, 능선회귀에는 두 개의 패러미터, \(\lambda\)와 솔루션 \(\beta\)가 있었다. 능선 회귀의 솔루션 \(\beta\)는 아래의 loss함수를 최소로 만드는 값이었으며, 여기서 \(\lambda\)는 패널티를 얼마나 부여할 것 인지를 조정하는 패러미터로서 작동하였다.

\[ L\left(\beta\right)=\frac{1}{2}\left(y-X\beta\right)^{T}\left(y-X\beta\right)+\frac{\lambda}{2}\left|\left|\beta\right|\right|_{2}^{2} \]

그리고 능선회귀(Ridge regression)의 손실함수를 최소로 만드는 해, \(\hat{\beta}_{Ridge}\),는 아래와 같았다.

\[ \hat{\beta}_{Ridge}=\left(X^{T}X+\lambda I\right)^{-1}X^{T}y \]

위의 식을 우리의 데이터에 적용해 보도록하자. 먼저 능선회귀의 \(\underline{y}\)는 평균으로 centered 된 MEDV변수가 될 것이다.

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

## [1]  8.883944  1.083944  5.483944 -1.116056 -6.616056  1.083944

다음으로 데이터 매트릭스 \(\mathbf{X}\) 역시 각 변수의 평균으로 centered되어 있는 것을 가정하고 있으므로 다음과 같이 설정해 줄 수 있다.

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

## [1] 355  13

kable(head(X), format = "markdown", padding = 2)

먼저 \(\lambda\)가 정해져 있을때 \(\hat{\beta}_{Ridge}\)의 값을 구해보자. 예를 들어 \(\lambda\)의 값이 1일 경우 \(\hat{\beta}_{Ridge}\)의 값은 다음과 같이 구할 수 있다.

lambda <- 1
a <- t(X) %*% X + lambda * diag(1, 13)
b <- t(X) %*% y
b_hat <-  as.numeric(solve(a,b))
sol_ridge <- b_hat
names(sol_ridge) <- colnames(train_Data)[-1]
sol_ridge

##        CRIM          ZN       INDUS        CHAS         NOX          RM 
## -0.08502257  0.05192394  0.01618251  2.30476609 -9.22756880  3.20167465 
##         AGE         DIS         RAD         TAX     PTRATIO           B 
## -0.01147562 -1.39125382  0.30277584 -0.01315725 -0.90526877  0.00918277 
##       LSTAT 
## -0.60098466

위의 R코드에서 solve() 함수의 경우는 선형 연립방정식을 풀어주는 내장함수로, \(\hat{\beta}_{Ridge}\)값을 구하는데에 있어 역행력을 구하여 곱해주는 대신 사용하였다. 위의 식을 다른 경우에도 사용할 수 있도록 함수로 정의해놓자.

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
}

만약 \(\lambda\)의 값이 0일때는 이론상 패널티가 부여되지 않는 모델인 OLS의 베타값과 같아야한다. 우리가 계산한 결과가 회귀분석의 결과와 같은지 확인해보도록 하자.

my_sol <- get_ridgeSol(0, X, y)
my_sol

##          CRIM            ZN         INDUS          CHAS           NOX 
##  -0.094852491   0.050420368   0.056713891   2.525930681 -17.739935159 
##            RM           AGE           DIS           RAD           TAX 
##   3.103207060  -0.003888088  -1.512454613   0.321826861  -0.012332120 
##       PTRATIO             B         LSTAT 
##  -1.006813202   0.008438354  -0.594691841

reg_X <- train_Data %>% as.data.frame()
reg_X <- sweep(reg_X, 2, colMeans(reg_X), "-")
r_sol <- lm(paste("MEDV ~", paste(c(colnames(train_Data)[-1], "0"), collapse = "+")), data = reg_X)
r_sol$coefficients

##          CRIM            ZN         INDUS          CHAS           NOX 
##  -0.094852491   0.050420368   0.056713891   2.525930681 -17.739935159 
##            RM           AGE           DIS           RAD           TAX 
##   3.103207060  -0.003888088  -1.512454613   0.321826861  -0.012332120 
##       PTRATIO             B         LSTAT 
##  -1.006813202   0.008438354  -0.594691841

위의 코드를 통하여 우리가 만든 능선회귀 패라미터 추정 함수의 값이 \(\lambda\)가 0일 경우의 이론값인 회귀분석의 값과 일치하는 것을 확인하였다. 그렇다면 능선회귀에서 해를 추정함에 있어서 유클리디안 놈을 기준으로 패널티를 부여하였으므로, \(\lambda\)값이 커짐에 따라서 능선회귀로 추정된 베타값의 유클리디안 놈의 크기 역시 줄어들 것이라 예상할 수 있을 것이다. 다음의 R코드를 통하여 이것을 확인해 보도록 하자.

먼저 0에서부터 100까지의 \(\lambda\)값에 대응하는 ridge regression의 값는 다음과 같이 구할 수 있다.

lambda <- seq(0, 100, by = 2)
beta_hat <- sapply(lambda, get_ridgeSol, x_data = X, y_data = y)
colnames(beta_hat) <- paste("lambda = ", lambda)
kable(beta_hat[,1:4],format = "markdown", paddind = 2)

이것들의 유클리디안 놈을 구하여 그래프를 그려보면 다음과 같다.

f <- function(x){sqrt(sum(x^2))}
beta_norm <- apply(beta_hat, 2, f)
plot(lambda, beta_norm)

예상대로 \(\lambda\)값이 커짐에 따라서 추정된 패라미터의 유클리디안 값이 줄어듬을 볼 수 있다. 또한 주어진 \(\lambda\)값에 따라서 능선회귀의 패러미터 값이 존재한다는 것을 알았다. 그렇다면 어떤 \(\lambda\)값을 베타를 추정하는데에 사용해야 할까? 필자가 아는 한 이론적으로 최적화된 \(\lambda\)값을 구하는 공식은 없는 것으로 알고 있다. 따라서 주어진 데이터를 가장 잘 설명할 수 있는 모델을 만들어내는 \(\lambda\)값을 사용하는데, 이러한 \(\lambda\)를 정하는 방법을 크로스 벨리데이션(Cross validation)이라고 부른다. 이에 대한 내용은 다음 포스팅에서 다뤄보도록 하겠다.