B1. Hồi quy bội

HỒI QUY TUYẾN TÍNH BỘI

1. ĐỌC DỮ LIỆU

library(readxl)
solieu <- read_excel("/Users/hotranhongnga/Downloads/Solieu.xlsx")
str(solieu)

## tibble [18 × 7] (S3: tbl_df/tbl/data.frame)
##  $ Nam: num [1:18] 2004 2005 2006 2007 2008 ...
##  $ Y  : num [1:18] 2523 2299 2641 2677 3706 ...
##  $ X2 : num [1:18] 266699 267923 279877 327492 326956 ...
##  $ X3 : num [1:18] 2.4 3.2 3.4 3 3.5 3.5 3.6 3.5 3.5 3.3 ...
##  $ X4 : num [1:18] 2.81 3.92 3.8 3.8 3.8 3.8 3.8 3.8 3.79 3.67 ...
##  $ X5 : num [1:18] 9.33 5.74 3.95 4.24 4 4 3.7 3.7 3.7 3.73 ...
##  $ X6 : num [1:18] 2.7 5.2 2.9 1.5 1.4 1.8 1.1 1.5 2.9 3.6 ...

2. MÔ HÌNH HỒI QUY TUYẾN TÍNH BỘI

Mô hình hồi quy tuyến tính bội được xây dựng có dạng như sau:

\[Y = \beta_{1} + \beta_{2}X_{2} + \beta_{3}X_{3} + \beta_{4}X_{4} + \beta_{5}X_{5} + \beta_{6}X_{6} + U\]

Trong đó,

• Y: là biến phụ thuộc.

• \(X_{i}\): Các biến độc lập.

• \(\beta_{i}\): Các hệ số hồi quy.

• U: Sai số ngẫu nhiên.

model_hoiquy <- lm(Y ~ X2 + X3 + X4 + X5 + X6, data = solieu)
summary(model_hoiquy)

## 
## Call:
## lm(formula = Y ~ X2 + X3 + X4 + X5 + X6, data = solieu)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -6225.1 -1550.3   368.3  2004.0  4240.9 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)  
## (Intercept)  3.142e+03  3.979e+04   0.079   0.9384  
## X2           2.606e-02  1.059e-02   2.462   0.0299 *
## X3          -2.944e+03  5.756e+03  -0.511   0.6183  
## X4           1.247e+03  6.225e+03   0.200   0.8446  
## X5          -8.905e+01  1.847e+03  -0.048   0.9623  
## X6          -9.374e+02  7.152e+02  -1.311   0.2145  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3303 on 12 degrees of freedom
## Multiple R-squared:  0.8661, Adjusted R-squared:  0.8102 
## F-statistic: 15.52 on 5 and 12 DF,  p-value: 7.055e-05

3. KẾT QUẢ MÔ HÌNH

3.1 Phương trình hồi quy mẫu (SRF)

\[\hat{Y} = 3142.0 + 0.02606 X_{2} - 2944.0 X_{3} + 1247.0 X_{4} - 89.05 X_{5} - 937.4 X_{6}\]

Ý nghĩa của các hệ số:

• \(\hat{\beta_1}\): 3142.0, cho biết giá trị trung bình của Y khi tất cả các biến độc lập bằng 0.

• \(\hat{\beta_2}\): 0.02606, khi \(X_{2}\) tăng 1 đơn vị, Y trung bình sẽ tăng 0.02606 đơn vị, các yếu tố khác không thay đổi.

• \(\hat{\beta_3}\): −2944.0, khi \(X_{3}\) tăng 1 đơn vị, Y trung bình sẽ giảm 2944.0 đơn vị, các yếu tố khác không thay đổi.

• \(\hat{\beta_4}\): 1247.0, khi \(X_{4}\) tăng 1 đơn vị, Y trung bình sẽ tăng 1247.0 đơn vị, các yếu tố khác không thay đổi.

• \(\hat{\beta_5}\): −89.05, khi \(X_{5}\) tăng 1 đơn vị, Y trung bình sẽ giảm 89.05 đơn vị, các yếu tố khác không thay đổi.

• \(\hat{\beta_6}\): −937.4, Khi \(X_{6}\) tăng 1 đơn vị, Y trung bình sẽ giảm 937.4 đơn vị, các yếu tố khác không thay đổi.

3.2 Mức độ phù hợp tổng thể của mô hình

• Hệ số Xác định (\(R^2\))

Multiple R-squared: 0.8661

Adjusted R-squared (Hiệu chỉnh): 0.8102

Kết luận: Mô hình có mức độ phù hợp rất cao. \(R^2=0.8661\) có nghĩa là 86.61% sự biến động của biến phụ thuộc được giải thích bởi sự biến động của các biến độc lập.

• Kiểm định F (F-statistic)

Giả thuyết kiểm định:

\(H_{0}: \beta_{i}=0\), với \(i∈{2,3,4,5,6}\) (Tất cả các hệ số hồi quy của biến độc lập đều bằng 0; Mô hình không phù hợp).

\({H_1}: \existsβ_{i} \ne 0\),với \(i∈{2,3,4,5,6}\) (Tồn tại ít nhất một hệ số hồi quy khác 0; Mô hình phù hợp ).

Dựa trên kết quả chạy mô hình:

F-statistic: 15.52

P-value: 0.00007055

Kết luận:

Với mức ý nghĩa \(α=0.05\), P-value của kiểm định F rất nhỏ (7.055e−05<0.05), do đó ta đủ cơ sở bác bỏ giả thuyết \(H_{0}\). Mô hình hồi quy phù hợp về mặt thống kê và ít nhất một biến độc lập có tác động đáng kể đến Y.

3.3 Ý nghĩa Thống kê của Từng Biến (Kiểm định t)

Với mức ý nghĩa \(α=0.05\)

Biến	Hệ số (Estimate)	P-value (Pr(>	t
X2	2.606e−02	0.0299*	Có ý nghĩa (Vì 0.0299<0.05).
X3	−2.944e+03	0.6183	Không có ý nghĩa
X4	1.247e+03	0.8446	Không có ý nghĩa
X5	−8.905e+01	0.9623	Không có ý nghĩa
X6	−9.374e+02	0.2145	Không có ý nghĩa

Chỉ có biến \(X_{2}\) là biến có ý nghĩa thống kê trong mô hình này (khi tất cả các biến khác cùng tồn tại). Các biến \(X_{3},X_{4},X_{5},X_{6}\) không thể hiện được tác động đáng kể lên Y ở mức ý nghĩa 5%.

4. KIỂM TRA ĐA CỘNG TUYẾN

library(car)

## Loading required package: carData

vif(model_hoiquy)

##        X2        X3        X4        X5        X6 
## 12.688458  4.934315  6.926831 11.992360  1.376346

Đánh giá mc độ đa cộng tuyến

*\(X_{2}\): 12.688 Rất nghiêm trọng (VIF > 10)

*\(X_{3}\): 4.934 Có thể chấp nhận (VIF < 5)

*\(X_{4}\): 6.927 Trung bình đến cao (VIF > 5)

*\(X_{5}\): 11.992 Rất nghiêm trọng (VIF > 10)

*\(X_{6}\): 1.376 Không nghiêm trọng (VIF < 5)

Kết luận:

Mô hình đang bị ảnh hưởng bởi đa cộng tuyến, đặc biệt ở hai biến X2 (VIF ≈12.69) và X5 (VIF ≈11.99). VIF > 10 là dấu hiệu cho thấy phương sai của các hệ số hồi quy bị phóng đại hơn 10 lần so với khi không có đa cộng tuyến.

Điều này giải thích \(R^2\) của mô hình rất (0.8661) nhưng hầu hết các biến lại không có ý nghĩa thống kê trong kiểm định t (P-value >0.05).

5. ƯỚC LƯỢNG KHOẢNG TIN CẬY CHO CÁC THAM SỐ (\(\beta_{i}\))

confint(model_hoiquy, level = 0.95)

##                     2.5 %       97.5 %
## (Intercept) -8.354565e+04 8.982909e+04
## X2           2.996023e-03 4.912197e-02
## X3          -1.548513e+04 9.597061e+03
## X4          -1.231636e+04 1.481062e+04
## X5          -4.112359e+03 3.934259e+03
## X6          -2.495639e+03 6.208415e+02

Quy tắc: Một tham số \(β_{i}\) được coi là có ý nghĩa thống kê tại mức \(α=0.05\) nếu khoảng tin cậy 95% của nó không chứa giá trị 0.

library(knitr)
library(dplyr)

## 
## Attaching package: 'dplyr'

## The following object is masked from 'package:car':
## 
##     recode

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

# 1. Tạo Data Frame
data_phantich <- data.frame(
  Biến = c("X2", "X3", "X4", "X5", "X6"),
  `KTC 95%` = c(
    "[0.0030, 0.0491]",
    "[-15,485.13, 9,597.06]",
    "[-12,316.36, 14,810.62]",
    "[-4,112.36, 3,934.26]",
    "[-2,495.64, 620.84]"
  ),
  `Chứa giá trị 0?` = c(
    "Không",
    "Có",
    "Có",
    "Có",
    "Có"
  ),
  `Kết luận` = c(
    "Có ý nghĩa thống kê (Bác bỏ H0: β2=0)",
    "Không có ý nghĩa thống kê (Chấp nhận H0: β3=0)",
    "Không có ý nghĩa thống kê",
    "Không có ý nghĩa thống kê",
    "Không có ý nghĩa thống kê"
  ),
  stringsAsFactors = FALSE
)

# 2. Hiển thị bảng
# Tham số align='c' (center) và caption (tiêu đề bảng)
kable(data_phantich, 
      align = 'lccc', 
      caption = "Phân tích Khoảng Tin cậy 95% cho các Hệ số Hồi quy")

Phân tích Khoảng Tin cậy 95% cho các Hệ số Hồi quy

Biến	KTC.95.	Chứa.giá.trị.0.	Kết.luận
X2	[0.0030, 0.0491]	Không	Có ý nghĩa thống kê (Bác bỏ H0: β2=0)
X3	[-15,485.13, 9,597.06]	Có	Không có ý nghĩa thống kê (Chấp nhận H0: β3=0)
X4	[-12,316.36, 14,810.62]	Có	Không có ý nghĩa thống kê
X5	[-4,112.36, 3,934.26]	Có	Không có ý nghĩa thống kê
X6	[-2,495.64, 620.84]	Có	Không có ý nghĩa thống kê

Kết luận:

Biến \(X_{2}\): Khoảng tin cậy 95% cho \(\beta_{2}\) là [0.0030,0.0491]. Vì cả hai giới hạn đều mang dấu dương và không chứa giá trị 0, ta kết luận biến \(X_{2}\) có tác động đáng kể lên biến phụ thuộc Y tại mức ý nghĩa 5%. Ta tin tưởng 95% rằng khi \(X_{2}\) tăng 1 đơn vị, Y sẽ tăng trong khoảng từ 0.0030 đến 0.0491 đơn vị.

Các biến \(X_{3},X_{4},X_{5},X_{6}\): Khoảng tin cậy 95% của các biến này đều chứa giá trị 0. Do đó, ta không thể bác bỏ giả thuyết \(H_0:β_{i}=0\). Các biến này không có tác động đáng kể lên Y tại mức ý nghĩa 5%.

6. KIỂM TRA GIẢ ĐỊNH MÔ HÌNH

6.1. Kiểm tra phương sai

Kiểm tra giả định phương sai đồng nhất (Homoscedasticity) của phần dư. Ta sử dụng Kiểm định Breusch-Pagan.

library(lmtest)

## Loading required package: zoo

## 
## Attaching package: 'zoo'

## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric

# Kiểm định Breusch-Pagan
bptest(model_hoiquy)

## 
##  studentized Breusch-Pagan test
## 
## data:  model_hoiquy
## BP = 9.314, df = 5, p-value = 0.09718

• Giả thuyết: \(H_{0}\): Phương sai đồng nhất

\(H_{1}\): Phương sai không đồng nhất

• Kết luận: Với mức ý nghĩa 5%, \(P-value=0.09718 ≥0.05\), chấp nhận \(H_{0}\). Mô hình đạt giả định phương sai đồng nhất.

6.2. Kiểm tra Tự Tương quan

Ta sử dụng Kiểm định Durbin-Watson.

dwtest(model_hoiquy)

## 
##  Durbin-Watson test
## 
## data:  model_hoiquy
## DW = 1.3097, p-value = 0.005148
## alternative hypothesis: true autocorrelation is greater than 0

• Giả thuyết:

\(H_{0}\): Không có tự tương quan bậc nhất

\(H_{1}\): Có tự tương quan bậc nhất

Kết luận: Với mức ý nghĩa 5%, \(P-value=0.005148 < 0.05\), ta bác bỏ \(H_{0}\). Phần dư có tự tương quan.

Đồ thị kiểm định giả định

par(mfrow = c(2, 2))
plot(model_hoiquy)

## Warning in sqrt(crit * p * (1 - hh)/hh): NaNs produced
## Warning in sqrt(crit * p * (1 - hh)/hh): NaNs produced

Đồ thị Normal Q-Q: Để kiểm tra phân phối chuẩn của phần dư

Đồ thị Residuals vs Fitted: Để kiểm tra phương sai đồng nhất.

7. KHẮC PHỤC KHUYẾT TẬT VÀ XÂY DỰNG MÔ HÌNH TỐI ƯU

7.1 Khắc phục Đa Cộng Tuyến

Mô hình gốc bị đa cộng tuyến nghiêm trọng (VIF của X2≈12.69 và X5≈11.99). Loại bỏ biến có VIF cao nhất là \(X_{2}\) và chạy lại mô hình (Mô hình Ràng buộc 1).

Mô hình Ràng buộc 1: Y∼\(X_{3}\)+\(X_{4}\)+\(X_{5}\)+\(X_{6}\)

# Mô hình Ràng buộc 1: Loại bỏ X2
model_toiuu1 <- lm(Y ~ X3 + X4 + X5 + X6, data = solieu)
print("VIF sau khi loại bỏ X2:")

## [1] "VIF sau khi loại bỏ X2:"

library(car)
vif(model_toiuu1)

##       X3       X4       X5       X6 
## 4.032487 2.486278 2.047758 1.118880

Kết luận VIF:

Tất cả các giá trị VIF trong mô hình mới đều nhỏ hơn 5. Điều này cho thấy vấn đề Đa cộng tuyến đã được khắc phục hoàn toàn. Ta chọn model_toiuu1 là mô hình cơ sở để tiếp tục khắc phục tự tương quan.

7.2 Khắc phục tự tương quan

Sử dụng mô hình mới model_toiuu1 để kiểm tra lại tự tương quan (Durbin-Watson).

# cài đặt gói sandwich
if (!requireNamespace("sandwich", quietly = TRUE)) {
  install.packages("sandwich")
}
library(sandwich)
library(lmtest)

# 1. Áp dụng Sai số chuẩn vững Newey-West
se_nw <- NeweyWest(model_toiuu1, lag = 1) 

# 2. kết quả hồi quy với sai số chuẩn
print("Kết quả Hồi quy cuối cùng (Sai số chuẩn vững Newey-West):")

## [1] "Kết quả Hồi quy cuối cùng (Sai số chuẩn vững Newey-West):"

coeftest(model_toiuu1, vcov = se_nw)

## 
## t test of coefficients:
## 
##              Estimate Std. Error t value  Pr(>|t|)    
## (Intercept)  95629.55   13566.54  7.0489 8.684e-06 ***
## X3           -9001.99    7629.09 -1.1800   0.25916    
## X4          -11023.45    5634.71 -1.9563   0.07226 .  
## X5           -4228.73     653.90 -6.4669 2.110e-05 ***
## X6            -175.89    1053.35 -0.1670   0.86995    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Mô hình model_toiuu1: Y∼\(X_{3}\)+\(X_{4}\)+\(X_{5}\)+\(X_{6}\) đã được ước lượng bằng sai số chuẩn vững Newey-West để khắc phục khuyết tật tự tương quan. Kết quả này là kết quả đáng tin cậy của mô hình.

Kết luận: sau khi loại bỏ đa cộng tuyến và sửa tự tương quan, với mức ý nghĩa 5%, biến \(X_{5}\) đã thành biến có tác động mạnh mẽ và ý nghĩa nhất trong mô hình mới, thay thế cho X2 ở mô hình gốc. Biến \(X_{4}\) cũng cho thấy tác động đáng chú ý (với mức ý nghĩa 1%).

7.3 Phương trình Hồi quy Cuối cùng (SRF)

\[\hat{Y} = 95629.55 - 9001.99 X_{3} - 11023.45 X_{4} - 4228.73 X_{5} - 175.89 X_{6}\]

Giải thích kinh tế (Biến X5):

Hệ số \(\hatβ_{5}}=−4228.73\)

Kết luận: Biến \(X_{5}\) có ý nghĩa thống kê và mang dấu âm.

Giải thích: Khi \(X_{5}\) tăng 1 đơn vị, thì biến phụ thuộc Y trung bình sẽ giảm 4228.73 đơn vị, trong khi các yếu tố độc lập khác được giữ không đổi.

8. DỰ BÁO VỚI MÔ HÌNH TỐI ƯU 1

Sử dụng mô hình tối ưu Y∼\(X_{3}\)+\(X_{4}\)+\(X_{5}\)+\(X_{6}\) để thực hiện dự báo.

# 1. Tạo dữ liệu mới
datamoi_dubao <- data.frame(
  X3 = 3.8,  # Giá trị X3 dự kiến
  X4 = 4.0,  # Giá trị X4 dự kiến
  X5 = 4.2,  # Giá trị X5 dự kiến
  X6 = 1.2   # Giá trị X6 dự kiến
)

# 2. Dự báo giá trị Y
dubao_final <- predict(model_toiuu1, 
                        newdata = datamoi_dubao, 
                        interval = "prediction", 
                        level = 0.95)

print("Dự báo giá trị Y (Khoảng tin cậy 95%):")

## [1] "Dự báo giá trị Y (Khoảng tin cậy 95%):"

dubao_final

##         fit       lwr      upr
## 1 -643.5758 -10708.47 9421.317

Giá trị Dự báo (\(\hat{Y}=−643.58\)):

Đây là giá trị trung bình tốt nhất mà mô hình ước tính cho biến phụ thuộc (Y), khi các biến độc lập \(X_{3}\),\(X_{4}\),\(X_{5}\),\(X_{6}\) mang các giá trị mới dự kiến.

Khoảng dự báo 95%:

Khoảng tin cậy 95% cho giá trị Y thực tế là [−10708.47,9421.32]. Ta tin tưởng 95% rằng giá trị thực tế của biến phụ thuộc (Y) tại các điều kiện mới (của \(X_{3}\),\(X_{4}\),\(X_{5}\),\(X_{6}\)) sẽ nằm trong khoảng từ −10708.47 đến 9421.32.