THỰC HÀNH TRÊN LỚP

library(tidyverse)
library(lmtest)
library(car)
library(forecast)
library(Ecdat)
BUK <- BudgetUK

1 Tạo function

Viết function để xử lí công việc:

• Mô hình hồi quy

• Kiểm định tự tương quan

• Kiểm định phương sai thay đổi

• Kiểm định có xảy ra đa cộng tuyến hay không

• Dự báo

Giải thích các biến truyền vào: data là dữ liệu đưa vào thực hiện tính toán, dependent_variable là biến phụ thuộc, independent_variables là biến độc lập.

Ở phần thực hiện dự báo đầu vào của em là trích xuất dữ liệu từ tập dữ liệu (data) dòng thứ 10. Dữ liệu này bao gồm cột của biến phụ thuộc (dependent_variable) và các cột biến độc lập (independent_variables).

hoiquy <- function(data, dependent_variable,independent_variables){
  # Tạo mô hình hồi quy 
  hq <- lm(formula = paste(dependent_variable, "~", paste0(independent_variables, collapse = " + ")), data = data)
  
  coefficients <- coef(hq)
  dependent_variable <- as.character(formula(hq))[2]
  independent_variables <- all.vars(formula(hq))[-1]
  
  # Tạo hàm hồi quy tuyến tính mẫu
  ham_hq <- paste(dependent_variable, "=", round(coefficients[1], 3))
  for (i in seq_along(independent_variables)) {
    ham_hq <- paste(ham_hq, "+", round(coefficients[i+1], 3), "*", independent_variables[i])
  }
  
  # Kiểm định phương sai thay đổi
  pstd <- ncvTest(hq)
  
  # Kiểm định tự tương quan
  dw_test <- durbinWatsonTest(hq)
  
  # Kiểm tra số lượng biến độc lập để xác định xem có vẽ đồ thị hay không, nếu là chỉ có một biến độc lập thì vẽ biểu đồ phân tán.
  num_independent_vars <- length(independent_variables)
  scatter_plot <- NULL
  if (num_independent_vars == 1) {
    scatter_plot <- ggplot(data, aes_string(x = independent_variables, y = dependent_variable)) +
      geom_point(color = "skyblue") +
      geom_smooth(method = "lm", col = "red") +
      ggtitle(paste("Scatter plot of", dependent_variable, "and", independent_variables))
  }
  
  # Kiểm tra đa cộng tuyến nếu có ít nhất 2 biến độc lập
  VIF <- NULL
  if (num_independent_vars >= 2){
    VIF <- vif(hq)
  }
  
  # Thực hiện dự báo
  forecast_data <- data[10, c(dependent_variable, independent_variables)]
  forecast_values <- forecast::forecast(hq, newdata = forecast_data)
  
  return(list(
    mohinhhoiquy = summary(hq),
    hamhoiquytuyentinhmau = ham_hq,
    kiemdinhphuongsaithaydoi = pstd,
    kiemdinhtutuongquan = dw_test,
    daconguyen = VIF,
    du_bao = forecast_values,
    scatter_plot = scatter_plot
  ))
}

2 Sử dụng function

Chỉ định biến phụ thuộc và biến độc lập

dependent_variable <- "totexp"
independent_variables <- c("wfood","wfuel","wcloth","walc","wtrans","wother","income")

Tính toán trên function đã tạo

hoiquy(BUK,dependent_variable, independent_variables=c("income"))

## $mohinhhoiquy
## 
## Call:
## lm(formula = paste(dependent_variable, "~", paste0(independent_variables, 
##     collapse = " + ")), data = data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -277.791  -23.539   -6.713   13.287  234.718 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 55.44803    2.42313   22.88   <2e-16 ***
## income       0.31743    0.01623   19.56   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 38.61 on 1517 degrees of freedom
## Multiple R-squared:  0.2014, Adjusted R-squared:  0.2008 
## F-statistic: 382.5 on 1 and 1517 DF,  p-value: < 2.2e-16
## 
## 
## $hamhoiquytuyentinhmau
## [1] "totexp = 55.448 + 0.317 * income"
## 
## $kiemdinhphuongsaithaydoi
## Non-constant Variance Score Test 
## Variance formula: ~ fitted.values 
## Chisquare = 890.1077, Df = 1, p = < 2.22e-16
## 
## $kiemdinhtutuongquan
##  lag Autocorrelation D-W Statistic p-value
##    1      0.02413827      1.950047   0.336
##  Alternative hypothesis: rho != 0
## 
## $daconguyen
## NULL
## 
## $du_bao
##   Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## 1       90.36487 40.84001 139.8897 14.59566 166.1341
## 
## $scatter_plot

Giải thích kết quả:

• Mô hình hồi quy tuyến tính: totexp = 55.44803 + 0.31743 income

Hệ số hồi quy: Hệ số cho “income” là 0.31743. Điều này nghĩa là, cho mỗi đơn vị tăng của “income”, giá trị dự kiến cho “totexp” tăng thêm 0.31743, giả sử tất cả các biến khác không thay đổi.

p-values cho cả hai hệ số đều rất nhỏ, điều này cho thấy cả hai hệ số đều có ý nghĩa thống kê.

R-squared: R-squared cho mô hình là 0.2014. Điều này có nghĩa là mô hình này giải thích được khoảng 20.14% sự biến động trong biến phụ thuộc “totexp”.

F-statistic (Thống kê F): Giá trị thống kê F lớn (382.5) và giá trị p-value liên quan đến nó rất nhỏ (< 2.2e-16), cho thấy mô hình hồi quy này có ý nghĩa thống kê, tức là nó phù hợp với dữ liệu tốt hơn so với mô hình không có bất kỳ biến nào.

Tóm lại, dựa trên kết quả này, chúng ta có thể kết luận rằng “income” có ảnh hưởng đáng kể đến “totexp”. Tuy nhiên, mô hình này chỉ giải thích được khoảng 20.14% sự biến động trong “totexp”, vì vậy có thể có các biến khác cũng ảnh hưởng đến “totexp”.

• Kiểm định phương sai thay đổi

Giá trị Chisquare: 890.1077, Độ tự do (Df): 1, Giá trị p: < 2.22e-16 (giá trị rất nhỏ, gần bằng 0). Giá trị p rất nhỏ, gần bằng 0, cho thấy rằng có thể bác bỏ giả thuyết H0:phương sai sai số không đổi. Có thể kết luận mô hình trên có phương sai sai số thay đổi.

• Kiểm định tự tương quan

Kết quả kiểm định tự tương quan có mục tiêu xem xét xem trong mô hình hồi quy có tự tương quan hay không. Trong trường hợp này, giá trị D-W Statistic gần 2, và giá trị p lớn hơn mức ý nghĩa 0.05, cho thấy không có đủ bằng chứng để bác bỏ giả thuyết rằng trong mô hình hồi quy không tự tương quan.

• Dự báo

Dự báo cho điểm cụ thể của biến phụ thuộc totexp có giá trị 90.36487. Ngoài ra, trong các cột còn lại, “Lo 80” và “Hi 80” đại diện cho giới hạn dưới và giới hạn trên của khoảng tin cậy 80% và “Lo 95” và “Hi 95” đại diện cho giới hạn dưới và giới hạn trên của khoảng tin cậy 95%

Tương tự, các câu lệnh dưới đây em thực hiện với function đã tạo, hàm hồi quy đa biến.

• Hồi quy đa biến với biến phụ thuộc là “totexp” và biến độc lập là “wfood”, “wfuel”, “wcloth”,“walc”

hoiquy(BUK,dependent_variable,independent_variables=c("wfood","wfuel","wcloth","walc"))

## $mohinhhoiquy
## 
## Call:
## lm(formula = paste(dependent_variable, "~", paste0(independent_variables, 
##     collapse = " + ")), data = data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -89.959 -22.290  -5.311  15.280 256.164 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  171.855      4.730  36.331  < 2e-16 ***
## wfood       -170.466      9.381 -18.172  < 2e-16 ***
## wfuel       -205.540     18.535 -11.089  < 2e-16 ***
## wcloth        50.267     10.701   4.697 2.88e-06 ***
## walc          15.224     14.972   1.017    0.309    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 35.84 on 1514 degrees of freedom
## Multiple R-squared:  0.3133, Adjusted R-squared:  0.3115 
## F-statistic: 172.7 on 4 and 1514 DF,  p-value: < 2.2e-16
## 
## 
## $hamhoiquytuyentinhmau
## [1] "totexp = 171.855 + -170.466 * wfood + -205.54 * wfuel + 50.267 * wcloth + 15.224 * walc"
## 
## $kiemdinhphuongsaithaydoi
## Non-constant Variance Score Test 
## Variance formula: ~ fitted.values 
## Chisquare = 254.7772, Df = 1, p = < 2.22e-16
## 
## $kiemdinhtutuongquan
##  lag Autocorrelation D-W Statistic p-value
##    1      0.03232575      1.935105   0.196
##  Alternative hypothesis: rho != 0
## 
## $daconguyen
##    wfood    wfuel   wcloth     walc 
## 1.149276 1.093700 1.216460 1.062243 
## 
## $du_bao
##   Point Forecast    Lo 80    Hi 80   Lo 95    Hi 95
## 1       82.39558 36.39441 128.3967 12.0173 152.7739
## 
## $scatter_plot
## NULL

• Hồi quy đa biến với biến phụ thuộc là “totexp” và biến độc lập là “wfood”, “wfuel”, “wcloth”,“walc”, “wother”

hoiquy(BUK,dependent_variable,independent_variables=c("wfood","wfuel","wcloth","walc","wother"))

## $mohinhhoiquy
## 
## Call:
## lm(formula = paste(dependent_variable, "~", paste0(independent_variables, 
##     collapse = " + ")), data = data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -89.889 -22.353  -5.259  15.162 256.069 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  170.180      7.698  22.107  < 2e-16 ***
## wfood       -168.939     10.895 -15.506  < 2e-16 ***
## wfuel       -204.158     19.206 -10.630  < 2e-16 ***
## wcloth        51.800     12.062   4.294 1.86e-05 ***
## walc          16.471     15.644   1.053    0.293    
## wother         3.030     10.986   0.276    0.783    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 35.85 on 1513 degrees of freedom
## Multiple R-squared:  0.3133, Adjusted R-squared:  0.3111 
## F-statistic: 138.1 on 5 and 1513 DF,  p-value: < 2.2e-16
## 
## 
## $hamhoiquytuyentinhmau
## [1] "totexp = 170.18 + -168.939 * wfood + -204.158 * wfuel + 51.8 * wcloth + 16.471 * walc + 3.03 * wother"
## 
## $kiemdinhphuongsaithaydoi
## Non-constant Variance Score Test 
## Variance formula: ~ fitted.values 
## Chisquare = 254.0881, Df = 1, p = < 2.22e-16
## 
## $kiemdinhtutuongquan
##  lag Autocorrelation D-W Statistic p-value
##    1      0.03230388      1.935148   0.178
##  Alternative hypothesis: rho != 0
## 
## $daconguyen
##    wfood    wfuel   wcloth     walc   wother 
## 1.549368 1.173654 1.544585 1.159041 1.520879 
## 
## $du_bao
##   Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## 1       82.28644 36.26842 128.3045 11.88236 152.6905
## 
## $scatter_plot
## NULL