Phần 1: Tìm hiểu và Chuẩn bị dữ liệu

1.1. Giới thiệu và mục tiêu

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

1.2. Đọc và làm sạch dữ liệu

# Đọ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", .)))

1.3. Cấu trúc tổng quan của dữ liệu

a) Số quan sát và số biến

dim(dam)
## [1] 32561    15

b) Danh sách và kiểu dữ liệu của các biến

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

c) Kiểm tra dữ liệu thiếu (NA):

any(is.na(dam))
## [1] FALSE

1.4. Chuyển đổi kiểu dữ liệu

Để 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"

Phần 2: Phân tích mô tả

2.1. Phân tích biến định tính

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.

2.2. Phân tích định lượng

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.

Phần 3: Ưóc lượng và kiểm định cho một biến

3.1. Ưóc lượng khoảng tin cậy 95% cho tuổi trung bình

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.

3.2. Ước lượng khoảng tin cậy 95% cho tỷ lệ người có thu nhập > 50K.

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

Phần 4: Phân tích mối quan hệ giữa hai biến

4.1. Mối quan hệ giữa hai biến định tính

Giới tính (sex) và Thu nhập (income)

Giả thuyết (H₀): Giới tính và Thu nhập là hai biến độc lập.

a) Bảng chéo và Biểu đồ:

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")

b) Kiểm định chi bình phương:

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.

c) So sánh tỷ lệ (Hiệu tỷ lệ):

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%).

Quốc gia gốc (native.country) và Thu nhập (income)

Giả thuyết (H₀): Quốc gia gốc (Mỹ/Nước ngoài) và Thu nhập là độc lập.

Bảng chéo, Tỷ số nguy cơ (RR) và Tỷ số chênh (OR):

# (Đâ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ê

4.2. Mối quan hệ giữa biến định tính và định lượng

Giới tính (sex) và Số giờ làm việc (hours.per.week)

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 5: Hồi quy Logistic/Probit (Dự đoán thu nhập)

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.

5.1. Xây dựng mô hình

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

Diễn giải kết quả Probit và Logistic, Cloglog

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

Phân biệt điểm khác biệt

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.

Kết luận chung

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.