한치록(2022). 계량경제학강의, 제4판, 박영사, 19장 부록 “머신러닝 예제”, 버전 0.2.2.
이 문서에 제시된 R 명령들은 주로 James, Witten, Hastie and Tibshirani (2013, An Introduction to Statistical Learning: with Applications in R, Springer)을 참고하였으며, 그 책으로 충분하지 않은 부분은 패키지 문서에 기초하여 직접 코딩하였거나 인터넷 검색으로부터 얻은 정보를 확인하여 작성하였다.
회귀(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 |
학습용 데이터 셋(z14
)의 전체 변수들을 예측변수로
사용하여 OLS 추정을 한 결과는 다음과 같다.
## OLS
<- lm(ynext~., data = z14)
ols summary(ols)$r.sq
# [1] 0.9809264
RMSE(z15$ynext, predict(ols, z15)) # RMSE() defined in index11.php
# [1] 51.38724
# random walk, defined in index11.php
rmspe.rw # [1] 53.24273
회귀의 R제곱은 0.981로 매우 높고, 시험용 자료(z15
)에
대한 예측 정확도(RMSE로 측정하며, 작을수록 정확함) 면에서 단순 임의보행
예측 RMSE (53.24273)보다 근소하게 더 낫다(RMSE = 51.38724).
이제 예측변수 부분집합을 선택하는 방법들(subset selection)을 살펴본다. 먼저 BIC에 의한 선택을 살펴보고, 교차검증(CV)은 코딩이 약간 복잡하므로 마지막에 살펴본다.
Best Subset Selection은 leaps 패키지에 구현되어 있다.
regsubsets
명령이 각 k에 해당하는 ’대표선수’를 선발하도록
회귀를 실행한다.
## Best Subset Selection
library(leaps) # install.packages("leaps") if necessary
<- regsubsets(ynext~., data=z14, nvmax=19) reg.full
‘대표선수’ 모형 중 BIC가 가장 작은 모형을 최적으로 간주하고 선택해
보자. 대표선수 모형 각각에 대하여 BIC 등 수치를 구하려면
summary(reg.full)
명령을 사용하면 된다. BIC는 그 결과 중
$bic
에 보관되어 있다. 다음 명령을 보라.
<- summary(reg.full)
reg.summ plot(reg.summ$bic, type='o', ylab = 'BIC')
<- which.min(reg.summ$bic))
(k.best # [1] 2
그림에서도 k = 2가 최적임을
확인할 수 있다. BIC를 최소화하는 k (k.best
)에 해당하는
계수 추정값들은 다음과 같다.
coef(reg.full, k.best)
# (Intercept) aged deathrate
# -4.8587439 9.2829695 0.8152651
선택된 변수는 aged
(고령인구비율)와
deathrate
(2014년 사망률)임을 알 수 있다.
이 두 변수를 사용한 모형이 훈련 데이터(train set)를 맞추는 정도(R제곱)은 다음과 같다.
$rsq[k.best]
reg.summ# [1] 0.9786799
전체 변수를 사용할 때(OLS R제곱 0.9809264)만은 못하지만, 직전 연도 사망률과 고령인구비율 2변수만으로도 학습용 데이터에 대하여 매우 높은 설명력(R제곱 0.9786799)을 제공한다.
이 모형선택 결과를 이용하여 test set (2015년 데이터)의 목표변수를
예측하는 방법으로는, 중복 작업이기는 하나 lm
을 이용하여
다시 회귀한 후 예측하거나, 아니면 reg.full
로부터 k = 2에 해당하는 계수를 추출하여
Xβ̂ 공식에 맞추어
계산하는 법이 있다. 첫 번째 방법은 중복 작업이므로, 불편하지만
regsubsets
에 의하여 이미 추정된 계수들을 사용하여 예측을
하고자 한다. 이 예측은 이 실습에서 반복적으로 등장할 것이므로 다음
함수를 만들어 사용하겠다.
<- function(x, newdata, id) {
predict.regsubsets <- coef(x, id)
bhat if (is.data.frame(newdata)) newdata <- model.matrix(eval(x$call[[2]]), newdata)
as.numeric(newdata[, names(bhat), drop = FALSE] %*% bhat)
}
이제 k.best
를 결정한 후 곧바로
predict(reg.full, z15, k.best)
처럼 사용할 수 있다.
RMSE(z15$ynext, predict(reg.full, z15, k.best)) # RMSE() defined in index11.php
# [1] 48.98381
위 코드가 맞는지 확인하기 위해서 ynext
를
aged
와 deathrate
에 대하여 다시 OLS 회귀를 하고
그 결과를 이용하여 예측하면 똑같은 결과를 얻는다.
<- lm(ynext~aged+deathrate, data=z14)
reg RMSE(z15$ynext, predict(reg, z15)) # RMSE() defined in index11.php
# [1] 48.98381
이상에서는 BIC를 이용하여 최적 변수를 선택하였다. Mallows’s Cp를 사용할 수도 있는데, 그러면 k = 5가 선택되고, 이 경우 test set에 적용한 결과는 다음과 같다.
plot(reg.summ$cp, type='o', ylab = 'Cp')
<- which.min(reg.summ$cp))
(k.cp # [1] 5
RMSE(z15$ynext, predict(reg.full, z15, k.cp)) # RMSE() defined in index11.php
# [1] 49.37972
Adjusted R-squared를 비교하면 다음과 같다.
plot(reg.summ$adjr2, type='o', ylab = 'Adj. R-squared')
<- which.max(reg.summ$adjr2)) # Maximize adjusted R-squared
(k.adjr2 # [1] 15
RMSE(z15$ynext, predict(reg.full, z15, k.adjr2)) # RMSE() defined in index11.php
# [1] 51.46263
Greedy한 Forward Stepwise Selection을 하려면 앞의
regsubsets
에 method = 'forward'
라는 옵션을
주면 된다. 그 나머지는 동일하다. 실제로 실행해 보면, 이 예에서는 선택된
모형이 똑같다(다른 예에서는 안 그럴 수도 있다).
## Forward Stepwise Selection
<- regsubsets(ynext~., data=z14, nvmax=19, method='forward')
reg.full.fwd <- summary(reg.full.fwd)
reg.summ.fwd coef(reg.full.fwd, which.min(reg.summ.fwd$bic))
# (Intercept) aged deathrate
# -4.8587439 9.2829695 0.8152651
Greedy한 Backward Stepwise Selection을 하려면 앞의
regsubsets
에 method = 'backward'
라는 옵션을
주면 된다. 그 나머지는 동일하다. 실제로 실행해 보면, 이 예에서는 선택된
모형이 똑같다(다른 예에서는 안 그럴 수도 있다).
## Backward Stepwise Selection
<- regsubsets(ynext~., data=z14, nvmax=19, method='backward')
reg.full.bwd <- summary(reg.full.bwd)
reg.summ.bwd coef(reg.full.bwd, which.min(reg.summ.bwd$bic))
# (Intercept) aged deathrate
# -4.8587439 9.2829695 0.8152651
참고로, k = 2인 경우에는 best subset selection, forward stepwise selection, backward stepwise selection에서 선택된 모형들이 동일하였다. k = 3에서는 그렇지 않다.
## Compare Best SS, Forward SS, Backward SS
coef(reg.full, 3)
# (Intercept) smoke aged deathrate
# -65.9705911 2.7054220 9.9856840 0.7939802
coef(reg.full.fwd, 3)
# (Intercept) smoke aged deathrate
# -65.9705911 2.7054220 9.9856840 0.7939802
coef(reg.full.bwd, 3)
# (Intercept) hdrink aged deathrate
# -39.3482555 1.9475844 9.7985234 0.8001178
k = 3에서 Best SS = Forward SS ≠ Backward SS 이다.
이제 10-fold CV를 이용하는 것을 실습해 본다. Best Subset Selection으로 실습할 것이다. 먼저 전체 학습용 데이터의 관측치들을 10개로 나눈다.
## Cross Validation
set.seed(1)
<- sample(1:10, nrow(z14), replace = TRUE) group
1~10 그룹별 표본크기는 다음과 같다.
table(group)
# group
# 1 2 3 4 5 6 7 8 9 10
# 21 17 22 21 23 24 22 21 28 24
다음으로 CV error들을 저장할 공간을 마련한다. 전체 변수 개수가
19개(dim(z14)
)에서 열 개수 20개 중 목표변수를 제외한 나머지
19개)이므로 19x10 행렬을 만들어 주자. k행은 변수 개수, j열은 j번째 그룹을 validation set으로
사용할 때의 SSR을 저장한다.
<- 19
kmax <- 10
nfolds <- matrix(NA, kmax, nfolds) # k = 1..19, fold = 1..10
cv.err head(cv.err)
# [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
# [1,] NA NA NA NA NA NA NA NA NA NA
# [2,] NA NA NA NA NA NA NA NA NA NA
# [3,] NA NA NA NA NA NA NA NA NA NA
# [4,] NA NA NA NA NA NA NA NA NA NA
# [5,] NA NA NA NA NA NA NA NA NA NA
# [6,] NA NA NA NA NA NA NA NA NA NA
이제 10-fold CV error를 실제로 구하자. j그룹 관측치들을 제외하고 모형을
추정한 후, 그 중 k개 변수를
사용한 추정 결과를 j그룹
관측치들에 적용한 예측오차 제곱합을 위에서 만든 cv.err
행렬의 (k,j) 원소로
저장한다.
for (j in 1:nfolds) { # 10-fold CV
<- regsubsets(ynext~., data=z14[group != j, ], nvmax=19) # change for forward/backward SS
regs <- model.matrix(ynext~., data=z14[group == j, ])
xval <- z14$ynext[group == j]
yval for (k in 1:kmax) { # number of selected variables
<- predict(regs, xval, k)
yhat <- sum((yval - yhat)^2)
cv.err[k,j]
} }
이제 각 행별로(즉, 각 변수 개수 k별로) 10개 CV error들을 합산하고 이것이 가장 작은 k를 선택하면 된다.
<- which.min(rowSums(cv.err)))
(k.cv.best # [1] 2
k = 2에서 CV 오차가 가장 작다. 다음 그림으로부터 이를 재확인할 수 있다.
plot(rowSums(cv.err), type='o')
k = 2를 결정하였으므로 이제 전체 학습용 데이터를 사용하여 k = 2에 해당하는 모형을 회귀하고 그 회귀결과를 이용하여 예측을 한다. 계수 추정값들과 test set 적용 시 RMSE는 다음과 같다.
<- regsubsets(ynext~., data=z14, nvmax = k.cv.best)
regs coef(regs, k.cv.best)
# (Intercept) aged deathrate
# -4.8587439 9.2829695 0.8152651
RMSE(z15$ynext, predict(regs, z15, k.cv.best)) # RMSE() defined in index11.php
# [1] 48.98381
수학적으로 보면 우변에 상수항만 있는 null 모형도 고려하여야 할 것이나, 실질적으로는 무의미하므로 고려하지 않고자 한다.
## -------------------------------------------------------------------- ## OLS ols <- lm(ynext~., data = z14) summary(ols)$r.sq RMSE(z15$ynext, predict(ols, z15)) rmspe.rw ## Best Subset Selection library(leaps) reg.full <- regsubsets(ynext~., data=z14, nvmax=19) reg.summ <- summary(reg.full) plot(reg.summ$bic, type='o', ylab = 'BIC') (k.best <- which.min(reg.summ$bic)) coef(reg.full, k.best) reg.summ$rsq[k.best] predict.regsubsets <- function(x, newdata, id) { bhat <- coef(x, id) if (is.data.frame(newdata)) newdata <- model.matrix(eval(x$call[[2]]), newdata) as.numeric(newdata[, names(bhat), drop = FALSE] %*% bhat) } RMSE(z15$ynext, predict(reg.full, z15, k.best)) reg <- lm(ynext~aged+deathrate, data=z14) RMSE(z15$ynext, predict(reg, z15)) plot(reg.summ$cp, type='o', ylab = 'Cp') (k.cp <- which.min(reg.summ$cp)) RMSE(z15$ynext, predict(reg.full, z15, k.cp)) plot(reg.summ$adjr2, type='o', ylab = 'Adj. R-squared') (k.adjr2 <- which.max(reg.summ$adjr2)) # Maximize adjusted R-squared RMSE(z15$ynext, predict(reg.full, z15, k.adjr2)) ## Forward Stepwise Selection reg.full.fwd <- regsubsets(ynext~., data=z14, nvmax=19, method='forward') reg.summ.fwd <- summary(reg.full.fwd) coef(reg.full.fwd, which.min(reg.summ.fwd$bic)) ## Backward Stepwise Selection reg.full.bwd <- regsubsets(ynext~., data=z14, nvmax=19, method='backward') reg.summ.bwd <- summary(reg.full.bwd) coef(reg.full.bwd, which.min(reg.summ.bwd$bic)) ## Compare Best SS, Forward SS, Backward SS coef(reg.full, 3) coef(reg.full.fwd, 3) coef(reg.full.bwd, 3) ## Cross Validation set.seed(1) group <- sample(1:10, nrow(z14), replace = TRUE) table(group) kmax <- 19 nfolds <- 10 cv.err <- matrix(NA, kmax, nfolds) # k = 1..19, fold = 1..10 head(cv.err) for (j in 1:nfolds) { # 10-fold CV regs <- regsubsets(ynext~., data=z14[group != j, ], nvmax=19) # change for forward/backward SS xval <- model.matrix(ynext~., data=z14[group == j, ]) yval <- z14$ynext[group == j] for (k in 1:kmax) { # number of selected variables yhat <- predict(regs, xval, k) cv.err[k,j] <- sum((yval - yhat)^2) } } (k.cv.best <- which.min(rowSums(cv.err))) plot(rowSums(cv.err), type='o') regs <- regsubsets(ynext~., data=z14, nvmax = k.cv.best) coef(regs, k.cv.best) RMSE(z15$ynext, predict(regs, z15, k.cv.best)) ## --------------------------------------------------------------------