Validation

Sử dụng cách tiếp cận validation để ước lượng tỷ lệ sai số từ các model linear khác nhau của tập dữ liệu Auto.

library(ISLR)
set.seed(1)
#Sample không replace dữ liệu thành các tập train test tỷ lệ 50:50
train = sample(392,196)
#Hồi qui model trên tập train. argument subset dung để lưu index của tập dữ liệu training model
lm.fit <- lm(mpg ~ horsepower, data = Auto, subset = train)
summary(lm.fit)
## 
## Call:
## lm(formula = mpg ~ horsepower, data = Auto, subset = train)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -13.698  -3.085  -0.216   2.680  16.770 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 40.340377   1.002269   40.25   <2e-16 ***
## horsepower  -0.161701   0.008809  -18.36   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.692 on 194 degrees of freedom
## Multiple R-squared:  0.6346, Adjusted R-squared:  0.6327 
## F-statistic: 336.9 on 1 and 194 DF,  p-value: < 2.2e-16
#Dùng lệnh attach để sử dụng trường mà không cần gọi dataframe
attach(Auto)
#Tính MSE trên tập test
mean((mpg-predict(lm.fit,Auto))[-train]^2)
## [1] 26.14142

Ước lượng cho các model đa thức bậc 2 và 3

#Bậc 2
lm.fit2 <- lm(mpg ~ poly(horsepower,2,raw = TRUE), data = Auto, subset = train)
summary(lm.fit2)
## 
## Call:
## lm(formula = mpg ~ poly(horsepower, 2, raw = TRUE), data = Auto, 
##     subset = train)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -14.8507  -2.6415  -0.0274   2.2932  14.7340 
## 
## Coefficients:
##                                    Estimate Std. Error t value Pr(>|t|)
## (Intercept)                      55.8303172  2.7105421  20.597  < 2e-16
## poly(horsepower, 2, raw = TRUE)1 -0.4416096  0.0467737  -9.441  < 2e-16
## poly(horsepower, 2, raw = TRUE)2  0.0011220  0.0001847   6.076 6.45e-09
##                                     
## (Intercept)                      ***
## poly(horsepower, 2, raw = TRUE)1 ***
## poly(horsepower, 2, raw = TRUE)2 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.31 on 193 degrees of freedom
## Multiple R-squared:  0.6933, Adjusted R-squared:  0.6901 
## F-statistic: 218.1 on 2 and 193 DF,  p-value: < 2.2e-16
#Tính MSE trên tập test
mean((mpg-predict(lm.fit2,Auto))[-train]^2)
## [1] 19.82259
#Bậc 3
lm.fit3 <- lm(mpg ~ poly(horsepower,3,raw = TRUE), data = Auto, subset = train)
summary(lm.fit3)
## 
## Call:
## lm(formula = mpg ~ poly(horsepower, 3, raw = TRUE), data = Auto, 
##     subset = train)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -14.8585  -2.5255  -0.0775   2.3203  14.6420 
## 
## Coefficients:
##                                    Estimate Std. Error t value Pr(>|t|)
## (Intercept)                       5.853e+01  7.068e+00   8.282 2.03e-14
## poly(horsepower, 3, raw = TRUE)1 -5.143e-01  1.817e-01  -2.830  0.00514
## poly(horsepower, 3, raw = TRUE)2  1.722e-03  1.461e-03   1.179  0.24000
## poly(horsepower, 3, raw = TRUE)3 -1.531e-06  3.695e-06  -0.414  0.67915
##                                     
## (Intercept)                      ***
## poly(horsepower, 3, raw = TRUE)1 ** 
## poly(horsepower, 3, raw = TRUE)2    
## poly(horsepower, 3, raw = TRUE)3    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.319 on 192 degrees of freedom
## Multiple R-squared:  0.6935, Adjusted R-squared:  0.6888 
## F-statistic: 144.8 on 3 and 192 DF,  p-value: < 2.2e-16
#Tính MSE trên tập test
mean((mpg-predict(lm.fit3,Auto))[-train]^2)
## [1] 19.78252

Kết quả cho thấy MSE bậc 2 lớn hơn MSE bậc 3 nhưng không đáng kể, MSE linear nhỏ hơn MSE của bậc 2. Thực hiện trên validation trên 1 lần chia mẫu khác.

set.seed(2)
train = sample(392,196)

#Linear
lm.fit <- lm(mpg ~ horsepower, data = Auto, subset = train)
attach(Auto)
mean((mpg-predict(lm.fit,Auto))[-train]^2)
## [1] 23.29559
#Bậc 2 (quadratic)
lm.fit2 <- lm(mpg ~ poly(horsepower,2,raw = TRUE), data = Auto, subset = train)
mean((mpg-predict(lm.fit2,Auto))[-train]^2)
## [1] 18.90124
#Bậc 3
lm.fit3 <- lm(mpg ~ poly(horsepower,3,raw = TRUE), data = Auto, subset = train)
mean((mpg-predict(lm.fit3,Auto))[-train]^2)
## [1] 19.2574

Cách chia mẫu thứ 2 lại đưa ra kết quả trái ngược khi bậc 2 có MSE thấp hơn bậc 3 tuy nhiên MSE linear vẫn là cao nhất.

Theo cách tiếp cận validation thông thường này thì kết quả trên tập test cho thấy giá trị MSE trả về của model quadratic là tốt hơn so với model linear. Tuy nhiên không có bằng chứng rõ ràng để khẳng định bậc 3 sẽ tốt hơn bậc 2.

Phương pháp LOOCV (leav-one-out cross-validation)

LOOCV được ước lượng sẵn trong mọi model generalized linear của hàm glm() và tính từ boot::cv.glm().

#Sử dụng glm() để tính được cv
glm.fit=glm(mpg~horsepower, data=Auto)
summary(glm.fit)
## 
## Call:
## glm(formula = mpg ~ horsepower, data = Auto)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -13.5710   -3.2592   -0.3435    2.7630   16.9240  
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 39.935861   0.717499   55.66   <2e-16 ***
## horsepower  -0.157845   0.006446  -24.49   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 24.06645)
## 
##     Null deviance: 23819.0  on 391  degrees of freedom
## Residual deviance:  9385.9  on 390  degrees of freedom
## AIC: 2363.3
## 
## Number of Fisher Scoring iterations: 2
cv.err <- boot::cv.glm(Auto, glm.fit)
cv.err$delta
## [1] 24.23151 24.23114

Kết quả trả về cho thấy ước lượng sai số bình quân của toàn bộ model là 24.231 theo phương pháp LOOCV. Có 2 giá trị trả về là đại điện cho trường hợp có điều chỉnh kích thước mẫu và không điều chỉnh kích thước mẫu.

Viết vòng lặp tính LOOCV cho các bậc từ 1 đến 5

cv.err <- c()
for (i in 1:5){
  glm.fit <- glm(mpg~poly(horsepower,i,raw=TRUE), data=Auto)
  cv.err[i] <- boot::cv.glm(Auto,glm.fit)$delta[1]
}
cv.err
## [1] 24.23151 19.24821 19.33498 19.42443 19.03321

Từ hệ số cv.err của phương pháp LOOCV ta có thể thấy MSE giảm dần khi thay đổi từ model linear sang model quadratic. Và khi bậc của đa thức tăng dần thì MSE của model không có dấu hiệu giảm mạnh đáng kể.

k-fold cross-validation

Chúng ta cũng hoàn toàn có thể tính CV cho k-fold trên chính hàm boot::cv.glm() bằng cách truyền thêm vào argument K số fold cần resample. Chằng hạn ta thực hiện k-fold CV cho các model bậc từ 1 đến 10 đối với dữ liệu Auto.

cv.err.10 <- rep(0,10)
for(i in 1:10){
  glm.fit <- glm(mpg~poly(horsepower,i,raw=TRUE),data=Auto)
  cv.err.10[i] <- boot::cv.glm(Auto, glm.fit, K=10)$delta[1]
}
cv.err.10
##  [1] 24.29029 19.15613 19.31680 19.67053 19.03129 19.03904 18.90680
##  [8] 19.14164 19.49465 19.26043

Chúng ta cũng thấy rằng kết quả trả về cho thấy MSE giảm mạnh khi chuyển từ model linear sang quadratic và không có dẫu hiệu cho thấy bậc cao hơn thì sẽ làm giảm MSE.

Thêm một lưu ý khác đó là mỗi lần chúng ta chạy lại validation model sẽ thu được 1 kết quả khác nhau do quá trình k-fold CV chia mẫu một cách ngẫu nhiên. Chạy lại validation 1 lần nữa để kiểm chứng:

cv.err.10 <- rep(0,10)
for(i in 1:10){
  glm.fit <- glm(mpg~poly(horsepower,i,raw=TRUE),data=Auto)
  cv.err.10[i] <- boot::cv.glm(Auto, glm.fit, K=10)$delta[1]
}
cv.err.10
##  [1] 24.17293 19.18780 19.49280 19.71524 18.87332 18.89430 18.96145
##  [8] 18.93144 18.89767 19.36071

Bootstrap

Một trong những ưu thế của Bootstrap đó là nó áp dụng linh hoạt trong mọi model và mọi tình huống mà không yêu cầu tính toán phức tạp. Chúng ta có thể tự viết các hàm bootstrap cho riêng mình hoặc tiện lợi nhất là sử dụng boot::boot(). Hàm số này sẽ lặp lại quá trình resample mẫu có replacement trên một tập data gốc.

Chúng ta sẽ sử dụng bộ dữ liệu Portfolio về lợi suất của 2 tài sản X, Y của package ISLR để tìm ra tỷ lệ phân bổ đầu tư \(\alpha\) như thế nào sẽ tối thiểu rủi ro danh mục. Theo lý thuyết thì giá trị này bằng:

\[\hat \alpha = \frac{\hat \sigma_{Y} ^2- \hat \sigma_{XY}}{\hat \sigma_{X}^2+\hat \sigma_{Y}^2-2\hat\sigma_{X}\hat\sigma_{Y}}\]

Trong đó \(\hat \sigma^2\) là giá trị ước lượng của các phương sai (variance) và \(\hat\sigma_{XY}\) là hiệp phương sai (covariance) của X và Y.

Đầu tiên ta sẽ viết hàm tính alpha:

alpha.fn <- function(data,index){
  X=data$X[index]
  Y=data$Y[index]
  return (var(Y) - cov(X,Y))/(var(X)+var(Y)-2*cov(X,Y))
}

Test thử hàm alpha.fn cho 100 quan sát được lựa chọn ngẫu nhiên có replacement từ Porfolio.

alpha.fn(Portfolio,sample(100,100,replace = TRUE))
## [1] 0.463758

Sử dụng phương pháp bootstrap để tính \(\alpha\) trên tập dữ liệu gốc là Portfolio 1000 lần:

boot::boot(Portfolio, alpha.fn, R=1000)
## 
## ORDINARY NONPARAMETRIC BOOTSTRAP
## 
## 
## Call:
## boot::boot(data = Portfolio, statistic = alpha.fn, R = 1000)
## 
## 
## Bootstrap Statistics :
##      original      bias    std. error
## t1* 0.6818792 -0.01139355   0.1251636
# alpha.fn là hàm cần tính, R là số lần bootstrap mẫu

Kết quả chỉ ra rằng sử dụng mẫu gốc có \(\hat\alpha\)=0.6819 phương pháp bootstrap đưa ra kết quả ước lượng cho phương sai là \(SE(\hat \alpha)\)=0.1217.

Ứng dụng của bootstrap trong ước lượng độc chính xác của Linear Regression Model

Bootstrap có thể được sử dụng để đánh giá mức độ biến động và dự báo các hệ số ước lượng \(\beta_{0}\)\(\beta_{1}\) của các model thống kê. Chẳng hạn như chúng ta sử dụng phương pháp bootstrap để ước lượng hệ số của chặn và hệ số góc của phương trình hồi qui horsepower theo mpg trong tập dữ liệu Auto.

Tạo phương trình boot.fn() trả về các hệ số hồi qui từ tập dữ liệu được simulation.

boot.fn=function(data,index) {
  return(coef(lm(mpg~horsepower,data = data,subset = index)))
}
#Test thử hàm:
boot.fn(Auto,sample(396,396,replace = TRUE))
## (Intercept)  horsepower 
##  39.3909503  -0.1492495

Bootstrap model 1000 lần kết quả thu được như sau:

boot::boot(Auto, boot.fn, 1000)
## 
## ORDINARY NONPARAMETRIC BOOTSTRAP
## 
## 
## Call:
## boot::boot(data = Auto, statistic = boot.fn, R = 1000)
## 
## 
## Bootstrap Statistics :
##       original        bias    std. error
## t1* 39.9358610  0.0699353317 0.906433806
## t2* -0.1578447 -0.0006407659 0.007901329

Kết quả ước lượng cho thấy mẫu gốc có các hệ số chặn là 39.9359 và hệ số góc là -0.1578. Sai số chuẩn của các hệ số ước lượng này lần lượt là 0.8679 và 0.0076. So sánh với kết quả ước lượng gốc:

summary(lm(data=Auto,mpg~horsepower))$coef
##               Estimate  Std. Error   t value      Pr(>|t|)
## (Intercept) 39.9358610 0.717498656  55.65984 1.220362e-187
## horsepower  -0.1578447 0.006445501 -24.48914  7.031989e-81

Ta thấy có một sự sai khác nhỏ trong phương sai các hệ số ước lượng từ phương pháp bootstrap và từ phương pháp OLS. Liệu rằng đây có phải là điểm yếu của phương pháp bootstrap. Thực tế cho thấy hệ số ước lượng của \(\beta_{0}\)\(\beta_{1}\) thu được từ OLS được ước lượng dựa trên một tham số chưa biết là bình phương sai số ngẫu nhiên của phương trình \(\sigma^2\). Giá trị này được ước lượng dựa trên RSS và công thức này phải giả định rằng mối quan hệ giữa các biến dự báo và biến được dự báo phải là tuyến tính. Chúng ta thực tế không thể biết dạng hàm cụ thể của mối quan hệ này. Do đó trong trường hợp hàm số không có dạng tuyến tính thì RSS sẽ bị phóng đại lên gấp nhiều lần dẫn đến ước lượng cho \(\beta\) sẽ bị chệch. Trong khi đó phương sai thu được từ phương pháp bootstrap không dựa trên những giả thuyết ràng buộc này. Do đó về mặt nào đó thì ước lượng phương sai của bootstrap còn chuẩn xác hơn so với OLS.