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

1 Thống kê mô tả dữ liệu

options(digits = 4) 
library(pastecs)
library(knitr)
kable(stat.desc(data))

	Nam	Y	X2	X3	X4	X5	X6
nbr.val	1.800e+01	1.800e+01	1.800e+01	18.0000	18.0000	18.0000	18.0000
nbr.null	0.000e+00	0.000e+00	0.000e+00	0.0000	0.0000	0.0000	0.0000
nbr.na	0.000e+00	0.000e+00	0.000e+00	0.0000	0.0000	0.0000	0.0000
min	2.004e+03	2.299e+03	2.667e+05	2.4000	2.8100	2.5600	0.6000
max	2.021e+03	2.295e+04	1.069e+06	3.7000	3.9200	9.3300	5.4000
range	1.700e+01	2.065e+04	8.026e+05	1.3000	1.1100	6.7700	4.8000
sum	3.622e+04	1.898e+05	1.054e+07	58.5000	63.0800	71.3600	44.4000
median	2.012e+03	9.430e+03	5.496e+05	3.3000	3.5950	3.7000	2.0500
mean	2.012e+03	1.054e+04	5.858e+05	3.2500	3.5044	3.9644	2.4667
SE.mean	1.258e+00	1.787e+03	6.354e+04	0.0729	0.0798	0.3541	0.3098
CI.mean.0.95	2.655e+00	3.771e+03	1.341e+05	0.1537	0.1684	0.7471	0.6535
var	2.850e+01	5.750e+07	7.268e+10	0.0956	0.1147	2.2573	1.7271
std.dev	5.338e+00	7.583e+03	2.696e+05	0.3092	0.3387	1.5024	1.3142
coef.var	2.700e-03	7.192e-01	4.602e-01	0.0951	0.0966	0.3790	0.5328

2 Mô hình pooled OLS

library(plm)

OLS <- plm(Y ~ X2 + X3 + X4 + X5 + X6,  data = data,  index = c("Nam"),  model = "pooling")

summary(OLS)

## Pooling Model
## 
## Call:
## plm(formula = Y ~ X2 + X3 + X4 + X5 + X6, data = data, model = "pooling", 
##     index = c("Nam"))
## 
## Balanced Panel: n = 18, T = 1, N = 18
## 
## Residuals:
##    Min. 1st Qu.  Median 3rd Qu.    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
## 
## Total Sum of Squares:    9.77e+08
## Residual Sum of Squares: 1.31e+08
## R-Squared:      0.866
## Adj. R-Squared: 0.81
## F-statistic: 15.5171 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\) ⇒ mô hình giải thích được khoảng 81.0% biến thiên 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.

3 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) -7.484e+04 8.112e+04
## X2           5.313e-03 4.681e-02
## X3          -1.423e+04 8.337e+03
## X4          -1.095e+04 1.345e+04
## X5          -3.708e+03 3.530e+03
## X6          -2.339e+03 4.643e+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 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.
Hệ số chặn (Intercept) cũng có khoảng tin cậy bao gồm 0, vì vậy không có ý nghĩa thống kê.

4 Ma trận tương quan giữa các biến

library(corrplot)
data2 <- na.omit(data)
data2 <- data.frame(sapply(data2, as.numeric))  # ép tất cả về dạng số

corr_matrix <- cor(data2)
kable(round(corr_matrix, 4), caption = "Bảng 2: Ma trận tương quan giữa các biến") %>%
  kable_styling(full_width = FALSE, position = "center")

Bảng 2: Ma trận tương quan giữa các biến
	Nam	Y	X2	X3	X4	X5	X6
Nam	1.0000	0.9125	0.9835	-0.0018	-0.5169	-0.6836	-0.0763
Y	0.9125	1.0000	0.9174	-0.1773	-0.5905	-0.5241	-0.1137
X2	0.9835	0.9174	1.0000	-0.1617	-0.6189	-0.5741	0.0259
X3	-0.0018	-0.1773	-0.1617	1.0000	0.6954	-0.6256	-0.2635
X4	-0.5169	-0.5905	-0.6189	0.6954	1.0000	-0.1859	-0.0458
X5	-0.6836	-0.5241	-0.5741	-0.6256	-0.1859	1.0000	0.2416
X6	-0.0763	-0.1137	0.0259	-0.2635	-0.0458	0.2416	1.0000

corrplot(corr_matrix, method = "color", addCoef.col = "black", tl.col = "black")

Nhận xét

Có dấu hiệu đa cộng tuyến trong nhóm {X2, X4, X5} và giữa các cặp {X3-X4, X3-X5} trong mô hình hồi quy, tuy nhiên để đảm bảo độ chính xác cao hơn thì ta sẽ thực hiện kiểm định hệ số phóng đại phương sai VIF

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

5.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 (Field, 2009), 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 <- plm(Y ~ X2 + X3 + X4 + X5 + X6,  data = data,  index = c("Nam"),  model = "pooling")
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 (đa cộng tuyến nghiêm trọng)
X3 = 4.934 , X4 = 6.927 mức vừa (sát ngưỡng)
X6 = 1.376 (không bị hiện tượng đa cộng tuyến)

5.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.

pbgtest(OLS)

## 
##  Breusch-Godfrey/Wooldridge test for serial correlation in panel models
## 
## data:  Y ~ X2 + X3 + X4 + X5 + X6
## chisq = 2.9, df = 1, p-value = 0.09
## alternative hypothesis: serial correlation in idiosyncratic errors

Nhận xét

Chisq = 2.9, 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.

5.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.

5.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ố.

# Ví dụ loại bỏ X5 vì VIF quá cao
OLS_fix1 <- lm(Y ~ X2 + X4 + X6 , data = data)
summary(OLS_fix1)

## 
## Call:
## lm(formula = Y ~ X2 + X4 + X6, data = data)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
##  -6949   -980    349   2023   3947 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.20e+03   1.15e+04    0.10     0.92    
## X2           2.51e-02   3.56e-03    7.06  5.7e-06 ***
## X4          -9.74e+02   2.84e+03   -0.34     0.74    
## X6          -8.01e+02   5.74e+02   -1.39     0.18    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3110 on 14 degrees of freedom
## Multiple R-squared:  0.862,  Adjusted R-squared:  0.832 
## F-statistic: 29.1 on 3 and 14 DF,  p-value: 2.86e-06

vif(OLS_fix1)

##    X2    X4    X6 
## 1.621 1.623 1.002

6 Khắc phục khuyết tật mô hình

OLS1 <- lm(Y ~ X2 + X3 + X6,  data = data)
summary(OLS1)

## 
## 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

Nhận xét

Mức độ phù hợp:
\(R^2 = 0.865\), \(R^2_{adj} = 0.836\)
⟹Tăng so với mô hình trước(≈ 0.810) và mô hình hiện tại giải thích tốt hơn dù dùng ít biến hơn.
Kiểm định F:
\(F = 29.9\) với \(p = 2.4 \times 10^{-6}\)
⟹ mô hình có ý nghĩa thống kê tổng thể.
Ý nghĩa từng biến:
- X2: hệ số dương, \(t = 9.14\), \(p \approx 2.8 \times 10^{-7}\) ⟹ rất có ý nghĩa thống kê.
- X3 (p = 0.50) và X6 (p = 0.15) ⟹ không có ý nghĩa ở mức 5%.
Kết luận:

Việc lược bỏ biến (bỏ X4, X5) giảm hiện tượng đa cộng tuyến và làm \(R^2_{adj}\) tăng, khẳng định X2 là biến giải thích chính của mô hình.

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

vif <- vif(OLS1)
print(vif)

##    X2    X3    X6 
## 1.027 1.103 1.075

Nhận xét

Như vậy, mô hình sau khi khắc phục đã loại bỏ được các khuyết tật tồn tại trước đó, đồng thời hệ số xác định hiệu chỉnh (Adjusted R²) cũng tăng thêm 0.026, cho thấy mức độ giải thích của mô hình được cải thiện.

Bước tiếp theo là tiến hành kiểm định lại các giả thiết OLS để đảm bảo mô hình không còn sai sót và đạt độ tin cậy cao.

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

bgtest(OLS1)

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

Nhận xét

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 có có hiện tượng tư tương quan.

6.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

Nhận xét

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.

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

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

7 Dự báo

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(OLS1, newdata = newdata, interval = "prediction")

##     fit  lwr   upr
## 1 19265 6290 32241
## 2 18349 1161 35536

Nhận xét

Với X₂ = 1,100,000, X₃ = 5, X₄ = 4, mô hình ước lượng giá trị trung bình của Y là 19,265, với khoảng tin cậy 95% nằm trong khoảng [6,290 ; 32,241].
Với X₂ = 1,090,246, X₃ = 6, X₄ = 2.8, giá trị trung bình dự báo của Y là 18,349, và khoảng tin cậy 95% tương ứng là [1,161 ; 35,536].

7.1 Dự báo cho toàn bộ mẫu dữ liệu hiện có

library(Metrics)

## Warning: package 'Metrics' was built under R version 4.3.3

# Dự báo cho toàn bộ mẫu hiện có
data$Y_hat <- predict(OLS1)

# Hiển thị vài dòng đầu tiên để so sánh
head(data[, c("Y", "Y_hat")])

## # A tibble: 6 × 2
##       Y  Y_hat
##   <dbl>  <dbl>
## 1  2523 3653. 
## 2  2299   36.5
## 3  2641 2063. 
## 4  2677 5241. 
## 5  3706 4443. 
## 6  5340 4882.

7.2 biểu đồ so sánh giá trị thực và dự báo

library(ggplot2)

ggplot(data, aes(x = Y, y = Y_hat)) +
  geom_point(color = "blue") +
  geom_abline(slope = 1, intercept = 0, color = "red", linetype = "dashed") +
  labs(title = "So sánh giá trị thực và giá trị dự báo",
       x = "Giá trị thực tế (Y)",
       y = "Giá trị dự báo (Y_hat)") +
  theme_minimal()

Nhận xét

Biểu đồ cho thấy các điểm dữ liệu phân bố khá gần với đường chéo \(Y = \hat{Y}\), chứng tỏ mô hình OLS1 dự báo tương đối chính xác.

Một số điểm nằm lệch khỏi đường chéo, chủ yếu ở vùng giá trị nhỏ và lớn của Y, cho thấy mô hình có xu hướng đánh giá thấp hoặc cao nhẹ trong một số trường hợp.

Tuy nhiên, nhìn chung, mối quan hệ tuyến tính giữa giá trị thực và dự báo được thể hiện rõ, khẳng định mô hình có khả năng dự báo tốt.

THỰC HÀNH 1 NGÀY 06/10/2025

Lâm Thảo My

2025-10

1 Thống kê mô tả dữ liệu

2 Mô hình pooled OLS

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

4 Ma trận tương quan giữa các biến

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

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

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

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

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

6 Khắc phục khuyết tật mô hình

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

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

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

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

7 Dự báo

7.1 Dự báo cho toàn bộ mẫu dữ liệu hiện có

7.2 biểu đồ so sánh giá trị thực và dự báo