Validation

모델의 성능을 평가할 때에는 항상 검증셋(test dataset)을 이용해야 합니다. training set에서는 항상 더 유연한 모델이 더 정확한 예측 결과를 보여줍니다. 그러나, 훈련셋에서 더 정확한 결과가 검정셋에서의 정확성을 보장해주지는 못합니다. 오히려, 유연한 모델일 수록 검정셋의 noise까지 학습하게 됩니다.즉, overfit 하게 되어 variance 가 커지게 됩니다. 우리는 bias-variance trade-off 를 고려하여 적절한 수준의 flexibility 를 갖는 모델을 선택(model selection)해야 합니다.

검증셋을 이용하여 모델의 성능을 평가하는 방법에는 아래의 세가지가 있습니다.

  1. Validation set approach
  2. LOOCV (Leave-one-out cross-validation)
  3. k-fold CV rule of thumb: k=5 or 10 (bias-variance trade-off)

Auto 데이터셋에서 연비(mpg)를 예측하는 회귀 모델을 만들고, 모델의 성능을 평가해봅시다.

library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
Auto = read.csv("Auto.csv")
str(Auto)
## 'data.frame':    397 obs. of  9 variables:
##  $ mpg         : num  18 15 18 16 17 15 14 14 14 15 ...
##  $ cylinders   : int  8 8 8 8 8 8 8 8 8 8 ...
##  $ displacement: num  307 350 318 304 302 429 454 440 455 390 ...
##  $ horsepower  : Factor w/ 94 levels "?","100","102",..: 17 35 29 29 24 42 47 46 48 40 ...
##  $ weight      : int  3504 3693 3436 3433 3449 4341 4354 4312 4425 3850 ...
##  $ acceleration: num  12 11.5 11 12 10.5 10 9 8.5 10 8.5 ...
##  $ year        : int  70 70 70 70 70 70 70 70 70 70 ...
##  $ origin      : int  1 1 1 1 1 1 1 1 1 1 ...
##  $ name        : Factor w/ 304 levels "amc ambassador brougham",..: 49 36 231 14 161 141 54 223 241 2 ...
Auto$horsepower = as.numeric(as.character(Auto$horsepower))
## Warning: NAs introduced by coercion
Auto = Auto[complete.cases(Auto),]
summary(Auto)
##       mpg          cylinders      displacement     horsepower   
##  Min.   : 9.00   Min.   :3.000   Min.   : 68.0   Min.   : 46.0  
##  1st Qu.:17.00   1st Qu.:4.000   1st Qu.:105.0   1st Qu.: 75.0  
##  Median :22.75   Median :4.000   Median :151.0   Median : 93.5  
##  Mean   :23.45   Mean   :5.472   Mean   :194.4   Mean   :104.5  
##  3rd Qu.:29.00   3rd Qu.:8.000   3rd Qu.:275.8   3rd Qu.:126.0  
##  Max.   :46.60   Max.   :8.000   Max.   :455.0   Max.   :230.0  
##                                                                 
##      weight      acceleration        year           origin     
##  Min.   :1613   Min.   : 8.00   Min.   :70.00   Min.   :1.000  
##  1st Qu.:2225   1st Qu.:13.78   1st Qu.:73.00   1st Qu.:1.000  
##  Median :2804   Median :15.50   Median :76.00   Median :1.000  
##  Mean   :2978   Mean   :15.54   Mean   :75.98   Mean   :1.577  
##  3rd Qu.:3615   3rd Qu.:17.02   3rd Qu.:79.00   3rd Qu.:2.000  
##  Max.   :5140   Max.   :24.80   Max.   :82.00   Max.   :3.000  
##                                                                
##                  name    
##  amc matador       :  5  
##  ford pinto        :  5  
##  toyota corolla    :  5  
##  amc gremlin       :  4  
##  amc hornet        :  4  
##  chevrolet chevette:  4  
##  (Other)           :365
set.seed(1)
train_index = sample(nrow(Auto), round(nrow(Auto)*2/3))
Auto = select(Auto, -name)
train_data = Auto[train_index, ]
test_data = Auto[-train_index, ]
lm.fit = lm(mpg ~ ., data = train_data)
pred = predict(lm.fit, test_data)
obs = test_data$mpg
mse = mean((pred - obs)^2)
mse
## [1] 10.56468

validation set approach 를 10번 반복해서 평균제곱오차의 추정치와 분산을 구해보자.

mse_vset = c()
for (i in 1:10){
  train_index = sample(nrow(Auto), round(nrow(Auto)*2/3), replace = F)
  train_data = Auto[train_index, ]
  test_data = Auto[-train_index, ]
  lm.fit = lm(mpg ~ ., data = train_data)
  pred = predict(lm.fit, test_data)
  obs = test_data$mpg
  mse_vset[i] = mean((pred - obs)^2)
}
mean(mse_vset)
## [1] 11.35968
var(mse_vset)
## [1] 3.006118

LOOCV 를 이용하여 평균제곱오차의 추정치와 분산을 구해보자.

mse_loocv = c()
for (i in 1:nrow(Auto)){
  lm.fit = lm(mpg ~ ., data = Auto[-i,])
  pred = predict(lm.fit, Auto[i,])
  obs = Auto$mpg[i]
  mse_loocv[i] = (pred - obs)^2
}
mean(mse_loocv)
## [1] 11.37113
var(mse_loocv)
## [1] 447.8021

K-fold CV

library(caret)
## Loading required package: lattice
## Loading required package: ggplot2
set.seed(1)
folds = createFolds(1:nrow(Auto), k=10)
mse_kcv = c()
for (i in 1:10){
  lm.fit = lm(mpg ~ ., data = Auto[-folds[[i]],])
  pred = predict(lm.fit, Auto[folds[[i]],])
  obs = Auto[folds[[i]],]$mpg
  mse_kcv[i] = mean((pred - obs)^2)
}
mean(mse_kcv)
## [1] 11.34605
var(mse_kcv)
## [1] 8.582353

위 세가지 검증 기법을 비교해봅시다.

Validation set approach는 training set의 크기가 적어 bias가 큰 경향이 있고, training vs. test dataset 분할의 임의성으로 인해 test error 에 대한 추정치의 변동성이 보통 크게 나타납니다. 즉, bias와 variance가 모두 큽니다.

반면에, LOOCV는 training dataset의 크기가 커서 bias를 줄일 수 있고, k-fold CV 에서의 bias 는 validation set approach와 LOOCV의 사이에 해당합니다.

variance 는 LOOCV 와 k-fold CV 중 어느 것이 더 클까요? 보통 예상과 달리 k-fold CV 이 LOOCV 보다 variance 가 더 작은 경향이 있습니다. LOOCV 는 training dataset 이 서로 매우 유사하므로 모델간의 상관 관계가 매우 높습니다. 모델 간의 상관 관계가 높을 수록 variance 가 더 크게 나타납니다.

Bootstrap

추정치의 정확도, 즉 변동성 평가

표본평균에서 모평균을 추정하는 문제를 기억해봅시다.

다음과 같은 모집단이 있다고 합시다.

population = rnorm(1000, 0, 10) # number=1000, mean=0, sd=10
mean(population)
## [1] -0.4106525
sd(population)
## [1] 10.51136

모집단에서 표본을 추출합니다.

sampl = sample(population, 100, replace = F) 

이제, 모집단의 평균과 표준편차를 모른다고 가정하고, 표본 데이터를 이용해서 모평균과 모표준편차를 추정해봅시다.

먼저, 모평균을 추청해봅시다.
중심극한정리에 따라서 모집단의 평균은 표본평균의 평균과 같습니다.
\[\mu = mean(\bar{X})\] \(\mu: population\ mean\)
\(\bar{X}: sample\ mean\)

즉, 표본평균의 평균이 모평균의 추정치가 되는데, 이 추정치는 얼마나 정확할까요? 우리는 먼저 모평균을 모른다고 가정하였으므로 추정치가 얼마나 정확한지 알수 없습니다. 따라서, 추정치가 얼마나 정확한가에 대한 질문보다는 추정치의 변동성이 얼마나 큰가라는 질문이 더 적절한 질문입니다. 추정치의 변동성을 나타내는 지표가 바로 표준 오차입니다.

standard error (of mean) = standard deviation of estimated population mean \[SE=\frac{\sigma}{\sqrt{n}}\]

\(\sigma = standard\ deviation\ of\ population\) \(n = sample\ size\)

참고로, 신뢰구간은 표준오차와 비슷한 개념인데(추정치의 변동성) 모집단의 모수값이 포함될 가능성이 있는 범위를 나타냅니다.

95% confidence interval of population mean \[\bar{X}-1.96\times\frac{\sigma}{\sqrt{n}} < \mu < \bar{X}+1.96\times\frac{\sigma}{\sqrt{n}}\]

다시, 표준오차로 돌아가서 우리는 모집단의 표준편차 \(\sigma\)를 모른다고 가정하였으므로, 표본의 표준편차 \(\hat{\sigma}\)를 대신 사용합니다.

\[\hat{\sigma} = s\] \(s = standard\ deviation\ of\ sample\)

표본의 표준편차는 다음과 같이 구할 수 있습니다.
\[s^2 = \frac{\sum(X_i-\bar{X})^2}{n-1}\] note: n-1 instead of n 왜, n 대신 n-1을 사용하는 걸까요?
위 식과 같이 n-1을 사용해서 구한 표준편차는 정확히 말하면 표본 표준편차는 아닙니다. 대신 우리는 이를 모집단의 표준편차에 대한 불편 추정량이라고 합니다. 모집단의 표준편차는 항상 표본의 표준편차보다 클 것입니다. 1을 빼주는 것은 우리는 표본의 평균을 알고 있고, 표본의 평균이 주어지면 잔차의 합은 항상 0이기 때문입니다. 즉, 모집단의 분산에 대한 불편 추정량으로서 표본의 분산에 대한 자유도는 n이 아니라 n-1 이 됩니다.
이제 위에서 설명한 식을 이용해서 다음과 같이 표본에서 모평균을 추정할 수 있습니다.

x_bar = mean(sampl)
n = length(sampl)
s = sqrt((sd(sampl)^2)*n/(n-1))
se = s/sqrt(n)
CI = c(x_bar - 1.96*s/sqrt(n), x_bar + 1.96*s/sqrt(n))
se; CI
## [1] 1.146679
## [1] -2.271601  2.223380

다음으로 bootstrapping 기법을 이용해 모평균을 추정해봅시다. bootstrapping 을 위해 반복적으로 적용한 function 을 먼저 정의해주어야 합니다. 여기서는 표본의 평균을 추정하는 함수가 되겠습니다.

fn = function(z, index){
  mean(z[index])
} # function need two arguments, the second one should be index

이제 표본에서 하위 표본을 추출하고(복원을 허락, resample), 해당 하위 표본에 평균을 구하는 함수를 적용합니다.

fn(sampl, sample(length(sampl), length(sampl), replace = T))
## [1] 0.09466536

위 과정을 n번 반복합니다.

# boot 
library(boot)
## 
## Attaching package: 'boot'
## The following object is masked from 'package:lattice':
## 
##     melanoma
res = boot(sampl, fn, R=1000) 
res$t0 
## [1] -0.02411041
# original = the mean of the original sample
mean(res$t)
## [1] -0.00940608
mean(res$t) - res$t0
## [1] 0.01470433
# bias  = the mean of resamples (the esimated mean of the population) - the mean of original sample 
sd(res$t)
## [1] 1.169339
# std. error = the standard deviation of the mean of resamples
boot.ci(res)
## Warning in boot.ci(res): bootstrap variances needed for studentized
## intervals
## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 1000 bootstrap replicates
## 
## CALL : 
## boot.ci(boot.out = res)
## 
## Intervals : 
## Level      Normal              Basic         
## 95%   (-2.3307,  2.2530 )   (-2.4767,  2.2449 )  
## 
## Level     Percentile            BCa          
## 95%   (-2.2931,  2.4285 )   (-2.2130,  2.5167 )  
## Calculations and Intervals on Original Scale

bootstrapping 방법으로 추정한 모평균과 신뢰구간을 위에서 모수적으로 추정한 모평균, 신뢰구간과 비교해봅시다.