Trong phần này, chúng ta sẽ bắt đầu bằng việc nạp bộ dữ liệu “Adult Income” và thực hiện các bước chuẩn bị ban đầu. Mục tiêu là đảm bảo dữ liệu sẵn sàng cho việc phân tích, bao gồm việc xử lý các giá trị bị thiếu (nếu có) và hiểu rõ cấu trúc của từng biến. Việc nắm vững các đặc điểm cơ bản của bộ dữ liệu là nền tảng để chúng ta có thể tiến hành các phân tích sâu hơn một cách chính xác và hiệu quả
Bộ dữ liệu “Adult Income” được thu thập từ cuộc điều tra dân số Hoa Kỳ, với mục tiêu là dự đoán xem một cá nhân có thu nhập trên 50.000 USD mỗi năm hay không dựa trên các đặc điểm cá nhân và nghề nghiệp.
Vì vậy trong phạm vi bộ dữ liệu, tôi sẽ phân tích theo 3 nội dung lớn:
Khám phá mối liên hệ giữa các biến định tính và mức thu nhập
Ước lượng khoản tin cậy cho một số biến định lượng
Thực hiện kiểm định giả thuyết thống kê.
# Đọc dữ liệu
dam <- read.csv("C:/Users/damde/Downloads/adult_clean.csv", header = T)
# Thay thế giá trị "?" bằng "Unknown"
dam <- dam %>%
mutate(across(where(is.character), ~ifelse(. == "?", "Unknown", .)))
dim(dam)
## [1] 32561 15
str(dam)
## 'data.frame': 32561 obs. of 15 variables:
## $ age : int 39 50 38 53 28 37 49 52 31 42 ...
## $ workclass : chr "State-gov" "Self-emp-not-inc" "Private" "Private" ...
## $ fnlwgt : int 77516 83311 215646 234721 338409 284582 160187 209642 45781 159449 ...
## $ education : chr "Bachelors" "Bachelors" "HS-grad" "11th" ...
## $ education.num : int 13 13 9 7 13 14 5 9 14 13 ...
## $ marital.status: chr "Never-married" "Married-civ-spouse" "Divorced" "Married-civ-spouse" ...
## $ occupation : chr "Adm-clerical" "Exec-managerial" "Handlers-cleaners" "Handlers-cleaners" ...
## $ relationship : chr "Not-in-family" "Husband" "Not-in-family" "Husband" ...
## $ race : chr "White" "White" "White" "Black" ...
## $ sex : chr "Male" "Male" "Male" "Male" ...
## $ capital.gain : int 2174 0 0 0 0 0 0 0 14084 5178 ...
## $ capital.loss : int 0 0 0 0 0 0 0 0 0 0 ...
## $ hours.per.week: int 40 13 40 40 40 40 16 45 50 40 ...
## $ native.country: chr "United-States" "United-States" "United-States" "United-States" ...
## $ income : chr "<=50K" "<=50K" "<=50K" "<=50K" ...
any(is.na(dam))
## [1] FALSE
Để phục vụ cho phân tích, chúng ta chuyển các biến ký tự thành factor và các biến số cần thiết thành numeric
# Chuyển các cột ký tự thành factor
dam <- data.frame(lapply(dam, function(x) if(is.character(x)) as.factor(x) else x))
# Chuyển các cột số về đúng định dạng numeric
dam$age <- as.numeric(dam$age)
dam$hours.per.week <- as.numeric(dam$hours.per.week)
# Kiểm tra lại cấu trúc
sapply(dam, class)
## age workclass fnlwgt education education.num
## "numeric" "factor" "integer" "factor" "integer"
## marital.status occupation relationship race sex
## "factor" "factor" "factor" "factor" "factor"
## capital.gain capital.loss hours.per.week native.country income
## "integer" "integer" "numeric" "factor" "factor"
Chúng ta sẽ xem xét sự phân bố của các đặc điểm như giới tính, trình độ học vấn, loại công việc,…
giới tính
library(dplyr)
dam %>%
count(education) %>% # Đếm số lượng cho mỗi nhóm education
ggplot(aes(x = n, y = reorder(education, n))) + # Vẽ số lượng (n) trên trục x và sắp xếp education theo n
geom_col(fill = "skyblue") +
labs(title = "Phân bố theo Trình độ học vấn", x = "Số lượng", y = "Học vấn") +
theme_minimal()
Dữ liệu có sự chênh lệch lớn về giới tính, với nam giới chiếm 66.9% và nữ giới chiếm 33.1%.
Trình độ học vấn
library(dplyr)
library(ggplot2)
dam %>%
count(education) %>%
ggplot(aes(x = n, y = reorder(education, n))) +
geom_col(fill = "skyblue") +
labs(title = "Phân bố theo Trình độ học vấn",
x = "Số lượng",
y = "Học vấn") +
theme_minimal()
Nhận xét: Trình độ ‘HS-grad’ (tốt nghiệp THPT), ‘Some-college’ (đang học CĐ/ĐH) và ‘Bachelors’ (Cử nhân) là các nhóm đông nhất.
Thu nhập
library(dplyr)
library(ggplot2)
dam %>%
count(income) %>%
ggplot(aes(x = income, y = n, fill = income)) +
geom_col() +
geom_text(aes(label = n), vjust = -0.5) + # Hiển thị số lượng trên mỗi cột
labs(title = "Phân bố thu nhập",
x = "Mức thu nhập",
y = "Số lượng") +
scale_fill_manual(values = c("<=50K" = "orange", ">50K" = "darkgreen")) +
theme_minimal()
Dữ liệu bị mất cân bằng rõ rệt, với 75.9% dân số có thu nhập <=50K.
Tuổi
summary(dam$age)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 17.00 28.00 37.00 38.58 48.00 90.00
hist(dam$age, col = "lightgreen", main = "Phân bố tuổi", xlab = "Tuổi")
Tuổi của những người tham gia khảo sát trải dài từ 17 đến 90, với độ tuổi trung bình là 38.6. Phân bố tuổi có xu hướng lệch phải.
Số giờ làm việc mỗi tuần
summary(dam$hours.per.week)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 1.00 40.00 40.00 40.44 45.00 99.00
hist(dam$hours.per.week, col = "lightblue", main = "Phân bố số giờ làm việc/tuần", xlab = "Số giờ làm việc")
Hầu hết mọi người làm việc khoảng 40 giờ/tuần.
Mục đích là ước tính độ tuổi trung bình của toàn bộ dân số dựa trên mẫu hiện có
age_conf <- t.test(dam$age, conf.level = 0.95)
age_conf
##
## One Sample t-test
##
## data: dam$age
## t = 510.39, df = 32560, p-value < 2.2e-16
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## 38.43348 38.72981
## sample estimates:
## mean of x
## 38.58165
Chúng ta có thể tin cậy 95% rằng tuổi trung bình thực sự của tổng thể nằm trong khoảng từ 38.43 đến 38.72 tuổi.
Mục tiêu là ước tính tỷ lệ người có thu nhập cao trong toàn bộ dân số.
n_high_income <- sum(dam$income == ">50K")
n_total <- nrow(dam)
prop.test(n_high_income, n_total, conf.level = 0.95)
##
## 1-sample proportions test with continuity correction
##
## data: n_high_income out of n_total, null probability 0.5
## X-squared = 8748.7, df = 1, p-value < 2.2e-16
## alternative hypothesis: true p is not equal to 0.5
## 95 percent confidence interval:
## 0.2361808 0.2454996
## sample estimates:
## p
## 0.2408096
Tỷ lệ người có thu nhập >50K trong mẫu là 24.1%. Chúng ta có thể tin cậy 95% rằng tỷ lệ thực sự trong tổng thể nằm trong khoảng từ 23.6% đến 24.6%.
Giả thuyết (H₀): Giới tính và Thu nhập là hai biến độc lập.
dam_sex_income <- table(dam$sex, dam$income)
addmargins(dam_sex_income)
##
## <=50K >50K Sum
## Female 9592 1179 10771
## Male 15128 6662 21790
## Sum 24720 7841 32561
ggplot(dam, aes(x = sex, fill = income)) +
geom_bar(position = "fill") +
labs(y = "Tỷ lệ", title = "Tỷ lệ thu nhập theo giới tính")
chisq.test(dam_sex_income)
##
## Pearson's Chi-squared test with Yates' continuity correction
##
## data: dam_sex_income
## X-squared = 1517.8, df = 1, p-value < 2.2e-16
p-value < 2.2e-16, rất nhỏ so với 0.05. Chúng ta bác bỏ giả thuyết H₀. Có mối liên hệ rất mạnh mẽ giữa giới tính và thu nhập.
Kiểm định giả thuyết H₁: Tỷ lệ nữ có thu nhập cao < Tỷ lệ nam có thu nhập cao
counts_high_income <- c(dam_sex_income["Female", ">50K"], dam_sex_income["Male", ">50K"])
totals_by_sex <- c(sum(dam_sex_income["Female", ]), sum(dam_sex_income["Male", ]))
prop.test(counts_high_income, totals_by_sex, alternative = "less", correct = FALSE)
##
## 2-sample test for equality of proportions without continuity correction
##
## data: counts_high_income out of totals_by_sex
## X-squared = 1518.9, df = 1, p-value < 2.2e-16
## alternative hypothesis: less
## 95 percent confidence interval:
## -1.0000000 -0.1891457
## sample estimates:
## prop 1 prop 2
## 0.1094606 0.3057366
Tỷ lệ nữ có thu nhập cao (11%) thấp hơn một cách có ý nghĩa thống kê so với nam (30.6%).
Giả thuyết (H₀): Quốc gia gốc (Mỹ/Nước ngoài) và Thu nhập là độc lập.
# (Đây là đoạn code bạn đã viết ở phần "Tính RR và OR")
dam$origin <- ifelse(dam$native.country == "United-States", "US", "Non-US")
tbl_origin_income <- table(dam$origin, dam$income)
# Tính RR
rr_origin <- riskratio(tbl_origin_income)
print("Tỷ số nguy cơ (Risk Ratio):")
## [1] "Tỷ số nguy cơ (Risk Ratio):"
rr_origin
## $data
##
## <=50K >50K Total
## Non-US 2721 670 3391
## US 21999 7171 29170
## Total 24720 7841 32561
##
## $measure
## risk ratio with 95% C.I.
## estimate lower upper
## Non-US 1.000000 NA NA
## US 1.244217 1.159238 1.335426
##
## $p.value
## two-sided
## midp.exact fisher.exact chi.square
## Non-US NA NA NA
## US 2.30395e-10 2.470941e-10 4.969598e-10
##
## $correction
## [1] FALSE
##
## attr(,"method")
## [1] "Unconditional MLE & normal approximation (Wald) CI"
# Tính OR
or_origin <- oddsratio(tbl_origin_income)
print("Tỷ số chênh (Odds Ratio):")
## [1] "Tỷ số chênh (Odds Ratio):"
or_origin
## $data
##
## <=50K >50K Total
## Non-US 2721 670 3391
## US 21999 7171 29170
## Total 24720 7841 32561
##
## $measure
## odds ratio with 95% C.I.
## estimate lower upper
## Non-US 1.00000 NA NA
## US 1.32359 1.212098 1.447236
##
## $p.value
## two-sided
## midp.exact fisher.exact chi.square
## Non-US NA NA NA
## US 2.30395e-10 2.470941e-10 4.969598e-10
##
## $correction
## [1] FALSE
##
## attr(,"method")
## [1] "median-unbiased estimate & mid-p exact CI"
RR = 1.24: Người ở Mỹ có nguy cơ thu nhập >50K cao hơn 1.24 lần so với người nước ngoài.
OR = 1.32: Cơ hội (odds) có thu nhập >50K của người ở Mỹ cao hơn 1.32 lần so với người nước ngoài.
Vì khoảng tin cậy của cả RR và OR đều không chứa 1, mối quan hệ này có ý nghĩa thống kê
Giả thuyết (H₀): Số giờ làm việc trung bình giữa nam và nữ là như nhau.
a) Biểu độ Boxplot so sánh:
ggplot(dam, aes(x = sex, y = hours.per.week, fill = sex)) +
geom_boxplot() +
labs(title = "So sánh số giờ làm việc/tuần giữa nam và nữ")
Phân bố chung: Hộp (box) của nhóm Nam được đặt ở vị trí cao hơn so với nhóm Nữ, cho thấy rằng 50% số nam giới ở khoảng giữa làm việc nhiều giờ hơn so với 50% số nữ giới tương ứng.
Giá trị trung vị (Median): Đường kẻ đậm bên trong hộp của nhóm Nam (khoảng 40 giờ) cao hơn một chút so với của nhóm Nữ.
Ngoại lệ (Outliers): Cả hai nhóm đều có rất nhiều điểm ngoại lệ, cho thấy sự biến động lớn trong số giờ làm việc. Nhiều người làm việc ít hơn hoặc nhiều hơn đáng kể so với mức trung bình của nhóm mình.
b) Kiểm định T-test hai mẫu độc lập
Trong kiểm định này, chúng ta muốn so sánh giá trị trung bình về số giờ làm việc của hai nhóm (Nam và Nữ).
H₀ (Giả thuyết không): Không có sự khác biệt về số giờ làm việc trung bình mỗi tuần giữa nam và nữ.
H₁ (Giả thuyết thay thế): Có sự khác biệt về số giờ làm việc trung bình mỗi tuần giữa nam và nữ.
hours_test <- t.test(hours.per.week ~ sex, data = dam)
hours_test
##
## Welch Two Sample t-test
##
## data: hours.per.week by sex
## t = -42.882, df = 21958, p-value < 2.2e-16
## alternative hypothesis: true difference in means between group Female and group Male is not equal to 0
## 95 percent confidence interval:
## -6.292787 -5.742664
## sample estimates:
## mean in group Female mean in group Male
## 36.41036 42.42809
p-value < 2.2e-16, rất nhỏ. Chúng ta bác bỏ giả thuyết H₀. Có sự khác biệt có ý nghĩa thống kê về số giờ làm việc trung bình giữa nam và nữ. Cụ thể, nam giới làm việc trung bình (42.4 giờ) nhiều hơn nữ giới (36.4 giờ).
Phần này là bước nâng cao, giống như phần cuối của bài mẫu. Mục tiêu là xây dựng mô hình dự đoán khả năng một người có thu nhập >50K hay không dựa trên các đặc điểm khác.
Chúng ta sẽ xây dựng mô hình dự đoán income dựa trên age, education.num, và hours.per.week.
# Tạo biến mục tiêu nhị phân (0/1)
dam$income_bin <- ifelse(dam$income == ">50K", 1, 0)
# Chạy mô hình hồi quy logistic
model_logit <- glm(income_bin ~ age + education.num + hours.per.week,
data = dam,
family = binomial(link = "logit"))
# Xem kết quả
summary(model_logit)
##
## Call:
## glm(formula = income_bin ~ age + education.num + hours.per.week,
## family = binomial(link = "logit"), data = dam)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -8.523161 0.110169 -77.36 <2e-16 ***
## age 0.046911 0.001160 40.45 <2e-16 ***
## education.num 0.345334 0.006480 53.29 <2e-16 ***
## hours.per.week 0.042834 0.001268 33.78 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 35948 on 32560 degrees of freedom
## Residual deviance: 28991 on 32557 degrees of freedom
## AIC: 28999
##
## Number of Fisher Scoring iterations: 5
Tiếp tục xây dựng mô hình Probit
# Mô hình hồi quy Probit
model_probit <- glm(income_bin ~ age + education.num + hours.per.week,
data = dam,
family = binomial(link = "probit"))
# Xem kết quả
summary(model_probit)
##
## Call:
## glm(formula = income_bin ~ age + education.num + hours.per.week,
## family = binomial(link = "probit"), data = dam)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -4.9488518 0.0595976 -83.04 <2e-16 ***
## age 0.0276700 0.0006545 42.27 <2e-16 ***
## education.num 0.1974943 0.0036587 53.98 <2e-16 ***
## hours.per.week 0.0248457 0.0007268 34.19 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 35948 on 32560 degrees of freedom
## Residual deviance: 28930 on 32557 degrees of freedom
## AIC: 28938
##
## Number of Fisher Scoring iterations: 5
Xây dựng mô hình Cloglog
# Xây dựng mô hình CLogLog
model_cloglog <- glm(income_bin ~ age + education.num + hours.per.week,
data = dam,
family = binomial(link = "cloglog"))
# Xem kết quả
summary(model_cloglog)
##
## Call:
## glm(formula = income_bin ~ age + education.num + hours.per.week,
## family = binomial(link = "cloglog"), data = dam)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -6.6503394 0.0800237 -83.11 <2e-16 ***
## age 0.0325705 0.0008993 36.22 <2e-16 ***
## education.num 0.2690124 0.0049211 54.66 <2e-16 ***
## hours.per.week 0.0271375 0.0009283 29.23 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 35948 on 32560 degrees of freedom
## Residual deviance: 29396 on 32557 degrees of freedom
## AIC: 29404
##
## Number of Fisher Scoring iterations: 10
Ý nghĩa thống kê: Trong cả ba mô hình, p-value (Pr(>|z|)) của các biến age, education.num, và hours.per.week đều cực kỳ nhỏ (< 2e-16). Điều này khẳng định cả ba biến đều là những yếu tố dự báo có ý nghĩa thống kê rất cao đến thu nhập.
Chiều hướng tác động: Dấu của các hệ số (Estimate) trong cả hai mô hình đều là dương (+). Điều này có nghĩa là khi tuổi, trình độ học vấn, hoặc số giờ làm việc tăng lên thì khả năng có thu nhập cao (>50K) cũng tăng lên.
Giá trị hệ số (Estimate) của ba mô hình là khác nhau (ví dụ, hệ số của age trong Logit là 0.0469, còn trong Probit là 0.0276). Điều này là bình thường vì chúng được tính trên hai thang đo toán học khác nhau (Logit dùng log-odds, Probit dùng z-score)
Mức độ phù hợp AIC. AIC của Logit là 28999, còn của Probit là 28938, Cloglog là 29404. Có thể thấy chỉ số AIC của Probit thấp hơn một chút vậy nên Probit phù hợp hơn với dữ liệu hơn 1 chút.
Cả hai mô hình đều là những mô hình tốt để dự báo thu nhập. Chúng cùng chỉ ra rằng học vấn (education.num) có tác động mạnh nhất, theo sau là tuổi tác (age) và số giờ làm việc (hours.per.week).
Mặc dù mô hình Probit có chỉ số AIC tốt hơn một chút, mô hình Logit thường được ưa chuộng hơn trong thực tế vì hệ số của nó có thể chuyển đổi thành Tỷ số chênh (Odds Ratio), giúp việc diễn giải trở nên trực quan hơn.