── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ plm::between() masks dplyr::between()
✖ dplyr::filter() masks stats::filter()
✖ plm::lag() masks dplyr::lag(), stats::lag()
✖ plm::lead() masks dplyr::lead()
✖ car::recode() masks dplyr::recode()
✖ car::some() masks purrr::some()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(margins)
Tạo biến
file_path <-"~/Downloads/PHÂN TÍCH ĐỊNH LƯỢNG/DATA MIDTERM PTĐL.xlsx"sheet_names <-excel_sheets(file_path)data <-map_dfr(sheet_names, function(sheet) { df <-read_excel(file_path, sheet = sheet)# Thêm năm từ tên sheet df$year <-as.numeric(sheet)return(df)})
# Phân bố GRDP (2010-2024)ggplot(datafinal, aes(x = grdp)) +geom_histogram(aes(y = ..density..), bins =25, fill ="#C57CAC", alpha =0.8) +geom_density(color ="#780062", linewidth =1) +facet_wrap(~year, scales ="free") +labs(title ="Phân bố GRDP (2010–2024)", x ="GRDP", y ="Mật độ xác suất") +theme_minimal()
Warning: The dot-dot notation (`..density..`) was deprecated in ggplot2 3.4.0.
ℹ Please use `after_stat(density)` instead.
-> Phân bố này lệch phải mạnh, phần lớn các tỉnh/thành phố có GRDP thấp (tập trung gần 0), trong khi có một số ít tỉnh/thành phố có GRDP rất cao (phần đuôi dài kéo sang phải), tạo ra các giá trị ngoại lai (outliers).
Phân bố ln_GRDP (2010-2024)
# Phân bố ln(GRDP) (2010-2024)ggplot(data, aes(x = ln_grdp)) +geom_histogram(aes(y = ..density..), bins =25, fill ="steelblue3", alpha =0.8) +geom_density(color ="darkblue", size =1) +facet_wrap(~year, scales ="free") +labs(title ="Phân bố ln(GRDP) (2010–2024)", x ="ln(GRDP)", y ="Mật độ xác suất") +theme_minimal()
Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
ℹ Please use `linewidth` instead.
-> phân phối của ln(GRDP) đã trở nên đối xứng hơn và gần với phân phối Chuẩn hơn hẳn.
-> Sử dụng ln_GRDP thay cho GRDP
Phân bố inc (2010-2024)
# Phân bố INC (2010-2024)ggplot(data, aes(x = inc)) +geom_histogram(aes(y = ..density..), bins =25, fill ="#FFA07A", alpha =0.8) +geom_density(color ="#CD5C5C", size =1) +facet_wrap(~year, scales ="free") +labs(title ="Phân bố INC (2010–2024)", x ="INC", y ="Mật độ xác suất") +theme_minimal()
-> phân phối của INC đều bị lệch phải rõ rệt, dữ liệu có phương sai lớn, và sự khác biệt tuyệt đối về thu nhập giữa các tỉnh giàu và nghèo là rất lớn
Phân bố ln_inc (2010-2024)
# Phân bố ln(INC) (2010-2024)ggplot(data, aes(x = ln_inc)) +geom_histogram(aes(y = ..density..), bins =25, fill ="#A2CD5A", alpha =0.8) +geom_density(color ="#458B00", size =1) +facet_wrap(~year, scales ="free") +labs(title ="Phân bố ln(INC) (2010–2024)", x ="ln(INC)", y ="Mật độ xác suất") +theme_minimal()
-> Sau khi lấy logarit, phân phối của ln_inc đã trở nên đối xứng hơn (hình chuông) và tập trung hơn quanh một giá trị trung bình. Phương sai ổn định, giảm rủi ro lỗi mô hình.
Phân bố lab (2010-2024)
# Phân bố LAB (2010-2024)ggplot(data, aes(x = lab)) +geom_histogram(aes(y = ..density..), bins =25, fill ="#BCD2EE", alpha =0.8) +geom_density(color ="#6E7B8B", size =1) +facet_wrap(~year, scales ="free") +labs(title ="Phân bố LAB (2010–2024)", x ="LAB", y ="Mật độ xác suất") +theme_minimal()
-> phân phối của LAB đều bị lệch mạnh sang phải. Đỉnh của biểu đồ (Mode) luôn nằm ở các giá trị rất thấp (thường là khoảng 10–20).
-> không phản ánh đúng giá trị phổ biến nhất. Điều này dễ vi phạm giả định về phân phối chuẩn của các phần dư trong mô hình hồi quy tuyến tính (OLS).
Phân bố ln_lab (2010-2024)
# Phân bố ln(lab)) (2010-2024)ggplot(data, aes(x = ln_lab)) +geom_histogram(aes(y = ..density..), bins =25, fill ="#00AE72", alpha =0.8) +geom_density(color ="#006241", size =1) +facet_wrap(~year, scales ="free") +labs(title ="Phân bố ln(LAB) (2010–2024)", x ="ln(LAB)", y ="Mật độ xác suất") +theme_minimal()
-> Sau khi logarit hóa, phân phối của ln_lab đã chuyển từ dạng lệch phải sang dạng tương đối đối xứng (hình chuông) và tập trung hơn.
-> Logarit hoá lab làm ổn định phương sai (giảm Heteroskedasticity) và giảm thiểu ảnh hưởng của các giá trị ngoại lai.
Mô tả thống kê
summary (datafinal)
province ln_grdp grdp gini
Length:693 Min. :0.4055 Min. : 1.5 Min. :0.2000
Class :character 1st Qu.:3.5496 1st Qu.: 34.8 1st Qu.:0.3300
Mode :character Median :3.9435 Median : 51.6 Median :0.3600
Mean :3.9322 Mean : 62.8 Mean :0.3613
3rd Qu.:4.2959 3rd Qu.: 73.4 3rd Qu.:0.3900
Max. :5.9525 Max. :384.7 Max. :0.5300
ginisq ln_lab educ ln_inc
Min. :0.0400 Min. :1.649 Min. :60.20 Min. :-0.5674
1st Qu.:0.1090 1st Qu.:2.632 1st Qu.:92.60 1st Qu.: 0.7237
Median :0.1296 Median :2.902 Median :95.40 Median : 1.1439
Mean :0.1328 Mean :2.913 Mean :93.52 Mean : 1.0480
3rd Qu.:0.1520 3rd Qu.:3.210 3rd Qu.:97.70 3rd Qu.: 1.4306
Max. :0.2809 Max. :3.987 Max. :99.45 Max. : 2.1902
povr year
Min. : 0.000 Min. :2010
1st Qu.: 2.400 1st Qu.:2014
Median : 6.000 Median :2019
Mean : 9.221 Mean :2018
3rd Qu.:12.800 3rd Qu.:2022
Max. :50.800 Max. :2024
# Chọn các biếnvars <-c("ln_grdp","gini","ginisq","ln_lab","educ","ln_inc","povr")# Tính ma trận tương quan, bỏ NAcor_mat <-cor(datafinal[, vars], use="complete.obs")# Chỉ in nửa dưới (lower triangle)cor_mat[upper.tri(cor_mat)] <-NAprint(cor_mat)
ln_grdp gini ginisq ln_lab educ ln_inc povr
ln_grdp 1.0000000 NA NA NA NA NA NA
gini -0.2141026 1.0000000 NA NA NA NA NA
ginisq -0.2063057 0.9941681 1.0000000 NA NA NA NA
ln_lab 0.6510270 -0.1359710 -0.1250014 1.0000000 NA NA NA
educ 0.4169517 -0.5426827 -0.5628749 0.3583830 1.0000000 NA NA
ln_inc 0.8600157 -0.3314782 -0.3210916 0.6002146 0.5326628 1.0000000 NA
povr -0.6532811 0.5169314 0.5295523 -0.3924475 -0.7889024 -0.8052968 1
-> gini và ginisq có tương quan rất cao (0.994), cảnh báo đa cộng tuyến nếu cùng đưa vào mô hình. Cần cân nhắc: Chỉ giữ gini hoặc ginisq
Ma trận tương quan có màu
library (corrplot)
corrplot 0.95 loaded
corrplot( cor_mat, method ="color", # tô màutype ="lower", # chỉ nửa dưới (bậc thang)tl.col ="#551A8B", # màu chữ tên biếntl.srt =45, # xoay tên biến 45 độaddCoef.col ="#FFBBFF", # hiển thị giá trị tương quannumber.cex =0.8, # cỡ chữ sốdiag =TRUE,col =colorRampPalette(c("#AB82FF", "#8A2BE2"))(200) )
Warning in ind1:ind2: numerical expression has 3 elements: only the first used
Ước lượng các mô hình
Chuyển sang dữ liệu bảng
pdata <-pdata.frame(datafinal, index =c("province", "year"))
VIF > 10 –> có đa cộng tuyến rất mạnh giữa gini và ginisq
Khắc phục đa cộng tuyến
# Chuẩn hóa Gini (z-score)datafinal$gini_c <-scale(datafinal$gini, center =TRUE, scale =TRUE)# Tạo bình phương của Gini chuẩn hóadatafinal$ginisq_c <- datafinal$gini_c^2
Chạy lại POOLED_OLS
pdata <-pdata.frame(datafinal, index =c("province", "year"))
ols_model <-plm(ln_grdp ~ gini_c + ginisq_c + ln_lab + educ + ln_inc +povr, data = pdata, model ="pooling")summary (ols_model)
Pooling Model
Call:
plm(formula = ln_grdp ~ gini_c + ginisq_c + ln_lab + educ + ln_inc +
povr, data = pdata, model = "pooling")
Balanced Panel: n = 63, T = 11, N = 693
Residuals:
Min. 1st Qu. Median 3rd Qu. Max.
-2.491162 -0.158408 -0.021663 0.142780 1.918783
Coefficients:
Estimate Std. Error t-value Pr(>|t|)
(Intercept) 2.6835242 0.3244929 8.2699 6.955e-16 ***
gini_c 0.0299779 0.0141965 2.1116 0.03508 *
ginisq_c -0.0158923 0.0078728 -2.0186 0.04391 *
ln_lab 0.3380472 0.0374925 9.0164 < 2.2e-16 ***
educ -0.0074873 0.0032528 -2.3018 0.02165 *
ln_inc 0.9522888 0.0503382 18.9178 < 2.2e-16 ***
povr -0.0019492 0.0033091 -0.5890 0.55603
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Total Sum of Squares: 284.66
Residual Sum of Squares: 64.145
R-Squared: 0.77466
Adj. R-Squared: 0.77269
F-statistic: 393.054 on 6 and 686 DF, p-value: < 2.22e-16
fem_ind <-plm(formula, data = pdata, model ="within", effect ="individual")fem_time <-plm(formula, data = pdata, model ="within", effect ="time")fem_tw <-plm(formula, data = pdata, model ="within", effect ="twoways")summary (fem_ind)
Oneway (individual) effect Within Model
Call:
plm(formula = formula, data = pdata, effect = "individual", model = "within")
Balanced Panel: n = 63, T = 11, N = 693
Residuals:
Min. 1st Qu. Median 3rd Qu. Max.
-1.8740991 -0.0727263 0.0032553 0.0773675 0.9731919
Coefficients:
Estimate Std. Error t-value Pr(>|t|)
gini_c -0.03578690 0.01305753 -2.7407 0.006306 **
ginisq_c 0.00594303 0.00658143 0.9030 0.366874
ln_lab 0.15039226 0.06494956 2.3155 0.020908 *
educ 0.01347172 0.00651760 2.0670 0.039148 *
ln_inc 0.92774882 0.04788535 19.3744 < 2.2e-16 ***
povr -0.00090084 0.00309647 -0.2909 0.771205
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Total Sum of Squares: 157.18
Residual Sum of Squares: 25.552
R-Squared: 0.83744
Adj. R-Squared: 0.81972
F-statistic: 535.753 on 6 and 624 DF, p-value: < 2.22e-16
summary (fem_time)
Oneway (time) effect Within Model
Call:
plm(formula = formula, data = pdata, effect = "time", model = "within")
Balanced Panel: n = 63, T = 11, N = 693
Residuals:
Min. 1st Qu. Median 3rd Qu. Max.
-2.521168 -0.154212 -0.014747 0.140847 1.737194
Coefficients:
Estimate Std. Error t-value Pr(>|t|)
gini_c 0.0605146 0.0152552 3.9668 8.061e-05 ***
ginisq_c -0.0131635 0.0081173 -1.6217 0.105342
ln_lab 0.3242944 0.0377726 8.5854 < 2.2e-16 ***
educ -0.0092089 0.0032616 -2.8234 0.004892 **
ln_inc 1.1295754 0.0822812 13.7282 < 2.2e-16 ***
povr -0.0008079 0.0037351 -0.2163 0.828820
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Total Sum of Squares: 151.02
Residual Sum of Squares: 61.341
R-Squared: 0.59381
Adj. R-Squared: 0.5842
F-statistic: 164.709 on 6 and 676 DF, p-value: < 2.22e-16
summary (fem_tw)
Twoways effects Within Model
Call:
plm(formula = formula, data = pdata, effect = "twoways", model = "within")
Balanced Panel: n = 63, T = 11, N = 693
Residuals:
Min. 1st Qu. Median 3rd Qu. Max.
-1.9462144 -0.0551185 -0.0032316 0.0610168 0.9530155
Coefficients:
Estimate Std. Error t-value Pr(>|t|)
gini_c -1.6630e-02 1.4112e-02 -1.1785 0.239056
ginisq_c 9.5002e-03 6.7423e-03 1.4090 0.159330
ln_lab -3.9942e-03 6.9063e-02 -0.0578 0.953899
educ -5.0281e-05 6.6828e-03 -0.0075 0.993999
ln_inc 3.2797e-01 1.0614e-01 3.0900 0.002092 **
povr -5.9627e-03 3.2141e-03 -1.8552 0.064049 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Total Sum of Squares: 23.535
Residual Sum of Squares: 22.636
R-Squared: 0.038199
Adj. R-Squared: -0.083985
F-statistic: 4.06423 on 6 and 614 DF, p-value: 0.00052422
Kiểm định Ftest xem mô hình FEM nào là lựa chọn tối ưu
# Kiểm định F để chọn FEM individual vs OLSpFtest(fem_ind, ols_model) # Individual
F test for individual effects
data: formula
F = 15.201, df1 = 62, df2 = 624, p-value < 2.2e-16
alternative hypothesis: significant effects
pFtest(fem_time, ols_model) # Time
F test for time effects
data: formula
F = 3.0904, df1 = 10, df2 = 676, p-value = 0.0007376
alternative hypothesis: significant effects
pFtest(fem_tw, ols_model) # Twoways
F test for twoways effects
data: formula
F = 15.639, df1 = 72, df2 = 614, p-value < 2.2e-16
alternative hypothesis: significant effects
Individual effect: F = 15.32, p < 0.05 -> Có hiệu ứng cố định theo tỉnh
Time effect: F = 3.0756, p < 0.05 -> Có hiệu ứng cố định theo năm
Twoways: F = 15.668, p < 0.01 -> Có hiệu ứng cố định theo cả tỉnh và năm
Mặc dù FEM twoway là phù hợp nhưng kết quả cho thấy R2 quá thấp, nhiều biến mất ý nghĩa thống kê. Trong khi đó FEM individual giải thích tốt sự biến thiên của ln_GRDP giữa các địa phương, R2 cao và các biến đều có ý nghĩa thống kê.
– Chọn mô hình FEM individual.
REM
rem_ind <-plm(formula, data = pdata, model ="random", effect ="individual")rem_time <-plm(formula, data = pdata, model ="random", effect ="time")rem_tw <-plm(formula, data = pdata, model ="random", effect ="twoways")summary (rem_ind)
Oneway (individual) effect Random Effect Model
(Swamy-Arora's transformation)
Call:
plm(formula = formula, data = pdata, effect = "individual", model = "random")
Balanced Panel: n = 63, T = 11, N = 693
Effects:
var std.dev share
idiosyncratic 0.04095 0.20236 0.468
individual 0.04651 0.21567 0.532
theta: 0.7278
Residuals:
Min. 1st Qu. Median 3rd Qu. Max.
-2.0780444 -0.0832551 -0.0019043 0.0817365 1.1273270
Coefficients:
Estimate Std. Error z-value Pr(>|z|)
(Intercept) 2.19037396 0.44080275 4.9691 6.728e-07 ***
gini_c -0.02101911 0.01249479 -1.6822 0.09252 .
ginisq_c 0.00514770 0.00649879 0.7921 0.42830
ln_lab 0.22788022 0.05384936 4.2318 2.318e-05 ***
educ 0.00093211 0.00444969 0.2095 0.83408
ln_inc 0.93441271 0.04404207 21.2164 < 2.2e-16 ***
povr 0.00070018 0.00294739 0.2376 0.81222
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Total Sum of Squares: 166.63
Residual Sum of Squares: 28.978
R-Squared: 0.82609
Adj. R-Squared: 0.82457
Chisq: 3258.56 on 6 DF, p-value: < 2.22e-16
summary (rem_time)
Oneway (time) effect Random Effect Model
(Swamy-Arora's transformation)
Call:
plm(formula = formula, data = pdata, effect = "time", model = "random")
Balanced Panel: n = 63, T = 11, N = 693
Effects:
var std.dev share
idiosyncratic 0.09074 0.30123 1
time 0.00000 0.00000 0
theta: 0
Residuals:
Min. 1st Qu. Median 3rd Qu. Max.
-2.491162 -0.158408 -0.021663 0.142780 1.918783
Coefficients:
Estimate Std. Error z-value Pr(>|z|)
(Intercept) 2.6835242 0.3244929 8.2699 < 2e-16 ***
gini_c 0.0299779 0.0141965 2.1116 0.03472 *
ginisq_c -0.0158923 0.0078728 -2.0186 0.04352 *
ln_lab 0.3380472 0.0374925 9.0164 < 2e-16 ***
educ -0.0074873 0.0032528 -2.3018 0.02135 *
ln_inc 0.9522888 0.0503382 18.9178 < 2e-16 ***
povr -0.0019492 0.0033091 -0.5890 0.55583
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Total Sum of Squares: 284.66
Residual Sum of Squares: 64.145
R-Squared: 0.77466
Adj. R-Squared: 0.77269
Chisq: 2358.32 on 6 DF, p-value: < 2.22e-16
summary (rem_tw)
Twoways effects Random Effect Model
(Swamy-Arora's transformation)
Call:
plm(formula = formula, data = pdata, effect = "twoways", model = "random")
Balanced Panel: n = 63, T = 11, N = 693
Effects:
var std.dev share
idiosyncratic 0.03687 0.19200 0.44
individual 0.04688 0.21652 0.56
time 0.00000 0.00000 0.00
theta: 0.7417 (id) 0 (time) 0 (total)
Residuals:
Min. 1st Qu. Median 3rd Qu. Max.
-2.0687105 -0.0812921 -0.0034819 0.0797874 1.1111791
Coefficients:
Estimate Std. Error z-value Pr(>|z|)
(Intercept) 2.16143359 0.44758983 4.8290 1.372e-06 ***
gini_c -0.02191886 0.01247834 -1.7566 0.07899 .
ginisq_c 0.00535612 0.00647822 0.8268 0.40836
ln_lab 0.22326331 0.05443736 4.1013 4.109e-05 ***
educ 0.00138263 0.00452583 0.3055 0.75999
ln_inc 0.93475676 0.04410035 21.1961 < 2.2e-16 ***
povr 0.00066645 0.00294508 0.2263 0.82097
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Total Sum of Squares: 165.69
Residual Sum of Squares: 28.662
R-Squared: 0.82701
Adj. R-Squared: 0.8255
Chisq: 3279.59 on 6 DF, p-value: < 2.22e-16
# Tạo một custom model object để đưa vào modelsummaryfem_ind_DK <- fem_indfem_ind_DK$vcov <- vcov_DK_ind# Gộp các modelModels <-list("POOLED"= ols_model,"FEM_ind"= fem_ind,"REM_ind"= rem_ind,"FEM_ind_DK"= fem_ind_DK)# Hiển thị bảngmodelsummary::modelsummary( Models,statistic ="{estimate} ({std.error})",stars =c("*"=0.1, "**"=0.05, "***"=0.01))
POOLED
FEM_ind
REM_ind
FEM_ind_DK
* p < 0.1, ** p < 0.05, *** p < 0.01
(Intercept)
2.684***
2.190***
2.684 (0.324)
2.190 (0.441)
gini_c
0.030**
-0.036***
-0.021*
-0.036**
0.030 (0.014)
-0.036 (0.013)
-0.021 (0.012)
-0.036 (0.015)
ginisq_c
-0.016**
0.006
0.005
0.006
-0.016 (0.008)
0.006 (0.007)
0.005 (0.006)
0.006 (0.004)
ln_lab
0.338***
0.150**
0.228***
0.150***
0.338 (0.037)
0.150 (0.065)
0.228 (0.054)
0.150 (0.027)
educ
-0.007**
0.013**
0.001
0.013***
-0.007 (0.003)
0.013 (0.007)
0.001 (0.004)
0.013 (0.005)
ln_inc
0.952***
0.928***
0.934***
0.928***
0.952 (0.050)
0.928 (0.048)
0.934 (0.044)
0.928 (0.021)
povr
-0.002
-0.001
0.001
-0.001
-0.002 (0.003)
-0.001 (0.003)
0.001 (0.003)
-0.001 (0.001)
Num.Obs.
693
693
693
693
R2
0.775
0.837
0.826
0.837
R2 Adj.
0.773
0.820
0.825
0.820
AIC
333.4
-306.5
-217.3
-306.5
BIC
369.7
-274.7
-180.9
-274.7
RMSE
0.30
0.19
0.20
0.19
Tính turning point (ngưỡng bất bình đẳng)
Đạo hàm theo gini để tìm gini cực trị
H=B1 + 2.B2
hay GINI* = - B1 / 2.B2
- Khi H < 0 hay GINI* > turning point, gia tăng bất bình đẳng làm giảm tăng trưởng kinh tế
- Khi H > 0 hay GINI* < turning point, gia tăng bất bình đẳng làm tăng trưởng kinh tế
GINI* : turning point trên thang đo gốc
GINI_c* : turning point trên biến chuẩn hoá
Chọn mô hình POOLED để tính ngưỡng bất bình đẳng hay GINI* vì chỉ có POOLED thì các biến gini_c và ginisq_c mới có ý nghĩa thống kê cao ở CI 95%.
-> Nhận thấy GINI_c* ~ 0.94 quá cao, cần đưa về GINI* của dữ liệu gốc (tức chưa chuẩn hoá gini và ginisq), nhóm sử dụng công thức sau để tính GINI_c*:
GINI_c* = B gini_c / 2.B ginisq_c, sau đó chuẩn hoá biến bằng công thức:
#Lấy hệ sốb1 <-coef(ols_model)["gini_c"]b2 <-coef(ols_model)["ginisq_c"]#Tính turning point trên biến chuẩn hoá (gini_c scale)gini_c_star <--b1 / (2* b2)gini_c_star
gini_c
0.9431567
#Chuyển turning point về gini gốc ---gini_star <- gini_c_star * sd_gini + mean_ginigini_star
gini_c
0.4059399
# Chia nhóm theo Ginidata$gini_group <-ifelse(data$gini <0.40, "Gini < 0.40", "Gini >= 0.40")# Tính số lượng tỉnh/thành theo nhóm mỗi nămlibrary(dplyr)count_table <- data %>%group_by(year, gini_group) %>%summarise(num_provinces =n(), .groups ="drop")print(count_table)
# Tính GRDP trung bình theo nhómmean_grdp <- data %>%group_by(gini_group) %>%summarise(mean_gini =mean(gini),mean_grdp_growth =mean(ln_grdp),.groups ="drop")print(mean_grdp)
