Giới thiệu
Bài này hướng dẫn các bạn thực hiện hồi quy Ridge một cách chi tiết. Nội dung của bài thuộc chương 6 cuốn An Introduction to Statistical Learning của James et al. (2013). Sách các bạn có thể download miễn phí ở đây.
#-----------------------
# Gọi các gói cần thiết:
#-----------------------
rm(list = ls())
library(ISLR)
data("Hitters")
library(tidyverse)
library(magrittr)
#------------------------
# Tiền xử lí số liệu
#------------------------
# Bỏ dữ liệu trống:
Hitters %<>% na.omit()
# Điều tra qua dữ liệu:
Hitters %>% str()## 'data.frame': 263 obs. of 20 variables:
## $ AtBat : int 315 479 496 321 594 185 298 323 401 574 ...
## $ Hits : int 81 130 141 87 169 37 73 81 92 159 ...
## $ HmRun : int 7 18 20 10 4 1 0 6 17 21 ...
## $ Runs : int 24 66 65 39 74 23 24 26 49 107 ...
## $ RBI : int 38 72 78 42 51 8 24 32 66 75 ...
## $ Walks : int 39 76 37 30 35 21 7 8 65 59 ...
## $ Years : int 14 3 11 2 11 2 3 2 13 10 ...
## $ CAtBat : int 3449 1624 5628 396 4408 214 509 341 5206 4631 ...
## $ CHits : int 835 457 1575 101 1133 42 108 86 1332 1300 ...
## $ CHmRun : int 69 63 225 12 19 1 0 6 253 90 ...
## $ CRuns : int 321 224 828 48 501 30 41 32 784 702 ...
## $ CRBI : int 414 266 838 46 336 9 37 34 890 504 ...
## $ CWalks : int 375 263 354 33 194 24 12 8 866 488 ...
## $ League : Factor w/ 2 levels "A","N": 2 1 2 2 1 2 1 2 1 1 ...
## $ Division : Factor w/ 2 levels "E","W": 2 2 1 1 2 1 2 2 1 1 ...
## $ PutOuts : int 632 880 200 805 282 76 121 143 0 238 ...
## $ Assists : int 43 82 11 40 421 127 283 290 0 445 ...
## $ Errors : int 10 14 3 4 25 7 9 19 0 22 ...
## $ Salary : num 475 480 500 91.5 750 ...
## $ NewLeague: Factor w/ 2 levels "A","N": 2 1 2 2 1 1 1 2 1 1 ...
## - attr(*, "na.action")=Class 'omit' Named int [1:59] 1 16 19 23 31 33 37 39 40 42 ...
## .. ..- attr(*, "names")= chr [1:59] "-Andy Allanson" "-Billy Beane" "-Bruce Bochte" "-Bob Boone" ...Hitters %>% dim()## [1] 263 20# Chuyển hóa thành ma trận thông tin đầu vào:
u <- Hitters %>%
model.matrix(Salary ~ ., data = .)
# u là một dạng "ma trận" với 263 dòng và 20 cột:
u %>% dim()## [1] 263 20u %>% str()## num [1:263, 1:20] 1 1 1 1 1 1 1 1 1 1 ...
## - attr(*, "dimnames")=List of 2
## ..$ : chr [1:263] "-Alan Ashby" "-Alvin Davis" "-Andre Dawson" "-Andres Galarraga" ...
## ..$ : chr [1:20] "(Intercept)" "AtBat" "Hits" "HmRun" ...
## - attr(*, "assign")= int [1:20] 0 1 2 3 4 5 6 7 8 9 ...
## - attr(*, "contrasts")=List of 3
## ..$ League : chr "contr.treatment"
## ..$ Division : chr "contr.treatment"
## ..$ NewLeague: chr "contr.treatment"# Loại ra cột đầu tiên của u để có ma trận thông
# tin đầu vào:
x <- u[, -1]
x %>% dim()## [1] 263 19# Biến phụ thuộc (tức Salary):
y <- Hitters %>% pull(Salary)
#----------------------------------
# Thực hiện Ridge Regression
#----------------------------------
# Tạo ra một khoảng lambda từ 10^10 đến 10^-4:
lambda_grid <- 10^seq(7, -4, length = 100)
# Thực hiện Ridge Regression:
library(glmnet)
ridge1 <- glmnet(x, y,
# chọn alpha = 0 cho Ridge.
# Nếu alpha = 1 sẽ là Lasso:
alpha = 0,
# Chuẩn hóa dữ liệu trước khi thực hiện:
standardize = TRUE,
# Với 100 giá trị lambda:
lambda = lambda_grid)
# Ứng với một lambda sẽ có 20 hệ số hồi quy (tính cả hệ số chặn):
he_so_hq <- coef(ridge1)
# Ví dụ, hệ số hồi quy ứng với lambda = 10^7 là:
he_so_hq[, 1]## (Intercept) AtBat Hits HmRun Runs
## 5.357373e+02 5.443222e-05 1.974423e-04 7.955694e-04 3.338681e-04
## RBI Walks Years CAtBat CHits
## 3.526612e-04 4.150630e-04 1.697604e-03 4.673161e-06 1.719785e-05
## CHmRun CRuns CRBI CWalks LeagueN
## 1.296898e-04 3.449997e-05 3.560316e-05 3.766565e-05 -5.789555e-04
## DivisionW PutOuts Assists Errors NewLeagueN
## -7.806370e-03 2.179882e-05 3.560422e-06 -1.661300e-05 -1.144934e-04# Lambda thứ 50 đó có giá trị:
ridge1$lambda[50]## [1] 35.93814# Các hệ số hồi quy ứng với lambda thứ 50:
hq <- he_so_hq[, 50]
hq## (Intercept) AtBat Hits HmRun Runs
## 6.407070e+01 -5.026628e-01 2.320295e+00 -1.396185e+00 1.101672e+00
## RBI Walks Years CAtBat CHits
## 7.703019e-01 3.016126e+00 -7.692944e+00 2.554182e-03 1.199480e-01
## CHmRun CRuns CRBI CWalks LeagueN
## 6.607760e-01 2.553854e-01 2.387831e-01 -2.098318e-01 4.973103e+01
## DivisionW PutOuts Assists Errors NewLeagueN
## -1.209431e+02 2.577556e-01 1.444918e-01 -3.506943e+00 -1.386494e+01# Có thể tính l2 norm tương ứng:
l2_norm <- sqrt(sum((hq[-1])^2))
l2_norm ## [1] 131.8458# Hàm predict() có thể được sử dụng để tính các hệ
# số hồi quy với một lambda bất kì. Như ta biết hồi quy
# OLS chỉ là một trường hợp riêng của Ridge khi lambda = 0:
predict(ridge1, s = 0, type = "coefficients", exact = TRUE)## 20 x 1 sparse Matrix of class "dgCMatrix"
## 1
## (Intercept) 163.93623646
## AtBat -1.97340819
## Hits 7.39040885
## HmRun 3.99960564
## Runs -2.22703202
## RBI -0.93669742
## Walks 6.20539846
## Years -3.56606365
## CAtBat -0.17721984
## CHits 0.20837873
## CHmRun 0.03122581
## CRuns 1.38009798
## CRBI 0.72134189
## CWalks -0.79770449
## LeagueN 63.39218815
## DivisionW -116.97962810
## PutOuts 0.28199310
## Assists 0.37409929
## Errors -3.41730413
## NewLeagueN -25.86532851# Có thể kiểm tra lại với chú ý rằng có một sự sai
# khác nhỏ bởi vì cách thức làm tròn số của hàm glmnet():
ols <- Hitters %>% lm(Salary ~., data = .)
ols$coefficients## (Intercept) AtBat Hits HmRun Runs
## 163.1035878 -1.9798729 7.5007675 4.3308829 -2.3762100
## RBI Walks Years CAtBat CHits
## -1.0449620 6.2312863 -3.4890543 -0.1713405 0.1339910
## CHmRun CRuns CRBI CWalks LeagueN
## -0.1728611 1.4543049 0.8077088 -0.8115709 62.5994230
## DivisionW PutOuts Assists Errors NewLeagueN
## -116.8492456 0.2818925 0.3710692 -3.3607605 -24.7623251# Có thể xem qua mức độ sai khác:
100*(predict(ridge1, s = 0, type = "coefficients", exact = TRUE) / ols$coefficients)## 20 x 1 sparse Matrix of class "dgCMatrix"
## 1
## (Intercept) 100.51050
## AtBat 99.67348
## Hits 98.52870
## HmRun 92.35081
## Runs 93.72202
## RBI 89.63938
## Walks 99.58455
## Years 102.20717
## CAtBat 103.43140
## CHits 155.51701
## CHmRun -18.06411
## CRuns 94.89743
## CRBI 89.30717
## CWalks 98.29141
## LeagueN 101.26641
## DivisionW 100.11158
## PutOuts 100.03568
## Assists 100.81658
## Errors 101.68247
## NewLeagueN 104.45436#--------------
# Tính MSE:
#--------------
# Tính các giá trị ước lượng khi lambda = 0:
ridge_dubao <- predict(ridge1, s = 0, newx = x)
# Viết hàm tính MSE:
mse <- function(dubao, thuc) {
return(mean((dubao - thuc)^2))
}
# Tính MSE:
mse(ridge_dubao, y)## [1] 92024.73# So với cách chúng ta thực hiện OLS thì
# có một sai khác nhỏ do cách thức làm
# tròn khác nhau:
mse(ols$fitted.values, Hitters$Salary)## [1] 92017.87#-------------------------------------------------
# Chúng ta có thể viết một hàm hình ảnh hóa sự
# Biến đổi của các hệ số hồi quy (trừ hệ số chặn)
# khi lambda thay đổi.
#-------------------------------------------------
plot_coef <- function(ridge_obj, ...) {
heso <- coef(ridge_obj)
cac_hs <- heso@Dimnames[[1]]
gia_tri <- heso@x
df1 <- data.frame(Coef = rep(cac_hs, length(lambda_grid)), Value = gia_tri)
tg <- expand.grid(1:nrow(heso), lambda_grid)
df1$Lambda <- tg$Var2
df1 %>% filter(Coef != "(Intercept)") %>%
ggplot(aes(log(Lambda), Value, color = Coef)) +
geom_line(show.legend = FALSE) +
theme_minimal() +
labs(title = "Coefficients V.s Lambda")
}
# Rồi sử dụng hàm vừa viết:
plot_coef(ridge1)
# Hoặc sử dụng hàm sẵn có (nhưng xấu):
plot(ridge1, xvar = "lambda", label = T)
#------------------------------------------------
# Khảo sát biến đổi của MSE khi lambda thay đổi
#------------------------------------------------
# Thực hiện 10 folds cross validation. Nhớ rằng có
# 10*100 = 1000 mô hình được chạy:
set.seed(1707)
k_folds <- cv.glmnet(x, y,
alpha = 0,
lambda = lambda_grid,
nfolds = 10)
# Viết hàm hình ảnh hóa biến đổi của MSE khi lambda thay đổi:
mse_lambda <- function(ridge_cv, ...) {
mse_lam <- data.frame(Lambda = ridge_cv$lambda,
MSE = ridge_cv$cvm)
mse_lam %>% ggplot(aes(log(Lambda), MSE)) +
geom_line() + geom_point(color = "blue", alpha = 0.3) +
geom_point(data = mse_lam %>% filter(MSE == min(MSE)), color = "red", size = 3) +
labs(title = "MSE V.s Lambda",
subtitle = "Note: The red point corresponds to the smallest value of the MSE") +
theme_minimal()
}
# Sử dụng hàm vừa viết:
mse_lambda(k_folds)
# MSE nhỏ nhất và lambda tương ứng là:
min(k_folds$cvm)## [1] 115907.2k_folds$lambda.min## [1] 5.994843# Các hệ số hồi quy Ridge lúc này (tương ứng với lambda thứ 57) là:
he_so_hq[, which.min(k_folds$cvm)]## (Intercept) AtBat Hits HmRun Runs
## 1.407607e+02 -1.466739e+00 5.118644e+00 1.203639e-01 -6.251881e-03
## RBI Walks Years CAtBat CHits
## 2.154588e-01 4.933163e+00 -1.083019e+01 -4.208203e-02 1.939445e-01
## CHmRun CRuns CRBI CWalks LeagueN
## 7.466971e-01 6.051086e-01 3.404316e-01 -5.562237e-01 6.074272e+01
## DivisionW PutOuts Assists Errors NewLeagueN
## -1.234518e+02 2.779227e-01 2.710992e-01 -3.814888e+00 -2.736373e+01# Với giá trị lambda tối ưu ở trên thì MSE cho train data là:
dubao <- predict(ridge1, s = k_folds$lambda.min, newx = x)
mse(dubao, y)## [1] 93931.71# So sánh với OLS:
mse(ols$fitted.values, y)## [1] 92017.87MSE của OLS thấp hơn cho thấy rằng hồi quy OLS tạo ra những giá trị dự báo tốt hơn so với Ridge. Tuy nhiên lưu ý rằng điều này chỉ mới đúng cho CHỈ MỘT bộ dữ liệu. Để so sánh khả năng dụ báo của hai mô hình chúng ta cần một cách tiếp cận khác chính xác hơn: lần lượt cho hai mô hình này dự báo với, chẳng hạn, 100 mẫu khác nhau rồi so sánh MSE trung bình của 100 lần chạy này:
MSE_ridge <- c()
MSE_ols <- c()
for (i in 1:100) {
set.seed(i)
train <- Hitters %>% sample_frac(0.5)
test <- dplyr::setdiff(Hitters, train)
# Chuẩn bị dữ liệu cho Ridge Rigression:
x_ridge_train <- train %>% model.matrix(Salary ~., data = .)
x_ridge_train <- x_ridge_train[, -1]
y_ridge_train <- train$Salary
x_ridge_test <- test %>% model.matrix(Salary ~., data = .)
x_ridge_test <- x_ridge_test[, -1]
y_ridge_test <- test$Salary
# Thực hiện Ridge Regression:
ridge2 <- glmnet(x_ridge_train, y_ridge_train,
alpha = 0, standardize = TRUE,
lambda = k_folds$lambda.min)
# Dự báo từ Ridge:
dubao <- predict(ridge2, s = k_folds$lambda.min, newx = x_ridge_test)
mse_ridge <- mse(dubao, y_ridge_test)
MSE_ridge <- c(MSE_ridge, mse_ridge)
# Thực hiện OLS:
ols2 <- train %>% lm(Salary ~., data = .)
# Dự báo từ OLS:
dubao <- predict(ols2, test %>% select(-Salary))
mse_ols <- mse(dubao, test$Salary)
MSE_ols <- c(MSE_ols, mse_ols)
}
df_ketqua <- data.frame(Model = c(rep("Ridge", 100), rep("OLS", 100)),
MSE = c(MSE_ridge, MSE_ols))
# Hình ảnh hóa dữ liệu:
df_ketqua %>% ggplot(aes(Model, MSE)) + geom_boxplot()
# Hình ảnh này chỉ ra dấu hiệu rằng Ridge tạo ra các giá trị dự báo
# sát hơn với giá trị thực tế. Chúng ta có thể đưa ra bằng chứng thống
# kê cho sự khác biệt trong dự báo của hai mô hình bằng kiểm định t:
t.test(df_ketqua$MSE ~ df_ketqua$Model)##
## Welch Two Sample t-test
##
## data: df_ketqua$MSE by df_ketqua$Model
## t = 2.6479, df = 195.8, p-value = 0.008761
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## 2036.413 13923.688
## sample estimates:
## mean in group OLS mean in group Ridge
## 133205.4 125225.3p-value = 0.008761 < 5% chỉ ra rằng có sự khác biệt về khả năng dự báo của hai cách tiếp cận. Cụ thể hơn, dự báo từ Ridge Regression tạo ra MSE trung bình trên 100 mẫu thử nghiệm là thấp hơn só với MSE trung bình từ cách tiếp cận truyền thống OLS.
