library(magrittr)
library(knitr)
오늘은 이전 포스팅에 이어 능선회귀의 패러미터를 실제 데이터를 가지고 추정하는 방법에 대하여 알아보도록 하겠다. 시작하기에 앞서 지난 포스팅에 쓰인 R 코드를 먼저 불러오도록 하자.
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")
kable(head(houseData), format = "markdown", padding = 2)
MEDV | CRIM | ZN | INDUS | CHAS | NOX | RM | AGE | DIS | RAD | TAX | PTRATIO | B | LSTAT |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|
24.0 | 0.00632 | 18 | 2.31 | 0 | 0.538 | 6.575 | 65.2 | 4.0900 | 1 | 296 | 15.3 | 396.90 | 4.98 |
21.6 | 0.02731 | 0 | 7.07 | 0 | 0.469 | 6.421 | 78.9 | 4.9671 | 2 | 242 | 17.8 | 396.90 | 9.14 |
34.7 | 0.02729 | 0 | 7.07 | 0 | 0.469 | 7.185 | 61.1 | 4.9671 | 2 | 242 | 17.8 | 392.83 | 4.03 |
33.4 | 0.03237 | 0 | 2.18 | 0 | 0.458 | 6.998 | 45.8 | 6.0622 | 3 | 222 | 18.7 | 394.63 | 2.94 |
36.2 | 0.06905 | 0 | 2.18 | 0 | 0.458 | 7.147 | 54.2 | 6.0622 | 3 | 222 | 18.7 | 396.90 | 5.33 |
28.7 | 0.02985 | 0 | 2.18 | 0 | 0.458 | 6.430 | 58.7 | 6.0622 | 3 | 222 | 18.7 | 394.12 | 5.21 |
데이터 셋 분할하기
능선회귀 모형을 앞에서 얻는 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)
CRIM | ZN | INDUS | CHAS | NOX | RM | AGE | DIS | RAD | TAX | PTRATIO | B | LSTAT |
---|---|---|---|---|---|---|---|---|---|---|---|---|
-3.011493 | 87.14225 | -9.5420282 | -0.0704225 | -0.1411437 | 0.5339127 | -28.324507 | 4.4585803 | -4.132394 | -145.19718 | -3.3329577 | 27.48392 | -8.5229014 |
-2.656613 | -12.85775 | -0.9620282 | -0.0704225 | -0.0081437 | 0.2849127 | 18.475493 | -0.2639197 | -5.132394 | -97.19718 | -0.0329577 | 30.27392 | -3.1929014 |
-2.976493 | 20.14225 | -8.6820282 | -0.0704225 | -0.0801437 | 0.5669127 | 1.475493 | -0.6835197 | -2.132394 | -179.19718 | -0.0329577 | 31.48392 | -4.9429014 |
-2.756433 | -12.85775 | -0.9620282 | -0.0704225 | -0.0081437 | -0.0160873 | 13.975493 | -0.6034197 | -5.132394 | -97.19718 | -0.0329577 | 27.97392 | -4.5729014 |
3.418237 | -12.85775 | 7.2379718 | -0.0704225 | 0.0318563 | 0.1429127 | 5.975493 | -1.6658197 | 14.867606 | 264.80282 | 1.7670423 | -267.46608 | -0.4429014 |
-2.858213 | -12.85775 | -3.4820282 | -0.0704225 | -0.0591437 | 0.1439127 | -16.524507 | 0.6741803 | -4.132394 | -114.19718 | 1.1670423 | 31.48392 | -5.2729014 |
먼저 \(\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)
lambda = 0 | lambda = 2 | lambda = 4 | lambda = 6 | |
---|---|---|---|---|
CRIM | -0.0948525 | -0.0817251 | -0.0792800 | -0.0784323 |
ZN | 0.0504204 | 0.0525327 | 0.0531733 | 0.0535817 |
INDUS | 0.0567139 | 0.0022061 | -0.0088664 | -0.0134199 |
CHAS | 2.5259307 | 2.1723990 | 1.9800356 | 1.8305849 |
NOX | -17.7399352 | -6.2344423 | -3.7806214 | -2.7126075 |
RM | 3.1032071 | 3.2157124 | 3.1909913 | 3.1470349 |
AGE | -0.0038881 | -0.0139280 | -0.0155682 | -0.0159501 |
DIS | -1.5124546 | -1.3484058 | -1.3127815 | -1.2967122 |
RAD | 0.3218269 | 0.2969896 | 0.2937430 | 0.2936018 |
TAX | -0.0123321 | -0.0134949 | -0.0138486 | -0.0140667 |
PTRATIO | -1.0068132 | -0.8715330 | -0.8470100 | -0.8389010 |
B | 0.0084384 | 0.0094381 | 0.0096347 | 0.0097074 |
LSTAT | -0.5946918 | -0.6046993 | -0.6103447 | -0.6151451 |
이것들의 유클리디안 놈을 구하여 그래프를 그려보면 다음과 같다.
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)이라고 부른다. 이에 대한 내용은 다음 포스팅에서 다뤄보도록 하겠다.
'Statistics > Linear models' 카테고리의 다른 글
R의 회귀분석계수 계산과정에 대하여 - gmp 패키지, 촐레스키(Cholesky) 분해, 그리고 QR 분해 (0) | 2017.08.09 |
---|---|
능선회귀분석 with R - 5. Cross Validation을 통한 람다값 정하기 (0) | 2017.06.23 |
능선회귀분석 with R - 4. 베이지안 MAP 추정량으로서의 능선회귀분석 (0) | 2017.06.18 |
능선회귀분석 with R - 3. Coefficient paths 그리기 (0) | 2017.06.17 |
능선회귀분석 with R - 2. 패라미터 추정하기 (0) | 2017.06.14 |
능선회귀분석 with R - 1. 데이터 마련하기 (0) | 2017.06.13 |
능선회귀(Ridge regression)에 대하여 2 - λ 와 β 의 관계 (0) | 2017.05.24 |
능선 회귀(Ridge regression)에 대하여 1 - 해(solution) 유도하기 (0) | 2017.05.24 |