Thực hiện Ridge Regression với R
R có nhiều packages thực hiện được Ridge Regression nhưng phổ biến là glmnet và caret. Lời khuyên của người viết bài này: nếu bạn là newbie và bắt đầu Machine Learning thì nên bắt đầu từ glmnet để hiểu rõ hơn chi tiết của cách tiếp cận này dù rằng glmnet rất là “thủ công”. Khi đã hiểu sâu hơn chì chuyển sang dùng caret - một gói khó sử dụng hơn glmnet và cần hiểu kĩ hơn về thống kê cũng như Machine Learning.
Ridge Regression với glmnet
Trong mục này chúng ta sẽ thực hiện hồi quy Ridge với gói glmnet:
library(ISLR)
library(tidyverse)
library(magrittr)
data("Hitters")
Hitters %>% dim()## [1] 322 20# Kiểm tra dữ liệu thiếu:
sapply(Hitters, function(x) {sum(is.na(x))})## AtBat Hits HmRun Runs RBI Walks Years
## 0 0 0 0 0 0 0
## CAtBat CHits CHmRun CRuns CRBI CWalks League
## 0 0 0 0 0 0 0
## Division PutOuts Assists Errors Salary NewLeague
## 0 0 0 0 59 0# Viết hàm chèn dữ liệu thiếu bằng mean:
chen_na <- function(x) {
x[is.na(x)] <- mean(x, na.rm = TRUE)
return(x)
}
Hitters %<>% mutate(Salary = chen_na(Salary))
#--------------------------------
# Chuẩn bị dữ liệu
#--------------------------------
# Chuyển hóa về dạng ma trận cho các Attributes
# bằng hàm model.matrix(). Hàm này nó cũng chuyển
# hóa bất kì biến số định tính nào về định lượng:
x <- model.matrix(Salary ~., Hitters )[, -1]
x %>% str()## num [1:322, 1:19] 293 315 479 496 321 594 185 298 323 401 ...
## - attr(*, "dimnames")=List of 2
## ..$ : chr [1:322] "1" "2" "3" "4" ...
## ..$ : chr [1:19] "AtBat" "Hits" "HmRun" "Runs" ...# Tách riêng ra biến cần dự báo:
y <- Hitters %>% pull(Salary)
#------------------------------
# Thực hiện Ridge Regression
#------------------------------
library(glmnet)
# alpha = 0 tương ứng với Ridge Regression còn
# alpha = 1 tương ứng với The Lasso. Trước hết
# chọn một dải lambda từ 10^10 đến 10^-2:
grid <- 10^seq(10, -2, length = 100)
# Thực hiện Ridge Regression:
ridge.mod <- glmnet(x, y,
alpha = 0,
lambda = grid,
standardize = FALSE)
# Ma trận hệ số hồi quy:
he_so <- coef(ridge.mod)
he_so %>% dim()## [1] 20 100# Hệ số hồi quy khi lambda là 11498 (làm tròn của 11497.57):
he_so[, 50]## (Intercept) AtBat Hits HmRun Runs
## 249.58740086 -1.62200239 4.16894757 0.29591978 0.89523586
## RBI Walks Years CAtBat CHits
## 0.99361304 3.93449997 -0.61182351 -0.25947505 0.58278410
## CHmRun CRuns CRBI CWalks LeagueN
## 0.03931028 1.12230749 0.41978581 -0.47062968 0.21049657
## DivisionW PutOuts Assists Errors NewLeagueN
## -0.91692200 0.23383606 0.30855730 -1.74421221 0.19189953# Có thể sử dụng hàm predict() cho nhiều mục đich, chẳng
# hạn, là tính các hệ số hồi quy khi lambda = 50:
predict(ridge.mod, s = 50, type = "coefficients")[1:20, ]## (Intercept) AtBat Hits HmRun Runs RBI
## 324.0856428 -2.0451418 6.1222913 3.7233693 -0.7940905 -0.1073457
## Walks Years CAtBat CHits CHmRun CRuns
## 5.1870410 -9.1212999 -0.1582074 0.2030202 -0.3960755 1.3118854
## CRBI CWalks LeagueN DivisionW PutOuts Assists
## 0.6036197 -0.6475475 9.9358003 -72.9237996 0.2339469 0.3846452
## Errors NewLeagueN
## -4.2738571 6.0114512#------------------------------------
# Kiểm tra test error
#------------------------------------
set.seed(29)
train <- sample(1:nrow(x), nrow(x)/2)
test <- (-train)
y.test <- y[test]
ridge.mod <- glmnet(x[train, ], y[train],
alpha = 0,
lambda = grid,
thresh = 1e-12)
# Khảo sát sự biến đổi của các hệ số hồi quy theo lambda:
k <- ridge.mod$beta
bien <- k@Dimnames
bien <- bien[[1]]
k1 <- k@x
coef_df <- data.frame(Value = k1, Coeff = rep(bien, 100))
lambda <- c()
for (i in grid) {
lambda <- c(lambda, rep(i, 19))
}
coef_df$Lambda <- lambda
theme_set(theme_minimal())
# Sự biến đổi của các hệ số hồi quy:
coef_df %>% ggplot(aes(log(Lambda), Value, color = Coeff)) +
geom_line(show.legend = FALSE, size = 0.9) +
facet_wrap(~ Coeff, scales = "free")
# Hoặc theo cách khác:
coef_df %>% ggplot(aes(log(Lambda), Value, color = Coeff)) +
geom_line() 
# Hình vừa rồi là tương tự với (nhanh nhưng xấu):
plot(ridge.mod, xvar = "lambda", label = T)
# Dự báo với lambda = 4:
ridge.pred <- predict(ridge.mod, s = 4, newx = x[test, ])
# Tính MSE cho Ridge Regression:
mean((ridge.pred - y.test)^2)## [1] 120374.4# So sánh với OLS:
ols.pred <- predict(ridge.mod, s = 0, newx = x[test,], exact = T)
mean((ols.pred - y.test)^2)## [1] 124136# Kêt quả này cho thấy MSE của Ridge Regression thấp
# hơn MSE của OLS. Nghĩa là Ridge Regression tốt hơn
# cho dự báo mức lương. Ở trên ta đã chọn một giá trị
# lambda bất kì 4. Chúng ta có thể tìm một lambda tối
# ưu như sau:
set.seed(1)
cv.out <- cv.glmnet(x[train, ],
y[train],
alpha = 0,
lambda = grid,
nfolds = 10)
u <- data.frame(lambda = cv.out$lambda, mse = cv.out$cvm)
u %<>% mutate(Min = mse == mse[which.min(mse)])
# Điểm màu xanh là giá trị tối ưu của lambda. Tại giá trị
# lambda này thì MSE là thấp nhất:
u %>%
ggplot(aes(log(lambda), mse, color = Min)) +
geom_line(show.legend = FALSE) +
geom_point(show.legend = FALSE) 
# giá trị lambda làm cho MSE thấp nhất cụ thể là:
lambda_totnhat <- u$lambda[which.min(u$mse)]
lambda_totnhat## [1] 4.641589# Hoặc cách khác:
cv.out$lambda.min## [1] 4.641589# Lúc đó MSE là:
u$mse[which.min(u$mse)]## [1] 89068.87# Tính toán lại khi lambda là 4.641589:
ridge.pred <- predict(ridge.mod, s = lambda_totnhat, newx = x[test, ])
mean((ridge.pred-y.test)^2)## [1] 120190# Dưới đây chúng ta so sánh MSE của hồi quy Ridge
# và OLS trên 100 mẫu thử nghiệm:
mse_ridge <- c()
mse_ols <- c()
for (i in 1:100) {
set.seed(i)
train <- sample(1:nrow(x), nrow(x)/2)
test <- (-train)
y.test <- y[test]
# Ridge Regression:
ridge.pred <- predict(ridge.mod, s = lambda_totnhat, newx = x[test, ])
mse1 <- mean((ridge.pred-y.test)^2)
# OLS:
ols.pred <- predict(ridge.mod, s = 0, newx = x[test,], exact = T)
mse2 <- mean((ols.pred - y.test)^2)
mse_ridge <- c(mse_ridge, mse1)
mse_ols <- c(mse_ols, mse2)
}
# Tạo data frame kết quả các MSE từ hai phương pháp:
ket_qua <- data.frame(MSE = c(mse_ridge, mse_ols), Model = c(rep("Ridge", 100),
rep("OLS", 100)))
# Phân tích hình ảnh:
p1 <- ket_qua %>% ggplot(aes(Model, MSE, color = Model, fill = Model)) +
geom_boxplot(show.legend = FALSE, alpha = 0.3)
p2 <- ket_qua %>% ggplot(aes(MSE)) +
geom_density(alpha = 0.3, color = "red", fill = "red") +
geom_histogram(aes(y = ..density..), fill = "blue", color = "blue", alpha = 0.3) +
facet_wrap(~ Model) +
labs(x = NULL, y = NULL)
library(gridExtra)
grid.arrange(p1, p2)
# Hình ảnh trên (box plot) cho thấy MSE của hồi quy Ridge
# là thấp hơn. Chúng ta có thể đưa ra bằng chứng thống kê
# cho sự khác biệt này về mức độ chính xác trong dự báo
# của hai cách tiếp cận:
# Thực hiện t test:
t.test(ket_qua$MSE ~ ket_qua$Model)##
## Welch Two Sample t-test
##
## data: ket_qua$MSE by ket_qua$Model
## t = 10.751, df = 191.06, p-value < 2.2e-16
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## 18990.36 27524.20
## sample estimates:
## mean in group OLS mean in group Ridge
## 119459.27 96201.99Ridge Regression với caret
Thực hiện hồi quy Ridge với gói caret trước hết các dự báo đầu vào phải ở dạng số và biến được dự báo (Salary) phải ở dạng một vector. Do vậy chuẩn bị số liệu trước khi chạy mô hình cũng có khác biệt:
# Chuyển hóa các biến ở dạng factor về dạng số:
in_put <- Hitters %>%
select(-Salary) %>%
mutate_if(is.factor, as.numeric)
out_put <- Hitters$Salary
# Thiết lập các điều kiện tinh chỉnh:
library(caret)
set.seed(1)
ctrl <- trainControl(method = "cv", number = 10)
# Thực hiện Ridge Regression với gói caret:
set.seed(2)
ridge_caret <- train(in_put, out_put,
method = "ridge",
trControl = ctrl,
preProc = c("center", "scale"),
metric = "RMSE",
tuneGrid = data.frame(.lambda = seq(0, 0.1, length = 15)))
# Nếu chọn RMSE làm tiêu chí lựa chọn mô hình thì
# lambda tối ưu là 0.02857143. Cần chú ý rằng giá
# trị này của lambda không tối ưu R2 (Rsquared):
ridge_caret ## Ridge Regression
##
## 322 samples
## 19 predictor
##
## Pre-processing: centered (19), scaled (19)
## Resampling: Cross-Validated (10 fold)
## Summary of sample sizes: 292, 289, 290, 290, 289, 289, ...
## Resampling results across tuning parameters:
##
## lambda RMSE Rsquared
## 0.000000000 327.8994 0.3863703
## 0.007142857 323.2033 0.4027872
## 0.014285714 322.3775 0.4062850
## 0.021428571 322.1052 0.4075957
## 0.028571429 322.0456 0.4080545
## 0.035714286 322.0851 0.4081229
## 0.042857143 322.1755 0.4080011
## 0.050000000 322.2938 0.4077865
## 0.057142857 322.4280 0.4075295
## 0.064285714 322.5717 0.4072576
## 0.071428571 322.7214 0.4069856
## 0.078571429 322.8751 0.4067218
## 0.085714286 323.0317 0.4064703
## 0.092857143 323.1907 0.4062333
## 0.100000000 323.3518 0.4060113
##
## RMSE was used to select the optimal model using the smallest value.
## The final value used for the model was lambda = 0.02857143.# Lấy các thông tin về kết quả:
u <- ridge_caret$results
# Có thể thấy rõ hơn về sự không trùng khớp về
# giá trị lambda khi tối ưu hóa MSE / RMSE và R2:
u %>% slice(which.min(RMSE))## # A tibble: 1 x 5
## lambda RMSE Rsquared RMSESD RsquaredSD
## <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 0.02857143 322.0456 0.4080545 71.56639 0.1914496u %>% slice(which.max(Rsquared))## # A tibble: 1 x 5
## lambda RMSE Rsquared RMSESD RsquaredSD
## <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 0.03571429 322.0851 0.4081229 72.99153 0.1908537# Hình ảnh hóa biến đổi của MSE, R2 khi
# lambda thay đổi:
u %>%
select(lambda, RMSE, Rsquared) %>%
mutate(MSE = RMSE^2) %>%
gather(Metric, Value, -lambda) %>%
filter(Metric %in% c("MSE", "Rsquared")) %>%
ggplot(aes(lambda, Value, color = Metric)) +
geom_line(show.legend = FALSE) + geom_point(show.legend = FALSE) +
facet_wrap(~ Metric, scale = "free") 