Khai báo thư viện

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
library(ggplot2)
library(corrplot)
## corrplot 0.95 loaded
library(MASS)
## 
## Attaching package: 'MASS'
## The following object is masked from 'package:dplyr':
## 
##     select
library(lmtest)
## Loading required package: zoo
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
library(car)
## Loading required package: carData
## 
## Attaching package: 'car'
## The following object is masked from 'package:dplyr':
## 
##     recode

Đọc dữ liệu

data <- read.csv("California_housing.csv")
head(data)
str(data)
## 'data.frame':    20640 obs. of  10 variables:
##  $ longitude         : num  -122 -122 -122 -122 -122 ...
##  $ latitude          : num  37.9 37.9 37.9 37.9 37.9 ...
##  $ housing_median_age: num  41 21 52 52 52 52 52 52 42 52 ...
##  $ total_rooms       : num  880 7099 1467 1274 1627 ...
##  $ total_bedrooms    : num  129 1106 190 235 280 ...
##  $ population        : num  322 2401 496 558 565 ...
##  $ households        : num  126 1138 177 219 259 ...
##  $ median_income     : num  8.33 8.3 7.26 5.64 3.85 ...
##  $ median_house_value: num  452600 358500 352100 341300 342200 ...
##  $ ocean_proximity   : chr  "NEAR BAY" "NEAR BAY" "NEAR BAY" "NEAR BAY" ...

Bộ dữ liệu có 9 biến định lượng và 1 biến định tính là ocean_proximit. Kiểm tra xem biến định tính này có bao nhiêu biến độc lập.

length(unique(data$ocean_proximity))
## [1] 5
unique(data$ocean_proximity)
## [1] "NEAR BAY"   "<1H OCEAN"  "INLAND"     "NEAR OCEAN" "ISLAND"

Biến ocean proximity có 5 giá trị riêng biệt, biến đổi biến này thành kiểu dữ liệu factor để hồi quy.

data$ocean_proximity <- as.factor(data$ocean_proximity)
str(data)
## 'data.frame':    20640 obs. of  10 variables:
##  $ longitude         : num  -122 -122 -122 -122 -122 ...
##  $ latitude          : num  37.9 37.9 37.9 37.9 37.9 ...
##  $ housing_median_age: num  41 21 52 52 52 52 52 52 42 52 ...
##  $ total_rooms       : num  880 7099 1467 1274 1627 ...
##  $ total_bedrooms    : num  129 1106 190 235 280 ...
##  $ population        : num  322 2401 496 558 565 ...
##  $ households        : num  126 1138 177 219 259 ...
##  $ median_income     : num  8.33 8.3 7.26 5.64 3.85 ...
##  $ median_house_value: num  452600 358500 352100 341300 342200 ...
##  $ ocean_proximity   : Factor w/ 5 levels "<1H OCEAN","INLAND",..: 4 4 4 4 4 4 4 4 4 4 ...

Kiểm tra giá trị khuyết trong 20.640 quang trắc:

colSums(is.na(data))
##          longitude           latitude housing_median_age        total_rooms 
##                  0                  0                  0                  0 
##     total_bedrooms         population         households      median_income 
##                207                  0                  0                  0 
## median_house_value    ocean_proximity 
##                  0                  0

Biến total_bedrooms có 207 diểm dữ liệu bị khuyết, chiếm 1% bộ dữ liệu, nhóm chọn bỏ các quang trắc này.

data_dropna <- na.omit(data)
colSums(!is.na(data_dropna))
##          longitude           latitude housing_median_age        total_rooms 
##              20433              20433              20433              20433 
##     total_bedrooms         population         households      median_income 
##              20433              20433              20433              20433 
## median_house_value    ocean_proximity 
##              20433              20433

###Thống kê mô tả của các biến

summary(data_dropna)
##    longitude         latitude     housing_median_age  total_rooms   
##  Min.   :-124.3   Min.   :32.54   Min.   : 1.00      Min.   :    2  
##  1st Qu.:-121.8   1st Qu.:33.93   1st Qu.:18.00      1st Qu.: 1450  
##  Median :-118.5   Median :34.26   Median :29.00      Median : 2127  
##  Mean   :-119.6   Mean   :35.63   Mean   :28.63      Mean   : 2637  
##  3rd Qu.:-118.0   3rd Qu.:37.72   3rd Qu.:37.00      3rd Qu.: 3143  
##  Max.   :-114.3   Max.   :41.95   Max.   :52.00      Max.   :39320  
##  total_bedrooms     population      households     median_income    
##  Min.   :   1.0   Min.   :    3   Min.   :   1.0   Min.   : 0.4999  
##  1st Qu.: 296.0   1st Qu.:  787   1st Qu.: 280.0   1st Qu.: 2.5637  
##  Median : 435.0   Median : 1166   Median : 409.0   Median : 3.5365  
##  Mean   : 537.9   Mean   : 1425   Mean   : 499.4   Mean   : 3.8712  
##  3rd Qu.: 647.0   3rd Qu.: 1722   3rd Qu.: 604.0   3rd Qu.: 4.7440  
##  Max.   :6445.0   Max.   :35682   Max.   :6082.0   Max.   :15.0001  
##  median_house_value   ocean_proximity
##  Min.   : 14999     <1H OCEAN :9034  
##  1st Qu.:119500     INLAND    :6496  
##  Median :179700     ISLAND    :   5  
##  Mean   :206864     NEAR BAY  :2270  
##  3rd Qu.:264700     NEAR OCEAN:2628  
##  Max.   :500001
numeric_cols <- sapply(data_dropna, is.numeric)

for (col in names(data_dropna)[numeric_cols]) {
  print(
    ggplot(data_dropna, aes(x = .data[[col]])) +
      geom_histogram(bins = 50, fill = "skyblue", color = "white") +
      labs(title = paste("Histogram of", col), x = col)
  )
}

Trong bộ dữ liệu này, ta có thể dễ dàng nhận thấy rằng 2 biến median_house_value và housing_median_age có giá trị max nhiều bất thườn. Kiểm tra thấy rằng có 2054 quan trắc, chiếm 10% bộ dữ liệu. 10% là con số tương đối đáng kể, nhưng vì bộ dữ liệu có hơn 20.000 quan trắc nên nhóm chọn cách loại bỏ các điểm dữ liệu này đi

max_val_1 <- max(data_dropna$median_house_value)
max_val_2 <- max(data_dropna$housing_median_age)
nrow(data_dropna[data_dropna$median_house_value == max_val_1 | data_dropna$housing_median_age == max_val_2, ])
## [1] 2054
data_dropoutlier <- data_dropna[data_dropna$median_house_value != max_val_1,]
data_dropoutlier <- data_dropoutlier[data_dropoutlier$housing_median_age != max_val_2,]
nrow(data_dropoutlier)
## [1] 18379
ggplot(data_dropoutlier, aes(x = median_house_value)) +
  geom_histogram(bins = 50, fill = "skyblue", color = "white") +
  labs(title = "Histogram of median_house_value")

# housing_median_age
ggplot(data_dropoutlier, aes(x = housing_median_age)) +
  geom_histogram(bins = 50, fill = "lightgreen", color = "white") +
  labs(title = "Histogram of housing_median_age")

set.seed(123) 

n <- nrow(data_dropoutlier)
train_index <- sample(1:n, size = 0.8 * n)

train_data <- data_dropoutlier[train_index, ]
test_data  <- data_dropoutlier[-train_index, ]

Vẽ biểu đồ tương quan giữa các biến, ta có thể dễ dàng nhận thấy rằng, biến longtitude và latitude có mối quan hệ tuyến tính khá mạnh với nhau. Bên cạnh đó các biến total_rooms, total_bedrooms, populations và households cũng có quan hệ tuyến tính với nhau.

num_vars <- sapply(train_data, is.numeric)
df_num <- train_data[, num_vars]
plot(df_num)

Vẽ ma trận tương quan để thấy các biến có quan hệ tuyến tính rõ hơn

numeric_vars <- train_data[sapply(train_data, is.numeric)]
corr_matrix <- cor(numeric_vars, use = "complete.obs")
corrplot(corr_matrix,
         method = "color",
         type = "upper",
         tl.col = "black",
         tl.srt = 45,
         addCoef.col = "black",   # <-- HIỂN THỊ GIÁ TRỊ SỐ
         number.cex = 0.7,        # kích thước chữ số
         col = colorRampPalette(c("blue","white","red"))(200))

Hồi quy biến median_house_value theo 9 biến còn lại

full_model <- lm(median_house_value ~.,data = train_data)
summary(full_model)
## 
## Call:
## lm(formula = median_house_value ~ ., data = train_data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -276615  -37328   -8751   26484  576611 
## 
## Coefficients:
##                             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)               -1.931e+06  8.649e+04 -22.331  < 2e-16 ***
## longitude                 -2.256e+04  1.002e+03 -22.523  < 2e-16 ***
## latitude                  -2.044e+04  9.861e+02 -20.724  < 2e-16 ***
## housing_median_age         7.002e+02  4.902e+01  14.284  < 2e-16 ***
## total_rooms               -7.123e+00  8.335e-01  -8.547  < 2e-16 ***
## total_bedrooms             8.766e+01  7.042e+00  12.448  < 2e-16 ***
## population                -2.869e+01  1.055e+00 -27.196  < 2e-16 ***
## households                 4.053e+01  7.483e+00   5.416 6.18e-08 ***
## median_income              3.844e+04  4.286e+02  89.691  < 2e-16 ***
## ocean_proximityINLAND     -4.185e+04  1.737e+03 -24.091  < 2e-16 ***
## ocean_proximityISLAND      1.832e+05  4.189e+04   4.373 1.23e-05 ***
## ocean_proximityNEAR BAY   -1.572e+04  2.069e+03  -7.598 3.19e-14 ***
## ocean_proximityNEAR OCEAN  5.741e+03  1.606e+03   3.574 0.000352 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 59200 on 14690 degrees of freedom
## Multiple R-squared:  0.6236, Adjusted R-squared:  0.6233 
## F-statistic:  2028 on 12 and 14690 DF,  p-value: < 2.2e-16

Kiểm tra vif của model

vif(full_model)
##                         GVIF Df GVIF^(1/(2*Df))
## longitude          16.516648  1        4.064068
## latitude           18.805012  1        4.336475
## housing_median_age  1.321904  1        1.149741
## total_rooms        14.386444  1        3.792947
## total_bedrooms     37.983306  1        6.163060
## population          6.380520  1        2.525969
## households         35.373700  1        5.947579
## median_income       1.908382  1        1.381442
## ocean_proximity     3.725273  4        1.178677

Vì trong ma trận tương quan, latitude có quan hệ tương quan với biến phụ thuộc cao hơn lontitude nên sẽ loại biến longitude. Biến total_bedrooms cũng có điểm VIF cao nhất nên cũng loại để xem model có được cải thiện không.Tiếp tục loại bỏ biến households vì VIF của biến này vẫn cao. Cuối cùng còn lại 6 biến có VIF khoảng từ 1-2, trong đó có total_rooms với population có VIF > 4.

model_1 <- lm(median_house_value ~. - longitude - total_bedrooms - households, data = train_data)
vif(model_1)
##                        GVIF Df GVIF^(1/(2*Df))
## latitude           1.467499  1        1.211404
## housing_median_age 1.291398  1        1.136397
## total_rooms        4.887500  1        2.210769
## population         4.525194  1        2.127250
## median_income      1.331484  1        1.153900
## ocean_proximity    1.734392  4        1.071256
summary(model_1)
## 
## Call:
## lm(formula = median_house_value ~ . - longitude - total_bedrooms - 
##     households, data = train_data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -237820  -39587  -10152   27064  372672 
## 
## Coefficients:
##                             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                4.301e+04  1.043e+04   4.123 3.76e-05 ***
## latitude                   7.123e+02  2.883e+02   2.471 0.013490 *  
## housing_median_age         6.722e+02  5.071e+01  13.256  < 2e-16 ***
## total_rooms                1.037e+01  5.084e-01  20.405  < 2e-16 ***
## population                -1.641e+01  9.296e-01 -17.650  < 2e-16 ***
## median_income              3.377e+04  3.746e+02  90.145  < 2e-16 ***
## ocean_proximityINLAND     -7.527e+04  1.399e+03 -53.788  < 2e-16 ***
## ocean_proximityISLAND      1.999e+05  4.382e+04   4.561 5.13e-06 ***
## ocean_proximityNEAR BAY   -7.585e+03  2.142e+03  -3.541 0.000399 ***
## ocean_proximityNEAR OCEAN  1.326e+04  1.641e+03   8.083 6.82e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 61950 on 14693 degrees of freedom
## Multiple R-squared:  0.5877, Adjusted R-squared:  0.5874 
## F-statistic:  2327 on 9 and 14693 DF,  p-value: < 2.2e-16

Chọn mô hình bằng phương pháp stepwise với chuẩn BIC.Sau khi chạy thử, ta thấy rằng mô hình cho BIC nhỏ nhất khi loại bỏ biến latitude đi. V mô hình tốt nhất là mô hình gồm 5 biến, housing_median_age, total_rooms, population, median_income và ocean_proximity

modZero <- lm(formula = median_house_value ~ 1, data = train_data)
best_model_1 <- MASS::stepAIC(model_1, direction = "both", scope = list(lower = modZero, upper = model_1), k = log(nrow(train_data)), trace = 10)
## Start:  AIC=324556.5
## median_house_value ~ (longitude + latitude + housing_median_age + 
##     total_rooms + total_bedrooms + population + households + 
##     median_income + ocean_proximity) - longitude - total_bedrooms - 
##     households
## 
##                      Df  Sum of Sq        RSS    AIC
## - latitude            1 2.3434e+10 5.6421e+13 324553
## <none>                             5.6398e+13 324556
## - housing_median_age  1 6.7445e+11 5.7072e+13 324722
## - population          1 1.1957e+12 5.7593e+13 324855
## - total_rooms         1 1.5982e+12 5.7996e+13 324958
## - ocean_proximity     4 1.4702e+13 7.1100e+13 327924
## - median_income       1 3.1191e+13 8.7589e+13 331020
## 
## Step:  AIC=324553
## median_house_value ~ housing_median_age + total_rooms + population + 
##     median_income + ocean_proximity
## 
##                      Df  Sum of Sq        RSS    AIC
## <none>                             5.6421e+13 324553
## + latitude            1 2.3434e+10 5.6398e+13 324556
## - housing_median_age  1 6.6929e+11 5.7090e+13 324717
## - population          1 1.2401e+12 5.7661e+13 324863
## - total_rooms         1 1.6256e+12 5.8047e+13 324961
## - ocean_proximity     4 1.6522e+13 7.2943e+13 328291
## - median_income       1 3.1169e+13 8.7590e+13 331010

###Kiểm định mô hình:

Kiểm định ANOVA: Với mức tin cậy 5%, ta có thể kết luận cả năm biến này đều có ý nghĩa thống kê, đóng góp vào mô hình

anova(best_model_1)

kiểm định xem phân phối của sai số có giống với phân phối chuẩn không. Vì shapiro test chỉ dùng cho dữ liệu < 5000, nên trong bài này nhóm dùng kiểm định Kolmogorov–Smirnov.Với p_value rất nhỏ, nên có thể kết luận rằng phân phối của biến phụ thuộc không giống với phân phối chuẩn.

ks.test(residuals(best_model_1),
        "pnorm",
        mean(residuals(best_model_1)),
        sd(residuals(best_model_1)))
## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  residuals(best_model_1)
## D = 0.092895, p-value < 2.2e-16
## alternative hypothesis: two-sided

Kiểm định phương sai đồng nhất Breusch–Pagan. P-value nhỏ -> phương sai không đồng nhất

bptest(best_model_1)
## 
##  studentized Breusch-Pagan test
## 
## data:  best_model_1
## BP = 428.79, df = 8, p-value < 2.2e-16
plot(best_model_1,1)

plot(best_model_1,2)

Biến đổi Box-cox.

boxcox(best_model_1,lambda = seq(0.15, 0.25, length = 10))

Với \(\lambda\) ~ 0.19 mô hình biến đổi thành

best_model_1 <- lm(median_house_value^0.19 ~ housing_median_age  + total_rooms + population + median_income + ocean_proximity, data = train_data)
summary(best_model_1)
## 
## Call:
## lm(formula = median_house_value^0.19 ~ housing_median_age + total_rooms + 
##     population + median_income + ocean_proximity, data = train_data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.3481 -0.4066 -0.0523  0.3469  2.8386 
## 
## Coefficients:
##                             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                8.791e+00  2.533e-02 347.049  < 2e-16 ***
## housing_median_age         4.162e-03  4.902e-04   8.490  < 2e-16 ***
## total_rooms                9.721e-05  4.908e-06  19.806  < 2e-16 ***
## population                -1.432e-04  8.946e-06 -16.011  < 2e-16 ***
## median_income              3.314e-01  3.621e-03  91.497  < 2e-16 ***
## ocean_proximityINLAND     -9.076e-01  1.226e-02 -74.054  < 2e-16 ***
## ocean_proximityISLAND      1.648e+00  4.237e-01   3.889 0.000101 ***
## ocean_proximityNEAR BAY   -7.947e-02  1.872e-02  -4.245 2.20e-05 ***
## ocean_proximityNEAR OCEAN  7.778e-02  1.586e-02   4.903 9.52e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.5991 on 14694 degrees of freedom
## Multiple R-squared:  0.6284, Adjusted R-squared:  0.6282 
## F-statistic:  3107 on 8 and 14694 DF,  p-value: < 2.2e-16
plot(best_model_1,1)

plot(best_model_1,2)

Loại outlier với cook distance:

cooks <- cooks.distance(best_model_1)
threshold <- 4/nrow(train_data)

influential_points <- which(cooks > threshold)
train_removed_cook <- train_data[-influential_points,]
best_model_1 <- lm(median_house_value^0.19 ~ housing_median_age  + total_rooms + population + median_income + ocean_proximity, data = train_removed_cook)

summary(best_model_1)
## 
## Call:
## lm(formula = median_house_value^0.19 ~ housing_median_age + total_rooms + 
##     population + median_income + ocean_proximity, data = train_removed_cook)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.68821 -0.37133 -0.03994  0.33253  2.10089 
## 
## Coefficients:
##                             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                8.812e+00  2.355e-02 374.108  < 2e-16 ***
## housing_median_age         3.410e-03  4.392e-04   7.765 8.74e-15 ***
## total_rooms                1.411e-04  5.276e-06  26.736  < 2e-16 ***
## population                -2.167e-04  9.808e-06 -22.089  < 2e-16 ***
## median_income              3.267e-01  3.368e-03  96.983  < 2e-16 ***
## ocean_proximityINLAND     -9.487e-01  1.093e-02 -86.807  < 2e-16 ***
## ocean_proximityNEAR BAY   -1.291e-01  1.698e-02  -7.601 3.13e-14 ***
## ocean_proximityNEAR OCEAN  9.171e-02  1.424e-02   6.442 1.21e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.5123 on 13969 degrees of freedom
## Multiple R-squared:  0.7088, Adjusted R-squared:  0.7087 
## F-statistic:  4858 on 7 and 13969 DF,  p-value: < 2.2e-16
plot(best_model_1,4)

kiểm định lại:

ks.test(residuals(best_model_1),
        "pnorm",
        mean(residuals(best_model_1)),
        sd(residuals(best_model_1)))
## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  residuals(best_model_1)
## D = 0.032902, p-value = 1.441e-13
## alternative hypothesis: two-sided
bptest(best_model_1)
## 
##  studentized Breusch-Pagan test
## 
## data:  best_model_1
## BP = 324.52, df = 7, p-value < 2.2e-16

Tuy p_value vẫn rất nhỏ nhưng D-statistic đã giảm 1 nửa từ 0.09 xuống còn 0.03 (dữ liệu đã gần phân phối chuẩn hơn) và BP giảm từ 428 xuống còn 324, có cải thiện phương sai đồng nhất. Nhìn chung biến đổi box-cox đã cải thiện phân phối của sai số giống với phân phối chuẩn hơn và tính đồng nhất phương sai đã được cải thiện.

Vậy với hồi quy bội, mô hình được chọn có \(R^2\) = 07088, mô hình giải thích 70 biến phụ thuộc. mô hình có dạng \(y=x^\frac{1}{0.19}\) và các biến độc lập đều có ý nghĩa thống kê.

summary(best_model_1)
## 
## Call:
## lm(formula = median_house_value^0.19 ~ housing_median_age + total_rooms + 
##     population + median_income + ocean_proximity, data = train_removed_cook)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.68821 -0.37133 -0.03994  0.33253  2.10089 
## 
## Coefficients:
##                             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                8.812e+00  2.355e-02 374.108  < 2e-16 ***
## housing_median_age         3.410e-03  4.392e-04   7.765 8.74e-15 ***
## total_rooms                1.411e-04  5.276e-06  26.736  < 2e-16 ***
## population                -2.167e-04  9.808e-06 -22.089  < 2e-16 ***
## median_income              3.267e-01  3.368e-03  96.983  < 2e-16 ***
## ocean_proximityINLAND     -9.487e-01  1.093e-02 -86.807  < 2e-16 ***
## ocean_proximityNEAR BAY   -1.291e-01  1.698e-02  -7.601 3.13e-14 ***
## ocean_proximityNEAR OCEAN  9.171e-02  1.424e-02   6.442 1.21e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.5123 on 13969 degrees of freedom
## Multiple R-squared:  0.7088, Adjusted R-squared:  0.7087 
## F-statistic:  4858 on 7 and 13969 DF,  p-value: < 2.2e-16

Hồi quy Ridge, hồi quy Lasso, hồi quy Elastic Net

library(glmnet)
## Loading required package: Matrix
## Loaded glmnet 4.1-10

khai báo biến y và biến x

y <- train_removed_cook$median_house_value^0.19
x <- model.matrix(
  median_house_value ~ .,
  data = train_removed_cook)[, -1]

y_test <- test_data$median_house_value^0.19
x_test <- model.matrix(
  median_house_value ~ .,
  data = test_data
)[, colnames(x)]

Hồi quy Ridge, cross validation để chọn ra lambda tối ưu

cv_ridge <- cv.glmnet(x, y, alpha = 0)
plot(cv_ridge)

y_pred_ridge <- predict(cv_ridge, newx = x, s = "lambda.min")
r2_ridge <- 1- sum((y-y_pred_ridge)^2) / sum((y-mean(y))^2)
mse_ridge <- mean((y - y_pred_ridge)^2)

mse_ridge
## [1] 0.2460247
r2_ridge
## [1] 0.7269046

Hồi quy Lasso, cross validation để tìm ra lambda tối ưu.

cv_lasso <- cv.glmnet(x, y, alpha = 1)
plot(cv_lasso)

y_pred_lasso <- predict(cv_lasso, newx = x, s = "lambda.min")
r2_lasso <- 1- sum((y-y_pred_lasso)^2) / sum((y-mean(y))^2)
mse_lasso <- mean((y - y_pred_lasso)^2)

mse_lasso
## [1] 0.2363514
r2_lasso
## [1] 0.7376423

Hồi quy Elastic Net

cv_elastic <- cv.glmnet(x, y, alpha = 0.5)
plot(cv_elastic)

y_pred_elastic <- predict(cv_elastic, newx = x, s = "lambda.min")
r2_elastic <- 1- sum((y-y_pred_elastic)^2) / sum((y-mean(y))^2)
mse_elastic <- mean((y - y_pred_elastic)^2)

mse_elastic
## [1] 0.2363743
r2_elastic
## [1] 0.7376168

##Chạy thử mô hình trên tập test

# Predict trên tập test
y_test_pred_trans <- as.numeric(
  predict(cv_ridge, newx = x_test, s = "lambda.1se")
)

#Tính RMSE
rmse_test_ridge_trans <- sqrt(mean((y_test - y_test_pred_trans)^2))

#Tính R^2
r2_test_ridge_trans <- 1 -
  sum((y_test - y_test_pred_trans)^2) /
  sum((y_test - mean(y_test))^2)

rmse_test_ridge_trans
## [1] 0.569859
r2_test_ridge_trans
## [1] 0.6440245
# Predict trên tập test
y_test_pred_trans <- as.numeric(
  predict(cv_lasso, newx = x_test, s = "lambda.1se")
)

#Tính RMSE
rmse_test_lasso_trans <- sqrt(mean((y_test - y_test_pred_trans)^2))

#R^2
r2_test_lasso_trans <- 1 -
  sum((y_test - y_test_pred_trans)^2) /
  sum((y_test - mean(y_test))^2)

rmse_test_lasso_trans
## [1] 0.5616713
r2_test_lasso_trans
## [1] 0.6541803
# Predict trên scale transform
y_test_pred_trans <- as.numeric(
  predict(cv_elastic, newx = x_test, s = "lambda.1se")
)
#Tính RMSE
rmse_test_elastic_trans <- sqrt(mean((y_test - y_test_pred_trans)^2))
#Tính R^2
r2_test_elastic_trans <- 1 -
  sum((y_test - y_test_pred_trans)^2) /
  sum((y_test - mean(y_test))^2)

rmse_test_elastic_trans
## [1] 0.5622278
r2_test_elastic_trans
## [1] 0.6534947
test_data$y_predicted <- predict(best_model_1, newdata = test_data)
test_data$y_predicted <- test_data$y_predicted
test_data
#Kiểm chứng với mô hình best_model
test_data$y_trains <- test_data$median_house_value^0.19
test_data
rmse_test_best_trans <- sqrt(mean((test_data$y_trains - test_data$y_predicted)^2))
r2_test_best_trans <- 1 -
  sum((test_data$y_trains - test_data$y_predicted)^2) /
  sum((test_data$y_trains - mean(test_data$y_predicted))^2)
rmse_test_best_trans
## [1] 0.591993
r2_test_best_trans
## [1] 0.615842
test_data$y_predicted <- predict(best_model_1, newdata = test_data)
test_data$median_house_value_predicted <- test_data$y_predicted^(1/0.19)
plot(test_data$median_house_value_predicted,test_data$median_house_value)
abline(0, 1, col = "red", lwd = 2)

Do mô hình được huấn luyện để dự đoán trên không gian đã biến đổi y^(0.19) nên các giá trị dự đoán thu được từ mô hình chưa phải là giá nhà thực tế. Để chuyển các dự đoán này trở về đơn vị gốc của giá nhà, cần áp dụng phép biến đổi ngược, tức là y^(1/0.19). Bước này đảm bảo rằng RMSE và R^2 được tính trên thang đo giá nhà thực. Sau khi áp dụng phép biến đổi ngược và đánh giá trên tập test, Elastic Net regression cho kết quả tốt nhất với RMSE thấp nhất và R^2 cao nhất. Điều này cho thấy Elastic Net không chỉ dự đoán chính xác hơn mà còn tổng quát hóa tốt, đặc biệt trong bối cảnh dữ liệu có đa cộng tuyến giữa các biến đầu vào. Do đó, Elastic Net được lựa chọn là mô hình cuối cùng cho bài toán dự đoán giá nhà.