library(readr)
library(ggplot2)
library(dplyr)
library(skimr)
library(psych)
library(DT)
library(pander)
library(formattable)
library(htmltools)
library(DescTools)
library(epitools)
library(data.table)
library(kableExtra)
library(epiR)
library(readxl)
library(lmtest)
data <- read_excel("So lieu.xlsx")
options(digits = 4)
datatable(data)

1 Mô hình hồi quy bội

Mô hình hồi quy: \[ Y = X_1 +X_2 + X_3+X_4+ X_5 \]

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

summary(OLS)
## 
## Call:
## lm(formula = Y ~ X2 + X3 + X4 + X5 + X6, data = data)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
##  -6225  -1550    368   2004   4241 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)  
## (Intercept)  3.14e+03   3.98e+04    0.08     0.94  
## X2           2.61e-02   1.06e-02    2.46     0.03 *
## X3          -2.94e+03   5.76e+03   -0.51     0.62  
## X4           1.25e+03   6.23e+03    0.20     0.84  
## X5          -8.90e+01   1.85e+03   -0.05     0.96  
## X6          -9.37e+02   7.15e+02   -1.31     0.21  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3300 on 12 degrees of freedom
## Multiple R-squared:  0.866,  Adjusted R-squared:  0.81 
## F-statistic: 15.5 on 5 and 12 DF,  p-value: 7.05e-05

Nhận xét

Mô hình hồi quy:

\[ Y = 3.14e+03 + (2.61e-02)X2 -(2.94e+03)X3 + (1.25e+03)X4 -(8.90e+01)X5 -(9.37e+02)X6 \]

  • \(R^2_{\text{hiệu chỉnh}} = 0.810\) cho biết được mô hình giải thích được khoảng 81.0% sự thay đổi của biến phụ thuộc \(Y\).

  • F-statistic = 15.517 với bậc tự do \((5;\ 12)\) và p-value = 7.05e-05 < 0.01 ⇒ mô hình có ý nghĩa thống kê ở mức tổng thể.

  • Xét riêng từng hệ số:

    • X2: hệ số 0.0261, \(t = 2.46\), \(p = 0.03\) ⇒ có ý nghĩa ở mức 5%. khi X2 tăng 1 đơn vị (giữ các biến khác không đổi), \(Y\) tăng trung bình 0.0261 đơn vị.

    • X3: hệ số −2940, \(p = 0.62\) ⇒ không ý nghĩa.

    • X4: hệ số 1250, \(p = 0.84\) ⇒ không ý nghĩa.

    • X5: hệ số −89, \(p = 0.96\) ⇒ không ý nghĩa.

    • X6: hệ số −937, \(p = 0.21\) ⇒ không ý nghĩa.

2 Khoảng tin cậy cho các tham số

Bài toán kiểm định:

• H₀: βᵢ = 0

• H₁: βᵢ ≠ 0

confint(OLS, level = 0.95)
##                  2.5 %    97.5 %
## (Intercept) -8.355e+04 8.983e+04
## X2           2.996e-03 4.912e-02
## X3          -1.549e+04 9.597e+03
## X4          -1.232e+04 1.481e+04
## X5          -4.112e+03 3.934e+03
## X6          -2.496e+03 6.208e+02

Nhận xét

Dựa vào bảng khoảng tin cậy 95% cho các tham số hồi quy:

  • Khoảng tin cậy của X2 nằm trong khoảng (0.00513; 0.04681) và không chứa 0, do đó hệ số\(\beta_2\) có ý nghĩa thống kê ở mức 5%. Biến X2 có ảnh hưởng rõ ràng đến biến phụ thuộc Y.

  • Khoảng tin cậy của hệ số chặn, X3, X4, X5 và X6 đều chứa giá trị 0, nên các hệ số này không có ý nghĩa thống kê ở mức 5%. Điều này cho thấy các biến X3, X4, X5, X6 không có tác động đáng kể đến Y trong mô hình hiện tại.

3 Kiểm định các giả thiết và mô hình

3.1 Kiểm định đa cộng tuyến

Theo kết quả kiểm định đa cộng tuyến bằng hệ số phóng đại phương sai (VIF), hệ số VIF có giá trị lớn hơn 10, sẽ có vấn đề về đa cộng tuyến, nếu VIF <2 thì ta sẽ kết luận là không có hiện tượng đa cộng tuyến giữa các biến độc lập.

OLS <- lm(Y ~ X2 + X3 + X4+X5 +X6,  data = data)
library(car)
vif <- vif(OLS)
print(vif)
##     X2     X3     X4     X5     X6 
## 12.688  4.934  6.927 11.992  1.376

Nhận xét

Kết quả VIF (từ bảng):

  • X2 = 12.688, X5 = 11.992. Cho thấy rằng mô hình có hai biến này sẽ gây ra hiện tượng đa cộng tuyến nghiêm trọng làm cho mô hình suy luận về hệ số sẽ bị sai lệch.

  • X3 = 4.934 , X4 = 6.927. Với chỉ số VIF có thể coi là có sự đa cộng tuyến nhẹ và có thể không ảnh hưởng đến kết quả suy luận.

  • X6 = 1.376

Nhưng vì X2 có ý nghĩa thống kê cao trong mô hình nên khi khắc phục mô hình ta sẽ giữ lại biến X2 và loại bỏ đi 2 biến có dấu hiệu là đa cộng tuyến là X4 và đa cộng tuyến nghiêm trọng là X5.

3.2 Kiểm định hiện tượng tự tương quan

Bài roán kiểm định

• H₀: mô hình không có hiện tượng tự tương quan giữa các phần dư.

• H₁: mô hình tồn tại hiện tượng tự tương quan.

bgtest(OLS)
## 
##  Breusch-Godfrey test for serial correlation of order up to 1
## 
## data:  OLS
## LM test = 2.9, df = 1, p-value = 0.09

Nhận xét

Giá trị p-value = 0.09 > 0,05 (mức ý nghĩa được chọn), chưa đủ cơ sở bác bỏ H₀. Vậy mô hình không có có hiện tượng tư tương quan.

3.3 Kiểm định phương sai sai số thay đổi

Bài roán kiểm định

• H₀: không có hiện tượng phương sai sai số thay đổi (phương sai sai số đồng nhất).

• H₁: xảy ra hiện tượng phương sai sai số thay đổi (phương sai sai số không đồng nhất).

library(lmtest)
bptest(Y ~ X2 + X3 + X4 +X5+ X6,  data = data, studentize = F)
## 
##  Breusch-Pagan test
## 
## data:  Y ~ X2 + X3 + X4 + X5 + X6
## BP = 8.6, df = 5, p-value = 0.1

Nhận xét

BP= 8.6, p-value= 0.1 > 0.05 (mức ý nghĩa được chọn), chưa đủ cơ sở bác bỏ H₀. Vậy mô hình không xảy ra hiện tượng phương sai sai số thay đổi.

3.4 Kiểm định phân phối chuẩn của phần dư

Bài roán kiểm định

• H₀: phần dư của mô hình tuân theo phân phối chuẩn (không vi phạm giả thiết OLS).

• H₁: phần dư của mô hình không tuân theo phân phối chuẩn.

library(tseries)
jarque.bera.test(residuals(OLS))
## 
##  Jarque Bera Test
## 
## data:  residuals(OLS)
## X-squared = 1.3, df = 2, p-value = 0.5

Nhận xét

Ta có p-value = 0.5 (> 0.05), do đó chưa đủ cơ sở bác bỏ giả thuyết H₀. Vậy phần dư của mô hình tuân theo phân phối chuẩn, mô hình đảm bảo giả thiết OLS về tính chuẩn của sai số.

4 Khắc phục mô hình

Như đã nói ở trên thì chúng ta sẽ bỏ biến X4 và X5 và tiến hành ước lượng mô hình lại

OLS <- lm(Y ~ X2 + X3  +  X6 ,  data = data)
summary(OLS)
## 
## Call:
## lm(formula = Y ~ X2 + X3 + X6, data = data)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
##  -6483  -1433    385   2112   4054 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  3.45e+03   9.14e+03    0.38     0.71    
## X2           2.56e-02   2.80e-03    9.14  2.8e-07 ***
## X3          -1.75e+03   2.53e+03   -0.69     0.50    
## X6          -9.00e+02   5.87e+02   -1.53     0.15    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3070 on 14 degrees of freedom
## Multiple R-squared:  0.865,  Adjusted R-squared:  0.836 
## F-statistic: 29.9 on 3 and 14 DF,  p-value: 2.4e-06

4.1 Kiếm định đa cộng tuyến

vif <- vif(OLS)
print(vif)
##    X2    X3    X6 
## 1.027 1.103 1.075

Như vậy thì mô hình đã được khắc phục hiện tượng và đồng thời cũng làm cho hệ số xác định hiệu chỉnh cũng đã được tăng lên 0.026. Tiếp đến cần phải kiểm định lại các kiểm định khác để xem mô hình có bị sai sót gì không.

4.2 Kiểm định hiện tượng tự tương quan

bgtest(OLS)
## 
##  Breusch-Godfrey test for serial correlation of order up to 1
## 
## data:  OLS
## LM test = 2.5, df = 1, p-value = 0.1

Với giá trị p là 0.1 thì ta thấy được rằng mô hình không có hiện tượng tự tương quan.

4.3 Kiểm định phương sai sai số thay đổi

bptest(Y ~ X2 + X3 + X6,  data = data, studentize = F)
## 
##  Breusch-Pagan test
## 
## data:  Y ~ X2 + X3 + X6
## BP = 6.1, df = 3, p-value = 0.1

Với giá trị p là 0.1 thì ta thấy được rằng mô hình không có hiện tượng phương sai sai số thay đổi.

5 Dự báo mô hình

Để dự báo thì ta sẽ tạo ra một dữ liệu mới như sau:

newdata <- data.frame(X2 = c(1100000,1090246), X3 = c(5,6), X6 = c(4,2.8))
newdata
##        X2 X3  X6
## 1 1100000  5 4.0
## 2 1090246  6 2.8
predict(OLS, newdata = newdata, interval = "prediction")
##     fit  lwr   upr
## 1 19265 6290 32241
## 2 18349 1161 35536

Đây là kết quả dự báo giá trị trung bình của Y và dự báo khoảng của Y:

  • Với X2 = 1100000 , X3 = 5, X4=4 thì ta có giá trị trung bình của Y là 19265 và khoảng ước lượng của Y là \([6290 , 32241]\).
  • Với X2 = 1090246 , X3 = 6, X4=2.8 thì ta có giá trị trung bình của Y là 18349 và khoảng ước lượng của Y là \([1161 , 35536]\).