library(haven)
D201 <- haven::read_dta("C:\\Users\\DELL\\Downloads\\2 - Data\\VHLSS2020_Household_Data\\MUC8.dta")
D202 <- haven::read_dta("C:\\Users\\DELL\\Downloads\\2 - Data\\VHLSS2020_Household_Data\\HO1.dta")
D203 <- haven::read_dta("C:\\Users\\DELL\\Downloads\\2 - Data\\VHLSS2020_Household_Data\\MUC1A.dta")
D204 <- haven::read_dta("C:\\Users\\DELL\\Downloads\\2 - Data\\VHLSS2020_Household_Data\\MUC2X.dta")
D205 <- haven::read_dta("C:\\Users\\DELL\\Downloads\\2 - Data\\VHLSS2020_Household_Data\\MUC2V.dta")
D206 <- haven::read_dta("C:\\Users\\DELL\\Downloads\\2 - Data\\VHLSS2020_Household_Data\\MUC4A.dta")
D207 <- haven::read_dta("C:\\Users\\DELL\\Downloads\\2 - Data\\VHLSS2020_Household_Data\\MUC82.dta")Cuối kỳ
Tải dữ liệu
Các biến
Biến phụ thuộc: m8c120 - Hộ nghèo 2020 (D201)
library(tidyverse)── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr 1.1.4 ✔ readr 2.1.5
✔ forcats 1.0.0 ✔ stringr 1.5.1
✔ ggplot2 3.5.1 ✔ tibble 3.2.1
✔ lubridate 1.9.3 ✔ tidyr 1.3.1
✔ purrr 1.0.2
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag() masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errorsD201 <- D201 %>%
rename( `tình trạng nghèo` = m8c120)D201$`tình trạng nghèo` <- ifelse(is.na(D201$`tình trạng nghèo`), NA,
ifelse(D201$`tình trạng nghèo` %in% c(2), 0, 1))Biến độc lập
- Khu vực cư trú: ttnt (D202)
D202$ttnt <- ifelse(is.na(D202$ttnt), NA,
ifelse(D202$ttnt %in% c(2), 1, 0))- Dân tộc chủ hộ: dantoc (D202)
D202$dantoc <- ifelse(D202$dantoc %in% c(1,4), 1, 0)- Giới tính chủ hộ (giới tính): biến m1ac2 với mã thành viên (m1ama) bằng 1 (D203)
Data <- D203 %>%
filter(m1ama == 1)Data <- Data %>%
rename( `giới tính` = m1ac2)Data$`giới tính` <- ifelse(is.na(Data$`giới tính`), NA,
ifelse(Data$`giới tính` %in% c(2), 0, 1))- Tuổi của chủ hộ (tuổi): m1ac5 (Data)
Data <- Data %>%
rename( `tuổi` = m1ac5)Quy mô hộ gia đình: tsnguoi (D202)
Trình độ học vấn của chủ hộ (trình độ học vấn): m2xc1 (D204) và m2vc1 (D205) với mã thành viên (m2xma và m2vma) bằng 1
Data1 <- merge(D204, D205, by = c("tinh", "huyen", "xa", "diaban", "hoso"), all = TRUE)Data1 <- Data1 %>%
filter(m2xma == 1 | m2vma == 1)Data1 <- Data1 %>%
mutate(`trình độ học vấn` = coalesce(m2vc1, m2xc1))table(Data1$`trình độ học vấn`)
0 00 01 02 03 04 05 06 07 08 09 1 10
247 6 24 65 66 64 150 73 69 38 390 640 1190
11 12 2 2. 3 4 4. 5 6 7 8 9 99
807 12625 1674 1 2158 1807 1 4786 2496 2268 1543 11401 2451 Data1 <- Data1 %>%
mutate(`trình độ học vấn` = case_when(
`trình độ học vấn` %in% c("0", "00","99") ~ "0",
`trình độ học vấn`%in% c("01", "1") ~ "1",
`trình độ học vấn` %in% c("02", "2") ~ "2",
`trình độ học vấn` %in% c("03", "3") ~ "3",
`trình độ học vấn` %in% c("04", "4") ~ "4",
`trình độ học vấn` %in% c("05", "5") ~ "5",
`trình độ học vấn` %in% c("06", "6") ~ "6",
`trình độ học vấn` %in% c("07", "7") ~ "7",
`trình độ học vấn` %in% c("08", "8") ~ "8",
`trình độ học vấn` %in% c("09", "9") ~ "9",
`trình độ học vấn` %in% c("10", "11", "12") ~ as.character(`trình độ học vấn`),
TRUE ~ NA_character_))table(Data1$`trình độ học vấn`)
0 1 10 11 12 2 3 4 5 6 7 8 9
2704 664 1190 807 12625 1739 2224 1871 4936 2569 2337 1581 11791 Data1$`trình độ học vấn` <- as.numeric(as.character(Data1$`trình độ học vấn`))- Tiếp cận tín dụng (tin_dung)
D207$` nguồn vay` <- ifelse(D207$m8c7 %in% 1:10, 1, 0)D207 <- D207 %>% filter(m8ma != 0)D207 <- D207 %>%
group_by(tinh, huyen, xa, diaban, hoso) %>%
mutate(tin_dung = ifelse(any(` nguồn vay` == 1), 1, 0)) %>%
ungroup()D207 <- D207 %>%
group_by(tinh, huyen, xa, diaban, hoso) %>%
slice(1) %>%
ungroup()- Việc làm phi nông nghiệp của chủ hộ (phi NN): m4ac1c (D206) với mã thành viên (m4ama) bằng 1
Data2 <- D206 %>%
filter(m4ama == 1)Data2 <- Data2 %>%
rename( `phi NN` = m4ac1c)Data2$`phi NN` <- ifelse(is.na(Data2$`phi NN`), NA,
ifelse(Data2$`phi NN` %in% c(2), 0, 1))- Tỷ lệ phụ thuộc (tỷ lệ phụ thuộc) = (số người phụ thuộc/tổng số người)*100%
Số người phụ thuộc: nam dưới 15 và trên 59 tuổi, nữ dưới 15 và trên 54 tuổi
Data3 <- D203 %>%
filter(!is.na(m1ac5) & !is.na(m1ac2)) %>%
mutate(
dependent = case_when(
(m1ac2 == 1 & (m1ac5 < 15 | m1ac5 > 59)) ~ 1,
(m1ac2 == 2 & (m1ac5 < 15 | m1ac5 > 54)) ~ 1,
TRUE ~ 0
)
) %>%
group_by(tinh, huyen, xa, diaban, hoso) %>%
summarise(total_dependents = sum(dependent))`summarise()` has grouped output by 'tinh', 'huyen', 'xa', 'diaban'. You can
override using the `.groups` argument.Data3$ID_match <- paste(Data3$tinh, Data3$huyen, Data3$xa, Data3$diaban, Data3$hoso, sep="_")
D202$ID_match <- paste(D202$tinh, D202$huyen, D202$xa, D202$diaban, D202$hoso, sep="_")
Data3$`tỷ lệ phụ thuộc` <- Data3$total_dependents / D202$tsnguoi[match(Data3$ID_match, D202$ID_match)]Data3$`tỷ lệ phụ thuộc` <- Data3$`tỷ lệ phụ thuộc`*100Ghép dữ liệu
DATA <- D201 %>%
full_join(D202, by = c("tinh", "huyen", "xa", "diaban", "hoso")) %>%
full_join(Data, by= c("tinh", "huyen", "xa", "diaban", "hoso")) %>%
full_join(Data1, by = c("tinh", "huyen", "xa", "diaban", "hoso")) %>%
full_join(Data2, by = c("tinh", "huyen", "xa", "diaban", "hoso")) %>%
full_join(Data3, by = c("tinh", "huyen", "xa", "diaban", "hoso")) %>%
full_join(D207, by = c("tinh", "huyen", "xa", "diaban", "hoso"))DATA <- DATA %>%
select(tinh, huyen, xa, diaban, hoso, `tình trạng nghèo`, `ttnt`, dantoc, tsnguoi, `giới tính`, `phi NN`, `tuổi`,`trình độ học vấn`, tin_dung, `tỷ lệ phụ thuộc`)DATA$tinh <- as_factor(DATA$tinh)levels(DATA$tinh) [1] "Thành phố Hà Nội" "Tỉnh Hà Giang" "Tỉnh Cao Bằng"
[4] "Tỉnh Bắc Kạn" "Tỉnh Tuyên Quang" "Tỉnh Lào Cai"
[7] "Tỉnh Điện Biên" "Tỉnh Lai Châu" "Tỉnh Sơn La"
[10] "Tỉnh Yên Bái" "Tỉnh Hoà Bình" "Tỉnh Thái Nguyên"
[13] "Tỉnh Lạng Sơn" "Tỉnh Quảng Ninh" "Tỉnh Bắc Giang"
[16] "Tỉnh Phú Thọ" "Tỉnh Vĩnh Phúc" "Tỉnh Bắc Ninh"
[19] "Tỉnh Hải Dương" "Thành phố Hải Phòng" "Tỉnh Hưng Yên"
[22] "Tỉnh Thái Bình" "Tỉnh Hà Nam" "Tỉnh Nam Định"
[25] "Tỉnh Ninh Bình" "Tỉnh Thanh Hoá" "Tỉnh Nghệ An"
[28] "Tỉnh Hà Tĩnh" "Tỉnh Quảng Bình" "Tỉnh Quảng Trị"
[31] "Tỉnh Thừa Thiên Huế" "Thành phố Đà Nẵng" "Tỉnh Quảng Nam"
[34] "Tỉnh Quảng Ngãi" "Tỉnh Bình Định" "Tỉnh Phú Yên"
[37] "Tỉnh Khánh Hoà" "Tỉnh Ninh Thuận" "Tỉnh Bình Thuận"
[40] "Tỉnh Kon Tum" "Tỉnh Gia Lai" "Tỉnh Đắk Lắk"
[43] "Tỉnh Đắk Nông" "Tỉnh Lâm Đồng" "Tỉnh Bình Phước"
[46] "Tỉnh Tây Ninh" "Tỉnh Bình Dương" "Tỉnh Đồng Nai"
[49] "Tỉnh Bà Rịa - Vũng Tàu" "Thành phố Hồ Chí Minh" "Tỉnh Long An"
[52] "Tỉnh Tiền Giang" "Tỉnh Bến Tre" "Tỉnh Trà Vinh"
[55] "Tỉnh Vĩnh Long" "Tỉnh Đồng Tháp" "Tỉnh An Giang"
[58] "Tỉnh Kiên Giang" "Thành phố Cần Thơ" "Tỉnh Hậu Giang"
[61] "Tỉnh Sóc Trăng" "Tỉnh Bạc Liêu" "Tỉnh Cà Mau" DATA$tinh <- as.numeric(DATA$tinh)DATA <- DATA %>%
mutate(TD_MNPB = ifelse(tinh %in% c(2:13, 15, 16), 1, 0),
BTB_DHMT = ifelse(tinh %in% c(26:39), 1, 0),
TN = ifelse(tinh %in% c(40:44), 1, 0),
DBSCL = ifelse(tinh %in% c(51:63), 1, 0),
DNB = ifelse(tinh %in% c(45:50), 1, 0))Thống kê mô tả
Bảng tỷ lệ hộ nghèo
DATA %>%
select(`tình trạng nghèo`) %>%
table() %>%
prop.table() %>%
round(digits = 3)*100tình trạng nghèo
0 1
94.6 5.4 Bảng tỷ lệ hộ nghèo theo dân tộc
DATA%>%
select(`tình trạng nghèo`, dantoc) %>%
table() %>%
prop.table(1) %>%
round(digits = 4)*100 dantoc
tình trạng nghèo 0 1
0 13.82 86.18
1 62.56 37.44Bảng tỷ lệ hộ nghèo theo khu vực
DATA%>%
select(`tình trạng nghèo`, ttnt) %>%
table() %>%
prop.table(1) %>%
round(digits = 4)*100 ttnt
tình trạng nghèo 0 1
0 34.00 66.00
1 12.35 87.65Bảng tỷ lệ hộ nghèo theo giới tính chủ hộ
DATA%>%
select(`tình trạng nghèo`, `giới tính`) %>%
table() %>%
prop.table(1) %>%
round(digits = 4)*100 giới tính
tình trạng nghèo 0 1
0 25.90 74.10
1 32.73 67.27Bảng thống kê định lượng theo tình trạng nghèo
DATA %>%
group_by(`tình trạng nghèo`) %>%
summarise(
TuoiTB = mean(`tuổi`),
QuymohoTB = mean(tsnguoi),
TrinhdohocvanTB = mean(`trình độ học vấn`),
TylephuthuocTB = mean(`tỷ lệ phụ thuộc`))# A tibble: 3 × 5
`tình trạng nghèo` TuoiTB QuymohoTB TrinhdohocvanTB TylephuthuocTB
<dbl> <dbl> <dbl> <dbl> <dbl>
1 0 51.1 3.67 NA 43.0
2 1 49.0 3.86 4.39 50.9
3 NA NA NA 8.65 NA Mô hình
Khai báo biến giả
DATA <- DATA %>%
mutate(
`tình trạng nghèo` = as.factor(`tình trạng nghèo`),
`khu vực` = as.factor(ttnt),
TD_MNPB= as.factor(TD_MNPB),
BTB_DHMT = as.factor(BTB_DHMT),
TN = as.factor(TN),
DBSCL = as.factor(DBSCL),
DNB = as.factor(DNB),
`dân tộc` = as.factor(dantoc),
tin_dung = as.factor(tin_dung),
`giới tính` = as.factor(`giới tính`),
`phi NN` = as.factor(`phi NN`))Mô hình logit
logit <- glm(`tình trạng nghèo` ~ `khu vực`+`dân tộc` + `tuổi` + `giới tính` + `tsnguoi` + `trình độ học vấn` + tin_dung + `tỷ lệ phụ thuộc` + `phi NN` + TD_MNPB + TN + DNB + DBSCL + BTB_DHMT, family = binomial(link = "logit"), data = DATA)Mô hình probit
probit <- glm( `tình trạng nghèo` ~ `khu vực`+`dân tộc` + `tuổi` + `giới tính` + `tsnguoi` + `trình độ học vấn` + tin_dung + `tỷ lệ phụ thuộc` + `phi NN` + TD_MNPB + TN + DNB + DBSCL + BTB_DHMT, family = binomial(link = "probit"), data = DATA)So sánh mô hình
Models <- list(logit, probit)
modelsummary::modelsummary(Models,
statistic = c('{statistic} '),
stars = c("*"=0.1, "**"=0.05, "***"=0.01))| (1) | (2) | |
|---|---|---|
| * p < 0.1, ** p < 0.05, *** p < 0.01 | ||
| (Intercept) | -0.858** | -0.405** |
| -2.430 | -2.363 | |
| khu vực1 | 0.370*** | 0.165*** |
| 3.003 | 2.786 | |
| dân tộc1 | -1.106*** | -0.582*** |
| -11.347 | -11.700 | |
| tuổi | -0.032*** | -0.015*** |
| -9.793 | -9.159 | |
| giới tính1 | -0.492*** | -0.261*** |
| -5.340 | -5.449 | |
| tsnguoi | -0.009 | -0.009 |
| -0.365 | -0.686 | |
| trình độ học vấn | -0.181*** | -0.097*** |
| -17.185 | -17.399 | |
| tin_dung1 | 0.894*** | 0.440*** |
| 5.183 | 5.323 | |
| tỷ lệ phụ thuộc | 0.015*** | 0.007*** |
| 9.484 | 9.098 | |
| phi NN1 | -0.836*** | -0.384*** |
| -6.256 | -6.152 | |
| TD_MNPB1 | 1.050*** | 0.472*** |
| 4.949 | 4.967 | |
| TN1 | 0.519** | 0.208** |
| 2.329 | 2.062 | |
| DNB1 | -0.857** | -0.451*** |
| -2.518 | -3.032 | |
| DBSCL1 | -0.099 | -0.126 |
| -0.443 | -1.280 | |
| BTB_DHMT1 | 0.422** | 0.135 |
| 1.988 | 1.445 | |
| Num.Obs. | 11175 | 11175 |
| AIC | 5278.5 | 5300.6 |
| BIC | 5388.3 | 5410.4 |
| Log.Lik. | -2624.228 | -2635.279 |
| F | 85.188 | 87.922 |
| RMSE | 0.26 | 0.26 |
Mô hình logit có AIC, BIC nhỏ hơn => chọn mô hình logit
Các kiểm định
Đa cộng tuyến
Kiểm tra đa cộng tuyến với các biến định lượng
logit1 <- glm(`tình trạng nghèo` ~ `tuổi` + `tsnguoi` +
`trình độ học vấn` + `tỷ lệ phụ thuộc`,
family = binomial(link = "logit"), data = DATA)car::vif(logit1) tuổi tsnguoi `trình độ học vấn` `tỷ lệ phụ thuộc`
1.439082 1.083404 1.091320 1.302598 VIF <5 => không có hiện tượng đa cộng tuyến
Kiểm định The Hosmer-Lemeshow GOF
performance::performance_hosmer(logit, n_bins=15)# Hosmer-Lemeshow Goodness-of-Fit Test
Chi-squared: 19.133
df: 13
p-value: 0.119Summary: model seems to fit well.p-value >= 0.05: Mô hình phù hợp tốt với dữ liệu
Chỉ số OR
cbind(Estimate=round(coef(logit),3),
OR=round(exp(coef(logit)),3)) Estimate OR
(Intercept) -0.858 0.424
`khu vực`1 0.370 1.447
`dân tộc`1 -1.106 0.331
tuổi -0.032 0.969
`giới tính`1 -0.492 0.611
tsnguoi -0.009 0.991
`trình độ học vấn` -0.181 0.835
tin_dung1 0.894 2.445
`tỷ lệ phụ thuộc` 0.015 1.015
`phi NN`1 -0.836 0.433
TD_MNPB1 1.050 2.857
TN1 0.519 1.681
DNB1 -0.857 0.425
DBSCL1 -0.099 0.906
BTB_DHMT1 0.422 1.525Đánh giá ảnh hưởng biên
mfx::logitmfx( `tình trạng nghèo`~ `khu vực`+`dân tộc` + `tuổi` + `giới tính` + `tsnguoi` + `trình độ học vấn` + tin_dung + `tỷ lệ phụ thuộc` + `phi NN` + TD_MNPB + TN + DNB + DBSCL + BTB_DHMT, data = DATA)Call:
mfx::logitmfx(formula = `tình trạng nghèo` ~ `khu vực` +
`dân tộc` + tuổi + `giới tính` + tsnguoi + `trình độ học vấn` +
tin_dung + `tỷ lệ phụ thuộc` + `phi NN` + TD_MNPB +
TN + DNB + DBSCL + BTB_DHMT, data = DATA)
Marginal Effects:
dF/dx Std. Err. z P>|z|
`khu vực`1 1.4144e-02 4.2571e-03 3.3225 0.0008922 ***
`dân tộc`1 -6.0557e-02 6.9089e-03 -8.7651 < 2.2e-16 ***
tuổi -1.3328e-03 1.4064e-04 -9.4767 < 2.2e-16 ***
`giới tính`1 -2.3760e-02 5.0601e-03 -4.6955 2.659e-06 ***
tsnguoi -3.5922e-04 9.8298e-04 -0.3654 0.7147872
`trình độ học vấn` -7.5976e-03 4.9063e-04 -15.4852 < 2.2e-16 ***
tin_dung1 2.7949e-02 3.9507e-03 7.0745 1.500e-12 ***
`tỷ lệ phụ thuộc` 6.2562e-04 6.6331e-05 9.4317 < 2.2e-16 ***
`phi NN`1 -2.8961e-02 3.7660e-03 -7.6902 1.469e-14 ***
TD_MNPB1 5.7990e-02 1.4924e-02 3.8858 0.0001020 ***
TN1 2.6390e-02 1.3417e-02 1.9669 0.0491939 *
DNB1 -2.6083e-02 7.2436e-03 -3.6008 0.0003173 ***
DBSCL1 -4.0454e-03 8.9270e-03 -0.4532 0.6504286
BTB_DHMT1 1.9583e-02 1.0777e-02 1.8172 0.0691905 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
dF/dx is for discrete change for the following variables:
[1] "`khu vực`1" "`dân tộc`1" "`giới tính`1" "tin_dung1" "`phi NN`1"
[6] "TD_MNPB1" "TN1" "DNB1" "DBSCL1" "BTB_DHMT1" Xem xét khả năng dự báo của mô hình
library(caret)Loading required package: lattice
Attaching package: 'caret'The following object is masked from 'package:purrr':
liftpredictions <- predict(logit, newdata = DATA, type = "response")
predicted_classes <- ifelse(predictions > 0.5, 1, 0)
conf_matrix <- confusionMatrix(factor(predicted_classes), factor(DATA$`tình trạng nghèo`))
print(conf_matrix)Confusion Matrix and Statistics
Reference
Prediction 0 1
0 10058 897
1 84 136
Accuracy : 0.9122
95% CI : (0.9068, 0.9174)
No Information Rate : 0.9076
P-Value [Acc > NIR] : 0.04548
Kappa : 0.1908
Mcnemar's Test P-Value : < 2e-16
Sensitivity : 0.9917
Specificity : 0.1317
Pos Pred Value : 0.9181
Neg Pred Value : 0.6182
Prevalence : 0.9076
Detection Rate : 0.9000
Detection Prevalence : 0.9803
Balanced Accuracy : 0.5617
'Positive' Class : 0