Phần 1: Giới thiệu

Thu nhập và nghề nghiệp luôn là hai chủ đề trung tâm trong các nghiên cứu về thị trường lao động, phản ánh trực tiếp mức sống, khả năng phát triển cá nhân và sự đóng góp vào tăng trưởng kinh tế - xã hội. Trong bối cảnh nền kinh tế hiện đại ngày càng cạnh tranh và biến động, việc nhận diện và phân tích các yếu tố tác động đến thu nhập và nghề nghiệp không chỉ có ý nghĩa về mặt học thuật mà còn mang tính thực tiễn cao trong việc hoạch định chính sách phát triển nguồn nhân lực, nâng cao chất lượng lao động và đảm bảo công bằng trong tiếp cận cơ hội nghề nghiệp.

Để thực hiện nghiên cứu này, em sử dụng Bộ dữ liệu Adult Census Income, một tập dữ liệu kinh điển, được trích xuất từ cuộc điều tra dân số năm 1994 do Cục Điều tra Dân số Hoa Kỳ thực hiện. Bộ dữ liệu bao gồm 30.162 quan sát với 15 biến đặc trưng, phản ánh đa dạng các đặc điểm về nhân khẩu học, kinh tế và xã hội của các cá nhân. Cụ thể như sau:

Age: Độ tuổi của cá nhân tham gia khảo sát.

Workclass: loại hình đơn vị làm việc của cá nhân.

Final weight: Trọng số mẫu

Education: Trình độ học vấn cao nhất mà cá nhân đạt được.

Education num: Phản ánh thứ bậc tương ứng với các trình độ học vấn.

Marital status: Tình trạng hôn nhân của cá nhân.

Occupation: Nghề nghiệp hiện tại của cá nhân.

Relationship: Mô tả vai trò của cá nhân trong mối quan hệ với các thành viên trong hộ gia đình, chẳng hạn như vợ/chồng, con cái, hay chủ hộ.

Race: Biểu thị chủng tộc của cá nhân.

Sex: Biến định tính phản ánh giới tính của cá nhân.

Review Rating: Đánh giá từ khách hàng.

Capital gain: Thể hiện tổng thu nhập từ lợi nhuận vốn (capital gain) mà cá nhân thu được trong năm. Đây là khoản thu nhập không thường xuyên nhưng có thể tạo sự khác biệt lớn trong tổng thu nhập cá nhân.

Capital loss: Ngược lại với capital gain, biến này phản ánh tổng mức thua lỗ từ các khoản đầu tư tài chính trong năm.

Hours per week: Cho biết số giờ làm việc trung bình mỗi tuần.

Native country: Quốc gia nơi cá nhân sinh ra.

Income: thể hiện mức thu nhập hàng năm của cá nhân.

#Đọc dữ liệu
a <- read.csv("C:/Users/Admin/Desktop/adultcensusincome.csv", 
                 header = TRUE, sep = ",")
#cấu trúc của dữ liệu 
str(a)
## 'data.frame':    30162 obs. of  15 variables:
##  $ age           : int  82 54 41 34 38 74 68 45 38 52 ...
##  $ workclass     : chr  "Private" "Private" "Private" "Private" ...
##  $ final.weight  : int  132870 140359 264663 216864 150601 88638 422013 172274 164526 129177 ...
##  $ education     : chr  "HS-grad" "7th-8th" "Some-college" "HS-grad" ...
##  $ education.num : int  9 4 10 9 6 16 9 16 15 13 ...
##  $ marital.status: chr  "Widowed" "Divorced" "Separated" "Divorced" ...
##  $ occupation    : chr  "Exec-managerial" "Machine-op-inspct" "Prof-specialty" "Other-service" ...
##  $ relationship  : chr  "Not-in-family" "Unmarried" "Own-child" "Unmarried" ...
##  $ race          : chr  "White" "White" "White" "White" ...
##  $ sex           : chr  "Female" "Female" "Female" "Female" ...
##  $ capital.gain  : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ capital.loss  : int  4356 3900 3900 3770 3770 3683 3683 3004 2824 2824 ...
##  $ hours.per.week: int  18 40 40 45 40 20 40 35 45 20 ...
##  $ native.country: chr  "United-States" "United-States" "United-States" "United-States" ...
##  $ income        : chr  "<=50K" "<=50K" "<=50K" "<=50K" ...
names(a)
##  [1] "age"            "workclass"      "final.weight"   "education"     
##  [5] "education.num"  "marital.status" "occupation"     "relationship"  
##  [9] "race"           "sex"            "capital.gain"   "capital.loss"  
## [13] "hours.per.week" "native.country" "income"

Trong phạm vi bài làm này, chỉ thực hiện thao tác với các biến định tính bao gồm: education(trình độ giáo dục), income(thu nhập), occupation(nghề nghiệp), sex(giới tính) .

# Chọn các biến định tính
a1 <- c("education","occupation", "sex","income")
print(a1)
## [1] "education"  "occupation" "sex"        "income"
# Tạo bộ dữ liệu mới chỉ chứa các biến định tính
b <- a[, a1]
str(b)
## 'data.frame':    30162 obs. of  4 variables:
##  $ education : chr  "HS-grad" "7th-8th" "Some-college" "HS-grad" ...
##  $ occupation: chr  "Exec-managerial" "Machine-op-inspct" "Prof-specialty" "Other-service" ...
##  $ sex       : chr  "Female" "Female" "Female" "Female" ...
##  $ income    : chr  "<=50K" "<=50K" "<=50K" "<=50K" ...

Tiếp theo, em tiến hành Kiểm tra xem có giá trị thiếu (NA) trong các cột định tính hay không.

#Kiểm tra xem có giá trị thiếu (NA) trong các cột định tính
colSums(is.na(b))
##  education occupation        sex     income 
##          0          0          0          0

Tất cả các cột định tính trong dataframe b (education, occupation, sex,income) đều có giá trị 0 cho số lượng giá trị thiếu (NA). Điều này có nghĩa là không có giá trị NA trong các cột định tính của bộ dữ liệu b.

Trước khi tiến hành các phân tích thống kê, việc chuyển đổi các biến định tính từ kiểu ký tự (character) sang kiểu phân loại (factor) là cần thiết nhằm đảm bảo rằng các hàm thống kê và mô hình phân tích sẽ xử lý chúng đúng cách. Trong R, các biến factor không chỉ giúp biểu diễn dữ liệu phân loại một cách chính xác hơn mà còn tối ưu hóa hiệu suất tính toán và cho phép áp dụng các phương pháp phân tích định tính như kiểm định Chi-bình phương hay mô hình hồi quy phân loại. Việc kiểm tra và xử lý kiểu dữ liệu phù hợp là bước tiền xử lý quan trọng, góp phần đảm bảo tính chính xác và hiệu quả của toàn bộ quy trình phân tích.

# Kiểm tra từng cột trong dataframe b
sapply(b, is.factor)
##  education occupation        sex     income 
##      FALSE      FALSE      FALSE      FALSE
# Chuyển đổi các cột character sang factor
b$education <- as.factor(b$education)
b$occupation <- as.factor(b$occupation)
b$sex <- as.factor(b$sex)
b$income <- as.factor(b$income)
#kiểm tra lại
sapply(b, is.factor)
##  education occupation        sex     income 
##       TRUE       TRUE       TRUE       TRUE

Phần 2: Các yếu tố tác động đến thu nhập

Thu nhập cá nhân không chỉ phản ánh kết quả của quá trình lao động mà còn chịu sự chi phối của nhiều yếu tố khác nhau, trong đó trình độ học vấn và nghề nghiệp là hai thành tố then chốt. Phân tích mối quan hệ giữa các yếu tố này với thu nhập giúp hiểu rõ hơn về tính phân tầng thu nhập trong xã hội và cơ hội kinh tế của các nhóm dân cư.

2.1 Ảnh hưởng của trình độ học vấn đến khả năng đạt được mức thu nhập cao

Trong nền kinh tế hiện đại, học vấn được xem là chìa khóa mở ra cơ hội việc làm tốt hơn và mức thu nhập cao hơn. Người có trình độ học vấn cao thường có khả năng đảm nhận các vị trí đòi hỏi kiến thức chuyên môn sâu, kỹ năng giải quyết vấn đề, tư duy phản biện và khả năng sáng tạo. Ngược lại, những người có trình độ học vấn thấp thường bị giới hạn trong việc tiếp cận các vị trí có mức lương cao hoặc nghề nghiệp ổn định. Do đó, việc phân tích mối liên hệ giữa trình độ học vấn và khả năng đạt được thu nhập cao là bước đầu tiên để kiểm chứng luận điểm này.

# Tạo bảng tần suất chéo giữa education và income
table1 <- table(b$education, b$income)

# Tính tỷ lệ phần trăm theo hàng (row percentages)
prop_table <- prop.table(table1, margin = 1) * 100
percent_50K <- round(prop_table[, ">50K"], 2)  # Lấy cột >50K, làm tròn 2 chữ số

# Tạo bảng kết hợp: tần suất, tổng, và tỷ lệ % >50K
table_combined <- cbind(table1, Tổng = rowSums(table1), `% >50K` = percent_50K)
# Sử dụng kableExtra để định dạng bảng HTML
library(kableExtra)

kable(table_combined, format = "html", 
      caption = "Bảng 2.1: Tần số và Tỷ lệ phần trăm theo trình độ học vấn và mức thu nhập") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"), 
                full_width = FALSE) %>%
  add_footnote("Ghi chú: Tỷ lệ % >50K được tính bằng số cá nhân có thu nhập >50K chia cho tổng số cá nhân trong mỗi nhóm trình độ học vấn.", 
               notation = "none") %>%
  column_spec(1, width = "15em") %>%  # Điều chỉnh độ rộng cột Education
  column_spec(2:5, width = "8em", extra_css = "text-align: center;") 
Bảng 2.1: Tần số và Tỷ lệ phần trăm theo trình độ học vấn và mức thu nhập
<=50K >50K Tổng % >50K
10th 761 59 820 7.20
11th 989 59 1048 5.63
12th 348 29 377 7.69
1st-4th 145 6 151 3.97
5th-6th 276 12 288 4.17
7th-8th 522 35 557 6.28
9th 430 25 455 5.49
Assoc-acdm 752 256 1008 25.40
Assoc-voc 963 344 1307 26.32
Bachelors 2918 2126 5044 42.15
Doctorate 95 280 375 74.67
HS-grad 8223 1617 9840 16.43
Masters 709 918 1627 56.42
Preschool 45 0 45 0.00
Prof-school 136 406 542 74.91
Some-college 5342 1336 6678 20.01
Ghi chú: Tỷ lệ % >50K được tính bằng số cá nhân có thu nhập >50K chia cho tổng số cá nhân trong mỗi nhóm trình độ học vấn.

Giải thích các phân loại trong biến Education:

  • Preschool: Mẫu giáo.

  • 1st-4th: Lớp 1-4.

  • 5th-6th: Lớp 5-6.

  • 7th-8th: Lớp 7-8.

  • 9th: Lớp 9.

  • 10th: Lớp 10.

  • 11th: Lớp 11.

  • 12th: Lớp 12.

  • HS-grad: Tốt nghiệp THPT.

  • Some-college: Học một số tín chỉ đại học.

  • Assoc-voc: Bằng cao đẳng nghề.

  • Assoc-acdm: Cao đẳng học thuật.

  • Bachelors: Bằng cử nhân.

  • Masters: Bằng thạc sĩ.

  • Prof-school: Bằng chuyên nghiệp (cao hơn thạc sĩ).

  • Doctorate: Bằng tiến sĩ.

Bảng tần suất chéo giữa trình độ học vấn và mức thu nhập cho thấy mối quan hệ rõ rệt giữa giáo dục và khả năng đạt thu nhập cao (>50K USD/năm). Dữ liệu chỉ ra rằng trình độ học vấn càng cao, tỷ lệ đạt thu nhập >50K càng lớn. Những người chỉ học đến mẫu giáo (Preschool) hoặc tiểu học (1st-4th, 5th-6th) có tỷ lệ >50K rất thấp, từ 0% đến 4.17%, phản ánh hạn chế trong cơ hội việc làm lương cao. Ngược lại, các nhóm có trình độ đại học trở lên, như Bachelors (42.15%), Masters (56.42%), Prof-school (74.91%), và Doctorate (74.67%), có tỷ lệ thu nhập cao vượt trội, gấp nhiều lần so với trung bình toàn mẫu (36.15%). Các bằng liên kết (Assoc-voc: 26.32%, Assoc-acdm: 25.40%) và tốt nghiệp trung học (HS-grad: 16.43%) cũng cải thiện thu nhập, nhưng không bằng trình độ cao hơn. Dữ liệu nhấn mạnh tầm quan trọng của giáo dục trong việc nâng cao thu nhập, đặc biệt là giáo dục đại học và chuyên môn.

library(ggplot2)
library(dplyr)
library(tidyr)

# Tạo bảng tần suất chéo giữa education và income
table1 <- table(b$education, b$income)

# Chuyển bảng tần suất chéo thành dataframe với các cột Less_Equal_50K, Greater_50K, và Total
data_df <- as.data.frame(table1) %>%
  dplyr::rename(education = Var1, income = Var2) %>%  # Đổi tên cột để rõ ràng
  tidyr::pivot_wider(names_from = income, values_from = Freq, values_fill = 0) %>%  # Chuyển sang dạng rộng
  dplyr::rename(Less_Equal_50K = `<=50K`, Greater_50K = `>50K`) %>%  # Đổi tên cột
  mutate(Total = Less_Equal_50K + Greater_50K)

# Tính % >50K tự động
library(dplyr)
data_df <- data_df %>%
  mutate(Percent_50K = (Greater_50K / Total) * 100)

# Sắp xếp dữ liệu theo Percent_50K tăng dần
data_df$education <- factor(data_df$education, 
                           levels = data_df$education[order(data_df$Percent_50K)])

# Tạo biểu đồ lollipop
library(ggplot2)

ggplot(data_df, aes(x = Percent_50K, y = education)) +
  geom_segment(aes(x = 0, xend = Percent_50K, y = education, yend = education), 
               color = "steelblue", size = 1.5) +  # Đường que kẹo
  geom_point(color = "steelblue", size = 4, shape = 21, fill = "white") +  # Điểm kẹo
  geom_text(aes(label = sprintf("%.1f%%", Percent_50K)), 
            hjust = -0.2, size = 3.5, color = "black") +  # Nhãn dữ liệu
  labs(title = "Tỷ lệ Thu nhập >50K theo trình độ Học vấn",
       x = "Tỷ lệ % >50K",
       y = "Trình độ Học vấn") +
  theme_minimal() +
  theme(axis.text.x = element_text(size = 10),
        axis.text.y = element_text(size = 10),
        plot.title = element_text(hjust = 0.5, size = 14, face = "bold")) +
  scale_x_continuous(limits = c(0, 85), breaks = seq(0, 85, by = 10)) +
  coord_cartesian(clip = "off")  # Đảm bảo nhãn dữ liệu không bị cắt

Kiểm định thống kê

Để kiểm tra xem liệu trình độ học vấn của cá nhân có mối quan hệ thống kê với khả năng đạt được mức thu nhập cao hay không, nghiên cứu tiến hành áp dụng kiểm định Chi-bình phương về tính độc lập giữa hai biến phân loại: trình độ học vấn và thu nhập. Kiểm định này cho phép đánh giá giả thuyết liệu sự khác biệt về thu nhập giữa các nhóm trình độ học vấn có xảy ra một cách ngẫu nhiên hay thực sự phản ánh một mối liên hệ có ý nghĩa thống kê trong tổng thể.

Đưa ra giả thuyết:

\[ \left\{ \begin{array}{ll} H_0: \text{Trình độ học vấn của cá nhân không có mối quan hệ với mức thu nhập.} \\ H_1: \text{Trình độ học vấn của cá nhân có mối quan hệ với mức thu nhập.} \end{array} \right. \]

# Thực hiện kiểm định Chi bình phương
chi_test <- chisq.test(table1)

# Hiển thị kết quả kiểm định
print(chi_test)
## 
##  Pearson's Chi-squared test
## 
## data:  table1
## X-squared = 4070.4, df = 15, p-value < 2.2e-16

Vì p-value < 2.2e-16 <0.05, ta bác bỏ giả thuyết H₀. Điều này có nghĩa là trình độ học vấn của cá nhân có mối quan hệ với mức thu nhập.

So sánh tỷ lệ thu nhập >50K giữa 2 nhóm trình độ học vấn

Sau khi đã xác định mối liên hệ tổng quát giữa trình độ học vấn và khả năng đạt mức thu nhập cao bằng kiểm định Chi-bình phương, bước tiếp theo cần thiết là tiến hành so sánh cụ thể tỷ lệ thu nhập >50K giữa hai nhóm trình độ học vấn tiêu biểu. Việc lựa chọn hai nhóm cụ thể để so sánh giúp phân tích sâu hơn về mức độ chênh lệch giữa các trình độ học vấn và tác động của chúng đến cơ hội thu nhập cao.

Để đơn giản hóa phân tích và phù hợp với mục tiêu kiểm định tỷ lệ, trình độ học vấn được phân nhóm thành hai loại, thông qua phép kiểm định hiệu hai tỷ lệ, em sẽ kiểm tra xem sự khác biệt về tỷ lệ đạt mức thu nhập >50K giữa hai nhóm này có mang ý nghĩa thống kê hay không.

  • Nhóm học vấn cao (Higher Education): Bao gồm các cá nhân có trình độ từ Cử nhân trở lên (Bachelors, Masters, Prof-school, Doctorate).

  • Nhóm học vấn thấp (Lower Education): Bao gồm các cá nhân còn lại.

Giả thuyết kiểm định được đặt ra như sau:

  • Giả thuyết rỗng (H₀): Tỷ lệ người có thu nhập cao giữa hai nhóm học vấn là bằng nhau.

\[ H_0: p_{Lower} = p_{Higher} \]

  • Giả thuyết đối (H₁): Tỷ lệ người có thu nhập cao giữa hai nhóm học vấn là khác nhau.

\[ H_1: p_{Lower} \neq p_{Higher} \]

Trong đó:
\(p_{Lower}\) là tỷ lệ người có thu nhập cao trong nhóm có trình độ học vấn thấp.
\(p_{Higher}\) là tỷ lệ người có thu nhập cao trong nhóm có trình độ học vấn cao.

# Nhóm trình độ học vấn thành 2 nhóm
b$edu_group <- ifelse(b$education %in% c("Bachelors", "Masters", "Doctorate", "Prof-school"), 
                      "Higher", "Lower")

# Tạo bảng tần suất chéo giữa nhóm học vấn và thu nhập
table_edu_income <- table(b$edu_group, b$income)
print(table_edu_income)
##         
##          <=50K  >50K
##   Higher  3858  3730
##   Lower  18796  3778
# Thực hiện kiểm định hiệu hai tỷ lệ
prop.test(table_edu_income[, ">50K"], rowSums(table_edu_income))
## 
##  2-sample test for equality of proportions with continuity correction
## 
## data:  table_edu_income[, ">50K"] out of rowSums(table_edu_income)
## X-squared = 3191, df = 1, p-value < 2.2e-16
## alternative hypothesis: two.sided
## 95 percent confidence interval:
##  0.3118596 0.3365503
## sample estimates:
##    prop 1    prop 2 
## 0.4915656 0.1673607

Kết quả kiểm định thu được cho thấy giá trị thống kê Chi bình phương (X-squared) đạt 3191 với một bậc tự do (df = 1). Đặc biệt, giá trị p-value < 2.2e-16, nhỏ hơn rất nhiều so với mức ý nghĩa 0.05, điều này cho phép bác bỏ giả thuyết H₀. Nói cách khác, với mức ý nghĩa 5%, có bằng chứng thống kê mạnh mẽ cho thấy trình độ học vấn có ảnh hưởng đến khả năng đạt mức thu nhập trên 50K USD/năm.

Relative Risk - RR

Bên cạnh kiểm định sự khác biệt về tỷ lệ, việc tính toán rủi ro tương đối (Relative Risk - RR) đóng vai trò quan trọng trong việc lượng hóa mức độ ảnh hưởng của trình độ học vấn đến khả năng đạt được mức thu nhập cao. Relative Risk giúp đo lường xem người có trình độ học vấn thấp có khả năng đạt mức thu nhập cao như thế nào so với người có trình độ học vấn cao hơn

library(epitools)

# Tính RR từ bảng table_edu_income
rr_result <- riskratio(table_edu_income)
print(rr_result)
## $data
##         
##          <=50K >50K Total
##   Higher  3858 3730  7588
##   Lower  18796 3778 22574
##   Total  22654 7508 30162
## 
## $measure
##         risk ratio with 95% C.I.
##           estimate    lower     upper
##   Higher 1.0000000       NA        NA
##   Lower  0.3404646 0.328092 0.3533037
## 
## $p.value
##         two-sided
##          midp.exact fisher.exact chi.square
##   Higher         NA           NA         NA
##   Lower           0            0          0
## 
## $correction
## [1] FALSE
## 
## attr(,"method")
## [1] "Unconditional MLE & normal approximation (Wald) CI"

Kết quả phân tích cho thấy tỷ lệ cá nhân có trình độ học vấn cao (Higher) đạt mức thu nhập trên 50K là 49.14% (3730/7588), trong khi tỷ lệ này ở nhóm có trình độ học vấn thấp (Lower) chỉ là 16.73% (3778/22574).

Giá trị Risk Ratio (RR) = 0.3405 cho thấy xác suất để người thuộc nhóm có trình độ học vấn thấp đạt mức thu nhập cao chỉ bằng khoảng 34.05% so với nhóm có trình độ học vấn cao. Nói cách khác, người có trình độ học vấn thấp có nguy cơ đạt mức thu nhập cao thấp hơn khoảng 66% so với người có trình độ học vấn cao.

Tỷ số Chênh (Odds Ratio - OR):

# Tính odds từng nhóm
odds_higher <- table_edu_income["Higher", ">50K"] / table_edu_income["Higher", "<=50K"]
odds_lower <- table_edu_income["Lower", ">50K"] / table_edu_income["Lower", "<=50K"]

odds_higher
## [1] 0.9668222
odds_lower
## [1] 0.2010002

Odds của nhóm Higher (trình độ học vấn cao) là 0.9668, nghĩa là cứ 1 người có thu nhập >50K thì có khoảng 1.03 người không đạt được mức đó. (1 / 0.9668 ≈ 1.03).

Odds của nhóm Lower (trình độ học vấn thấp) là 0.2010, nghĩa là cứ 1 người có thu nhập >50K thì có khoảng 4.81 người không đạt mức đó. Điều này cho thấy khả năng đạt thu nhập cao ở nhóm học vấn thấp là rất thấp.

# Tính Odds Ratio bằng hàm oddsratio()
library(epitools)
result_or <- oddsratio(table_edu_income, method = "wald", conf.level = 0.95)

# In kết quả Odds Ratio
print(result_or)
## $data
##         
##          <=50K >50K Total
##   Higher  3858 3730  7588
##   Lower  18796 3778 22574
##   Total  22654 7508 30162
## 
## $measure
##         odds ratio with 95% C.I.
##           estimate     lower     upper
##   Higher 1.0000000        NA        NA
##   Lower  0.2078978 0.1963829 0.2200879
## 
## $p.value
##         two-sided
##          midp.exact fisher.exact chi.square
##   Higher         NA           NA         NA
##   Lower           0            0          0
## 
## $correction
## [1] FALSE
## 
## attr(,"method")
## [1] "Unconditional MLE & normal approximation (Wald) CI"

Kết quả phân tích Odds Ratio cho thấy sự khác biệt đáng kể về khả năng đạt mức thu nhập cao giữa hai nhóm trình độ học vấn. Giá trị Odds Ratio = 0.2079 nghĩa là tỷ lệ giữa số người đạt thu nhập cao so với số người không đạt thu nhập cao trong nhóm có trình độ học vấn thấp chỉ bằng khoảng 20.8% so với tỷ lệ đó ở nhóm có trình độ học vấn cao.

2.2 Ảnh hưởng của nghề nghiệp đến khả năng đạt được thu nhập cao

Không chỉ dừng lại ở học vấn, bản thân nghề nghiệp là yếu tố trực tiếp quyết định mức thu nhập của cá nhân. Những người làm việc trong các ngành nghề yêu cầu trình độ chuyên môn cao, vị trí quản lý, hoặc có vai trò quan trọng trong tổ chức thường nhận được mức lương cao hơn đáng kể so với những người làm việc trong các ngành nghề phổ thông, lao động chân tay hoặc dịch vụ đơn giản. Việc kiểm tra mối liên hệ giữa nghề nghiệp và khả năng đạt được mức thu nhập cao sẽ giúp xác định rõ hơn tính phân hóa thu nhập giữa các nhóm nghề nghiệp.

# Đặt mã hóa UTF-8 để hỗ trợ tiếng Việt
Sys.setlocale("LC_ALL", "en_US.UTF-8")
## [1] "LC_COLLATE=en_US.UTF-8;LC_CTYPE=en_US.UTF-8;LC_MONETARY=en_US.UTF-8;LC_NUMERIC=C;LC_TIME=en_US.UTF-8"
# Tạo bảng tần suất chéo giữa education và income
table1 <- table(b$occupation, b$income)

# Tính tỷ lệ phần trăm theo hàng (row percentages)
prop_table <- prop.table(table1, margin = 1) * 100
percent_50K <- round(prop_table[, ">50K"], 2)  # Lấy cột >50K, làm tròn 2 chữ số

# Tạo bảng kết hợp: tần suất, tổng, và tỷ lệ % >50K
table_combined <- cbind(table1, Tổng = rowSums(table1), `% >50K` = percent_50K)
# Sử dụng kableExtra để định dạng bảng HTML
library(kableExtra)

kable(table_combined, format = "html", 
      caption = "Bảng 2.2: Tần số và Tỷ lệ phần trăm theo nghề nghiệp và mức thu nhập") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"), 
                full_width = FALSE) %>%
  add_footnote("Ghi chú: Tỷ lệ % >50K được tính bằng số cá nhân có thu nhập >50K chia cho tổng số cá nhân trong mỗi nhóm nghề nghiệp.", 
               notation = "none") %>%
  column_spec(1, width = "15em") %>%  # Điều chỉnh độ rộng cột Education
  column_spec(2:5, width = "8em", extra_css = "text-align: center;") 
Bảng 2.2: Tần số và Tỷ lệ phần trăm theo nghề nghiệp và mức thu nhập
<=50K >50K Tổng % >50K
Adm-clerical 3223 498 3721 13.38
Armed-Forces 8 1 9 11.11
Craft-repair 3122 908 4030 22.53
Exec-managerial 2055 1937 3992 48.52
Farming-fishing 874 115 989 11.63
Handlers-cleaners 1267 83 1350 6.15
Machine-op-inspct 1721 245 1966 12.46
Other-service 3080 132 3212 4.11
Priv-house-serv 142 1 143 0.70
Prof-specialty 2227 1811 4038 44.85
Protective-serv 434 210 644 32.61
Sales 2614 970 3584 27.06
Tech-support 634 278 912 30.48
Transport-moving 1253 319 1572 20.29
Ghi chú: Tỷ lệ % >50K được tính bằng số cá nhân có thu nhập >50K chia cho tổng số cá nhân trong mỗi nhóm nghề nghiệp.

Giải thích các phân loại trong biến :

  • Adm-clerical: Hành chính - Văn phòng.

  • Armed-Forces: Lực lượng vũ trang.

  • Craft-repair: Thợ thủ công - Sửa chữa.

  • Exec-managerial: Quản lý - Điều hành.

  • Farming-fishing: Nông nghiệp - Ngư nghiệp.

  • Handlers-cleaners: Bốc xếp - Dọn dẹp.

  • Machine-op-inspct: Vận hành máy móc - Kiểm tra.

  • Other-service: Dịch vụ khác.

  • Priv-house-serv: Dịch vụ gia đình tư nhân.

  • Prof-specialty: Chuyên gia chuyên môn.

  • Protective-serv: Dịch vụ bảo vệ.

  • Sales: Bán hàng.

  • Tech-support: Hỗ trợ kỹ thuật.

  • Transport-moving: Vận tải - Di chuyển.

Dựa trên Bảng 2.2, có thể nhận thấy sự khác biệt rõ rệt về mức thu nhập giữa các nhóm nghề nghiệp. Các nhóm nghề như Exec-managerial (Quản lý điều hành) và Prof-specialty (Chuyên môn nghiệp vụ chuyên sâu) có tỷ lệ cá nhân có thu nhập trên 50K cao nổi bật, lần lượt đạt 48.52% và 44.85%, phản ánh đặc thù của các nghề này yêu cầu trình độ chuyên môn cao và thường đi kèm với mức thu nhập tốt. Ngược lại, các nhóm nghề lao động phổ thông như Other-service (Dịch vụ khác), Priv-house-serv (Dịch vụ gia đình tư nhân), hay Handlers-cleaners (Bốc xếp - Dọn dẹp) có tỷ lệ >50K rất thấp, chỉ dao động từ 0.7% đến 6.15%, cho thấy mức thu nhập của nhóm nghề này còn khá hạn chế. Đáng chú ý, nhóm Protective-serv (Dịch vụ bảo vệ) và Tech-support cũng có tỷ lệ trên 50K tương đối cao (32.61% và 30.48%), phản ánh mức đãi ngộ tốt hơn so với các nhóm nghề nghiệp mang tính lao động chân tay. Nhìn chung, bảng số liệu cho thấy thu nhập cao thường gắn liền với những nghề nghiệp đòi hỏi kỹ năng chuyên môn, quản lý hoặc kỹ thuật cao.

Trực quan hoá

library(ggplot2)
library(dplyr)
library(tidyr)

# Tạo bảng tần suất chéo giữa occupation và income
table_occ_income <- table(b$occupation, b$income)

# Chuyển bảng tần suất chéo thành dataframe với các cột Less_Equal_50K, Greater_50K, và Total
data_occ <- as.data.frame(table_occ_income) %>%
  dplyr::rename(occupation = Var1, income = Var2) %>%
  tidyr::pivot_wider(names_from = income, values_from = Freq, values_fill = 0) %>%
  dplyr::rename(Less_Equal_50K = `<=50K`, Greater_50K = `>50K`) %>%
  mutate(Total = Less_Equal_50K + Greater_50K)

# Tính % >50K
data_occ <- data_occ %>%
  mutate(Percent_50K = (Greater_50K / Total) * 100)

# Sắp xếp occupation theo thứ tự tăng dần của tỷ lệ >50K
data_occ$occupation <- factor(data_occ$occupation, 
                              levels = data_occ$occupation[order(data_occ$Percent_50K)])

# Trực quan hóa bằng biểu đồ lollipop
ggplot(data_occ, aes(x = Percent_50K, y = occupation)) +
  geom_segment(aes(x = 0, xend = Percent_50K, y = occupation, yend = occupation), 
               color = "steelblue", size = 1.5) +
  geom_point(color = "steelblue", size = 4, shape = 21, fill = "white") +
  geom_text(aes(label = sprintf("%.1f%%", Percent_50K)), 
            hjust = -0.2, size = 3.5, color = "black") +
  labs(title = "Tỷ lệ Thu nhập >50K theo Nghề nghiệp",
       x = "Tỷ lệ % >50K",
       y = "Nghề nghiệp") +
  theme_minimal() +
  theme(axis.text.x = element_text(size = 10),
        axis.text.y = element_text(size = 10),
        plot.title = element_text(hjust = 0.5, size = 14, face = "bold")) +
  scale_x_continuous(limits = c(0, max(data_occ$Percent_50K) + 5), breaks = seq(0, 80, by = 10)) +
  coord_cartesian(clip = "off")

Kiểm định thống kê

Trong nền kinh tế hiện đại, nghề nghiệp là một trong những yếu tố then chốt quyết định đến mức thu nhập của mỗi cá nhân. Tuy nhiên, để xác định xem sự khác biệt về thu nhập giữa các nhóm nghề nghiệp chỉ là do ngẫu nhiên hay thực sự phản ánh một mối liên hệ có ý nghĩa thống kê, cần phải tiến hành kiểm định thống kê một cách khách quan.

Để làm rõ vấn đề này, nghiên cứu sử dụng kiểm định Chi-bình phương độc lập (Chi-squared Test of Independence) nhằm kiểm tra giả thuyết:

\[ \left\{ \begin{array}{ll} H_0: \text{Nghề nghiệp của cá nhân không có mối quan hệ với mức thu nhập.} \\ H_1: \text{Nghề nghiệp của cá nhân có mối quan hệ với khả năng đạt thu nhập cao hay thấp.} \end{array} \right. \]

# Thực hiện kiểm định Chi bình phương
chi_test <- chisq.test(table1)

# Hiển thị kết quả kiểm định
print(chi_test)
## 
##  Pearson's Chi-squared test
## 
## data:  table1
## X-squared = 3687.6, df = 13, p-value < 2.2e-16

Vì p-value < 2.2e-16 <0.05, ta bác bỏ giả thuyết H₀. Điều này có nghĩa là nghề nghiệp của cá nhân có mối quan hệ với khả năng đạt thu nhập cao hay thấp.

so sánh tỷ lệ thu nhập cao (>50K) giữa 2 nhóm nghề nghiệp “Exec-managerial”(Quản lý và điều hành) và “Craft-repair”(Thủ công và sửa chữa)

Trong cấu trúc nghề nghiệp, sự phân hóa về thu nhập thường phản ánh sự khác biệt về yêu cầu kỹ năng, trình độ chuyên môn và trách nhiệm công việc. Nhằm làm rõ mức độ chênh lệch thu nhập giữa các nhóm nghề nghiệp, nghiên cứu tiến hành so sánh tỷ lệ thu nhập cao (>50K) giữa hai nhóm nghề tiêu biểu:

  • Exec-managerial (Quản lý và điều hành): Đại diện cho nhóm nghề nghiệp có tính chất quản lý, đòi hỏi kỹ năng điều hành và chịu trách nhiệm cao.

  • Craft-repair (Thủ công và sửa chữa): Đại diện cho nhóm nghề có tính chất lao động kỹ thuật, thường yêu cầu tay nghề nhưng mức độ quản lý thấp hơn.

Đặt giả thuyết: \[ \left\{ \begin{array}{ll} H_0: & p_1 = p_2 \quad (\text{Tỷ lệ người có thu nhập cao giữa hai nhóm nghề là bằng nhau}) \\ H_1: & p_1 \neq p_2 \quad (\text{Tỷ lệ người có thu nhập cao giữa hai nhóm nghề là khác nhau}) \end{array} \right. \]

table_occ_income_2groups <- table(b$occupation, b$income)
table(b$occupation)
## 
##      Adm-clerical      Armed-Forces      Craft-repair   Exec-managerial 
##              3721                 9              4030              3992 
##   Farming-fishing Handlers-cleaners Machine-op-inspct     Other-service 
##               989              1350              1966              3212 
##   Priv-house-serv    Prof-specialty   Protective-serv             Sales 
##               143              4038               644              3584 
##      Tech-support  Transport-moving 
##               912              1572
table(b$occupation, b$income)
##                    
##                     <=50K >50K
##   Adm-clerical       3223  498
##   Armed-Forces          8    1
##   Craft-repair       3122  908
##   Exec-managerial    2055 1937
##   Farming-fishing     874  115
##   Handlers-cleaners  1267   83
##   Machine-op-inspct  1721  245
##   Other-service      3080  132
##   Priv-house-serv     142    1
##   Prof-specialty     2227 1811
##   Protective-serv     434  210
##   Sales              2614  970
##   Tech-support        634  278
##   Transport-moving   1253  319
library(dplyr)

# Bước 1: Lọc 2 nhóm nghề
b_subset <- b %>%
  filter(occupation %in% c("Exec-managerial", "Craft-repair")) %>%
  droplevels()  # Loại bỏ các levels cũ không cần thiết

# Bước 2: Tạo bảng tần suất đúng, chỉ có 2 nghề
table_occ_income <- table(b_subset$occupation, b_subset$income)

# Kiểm tra lại bảng đã đúng chưa
print(table_occ_income)
##                  
##                   <=50K >50K
##   Craft-repair     3122  908
##   Exec-managerial  2055 1937
# Bước 3: Thực hiện kiểm định hiệu hai tỷ lệ
prop.test(x = table_occ_income[, ">50K"],
          n = rowSums(table_occ_income),
          alternative = "two.sided",
          correct = TRUE)
## 
##  2-sample test for equality of proportions with continuity correction
## 
## data:  table_occ_income[, ">50K"] out of rowSums(table_occ_income)
## X-squared = 590.79, df = 1, p-value < 2.2e-16
## alternative hypothesis: two.sided
## 95 percent confidence interval:
##  -0.2803274 -0.2394931
## sample estimates:
##    prop 1    prop 2 
## 0.2253102 0.4852204

Vì < 2.2e-16 0.05, ta bác bỏ H₀. Điều này có nghĩa là có bằng chứng thống kê khẳng định rằng tỷ lệ người có thu nhập cao giữa hai nhóm nghề này là khác nhau.

Relative Risk - RR

Để đo lường mức độ khác biệt về khả năng đạt được mức thu nhập cao giữa hai nhóm nghề nghiệp “Quản lý và điều hành” (Exec-managerial) và “Thủ công và sửa chữa” (Craft-repair), nghiên cứu sử dụng chỉ số (Relative Risk - RR).

# Đảm bảo đã load thư viện cần thiết
library(dplyr)
library(epitools)

# Giả sử bạn đã đọc dữ liệu vào biến b
# b <- read.csv("adultcensusincome.csv")  # nếu chưa load lại thì bỏ comment

# Tạo bảng tần suất giữa occupation và income cho 2 nghề cần so sánh
table_occ_income <- table(b$occupation, b$income)

# Lọc 2 nghề: "Craft-repair" và "Exec-managerial"
table_2occ <- table_occ_income[c("Craft-repair", "Exec-managerial"), c(">50K", "<=50K")]

# Kiểm tra bảng 2x2 đã đúng chưa
print(table_2occ)
##                  
##                   >50K <=50K
##   Craft-repair     908  3122
##   Exec-managerial 1937  2055
# Tính Relative Risk (RR) dùng hàm riskratio
result_rr <- riskratio(table_2occ, rev = "columns")

# Xem kết quả
print(result_rr)
## $data
##                  
##                   <=50K >50K Total
##   Craft-repair     3122  908  4030
##   Exec-managerial  2055 1937  3992
##   Total            5177 2845  8022
## 
## $measure
##                  risk ratio with 95% C.I.
##                   estimate    lower   upper
##   Craft-repair    1.000000       NA      NA
##   Exec-managerial 2.153566 2.016903 2.29949
## 
## $p.value
##                  two-sided
##                   midp.exact  fisher.exact    chi.square
##   Craft-repair            NA            NA            NA
##   Exec-managerial          0 7.477099e-133 9.567349e-131
## 
## $correction
## [1] FALSE
## 
## attr(,"method")
## [1] "Unconditional MLE & normal approximation (Wald) CI"

Kết quả phân tích cho thấy sự khác biệt đáng kể về khả năng đạt được mức thu nhập cao giữa hai nhóm nghề nghiệp. Tỷ lệ cá nhân có thu nhập >50K trong nhóm “Quản lý và điều hành” (Exec-managerial) là 48.52% (1937/3992), trong khi tỷ lệ này ở nhóm “Thủ công và sửa chữa” (Craft-repair) chỉ đạt 22.53% (908/4030).

Giá trị Relative Risk (RR) = 2.154 cho biết xác suất đạt được mức thu nhập cao của những người làm nghề “Quản lý và điều hành” cao hơn khoảng 2.15 lần so với những người làm nghề “Thủ công và sửa chữa”.

Odds và tỷ số Chênh (Odds Ratio - OR)

# Tính odds cho từng nhóm
odds_craft <- table_2occ["Craft-repair", ">50K"] / table_2occ["Craft-repair", "<=50K"]
odds_exec <- table_2occ["Exec-managerial", ">50K"] / table_2occ["Exec-managerial", "<=50K"]

# In kết quả
cat("Odds - Craft-repair: ", odds_craft, "\n")
## Odds - Craft-repair:  0.2908392
cat("Odds - Exec-managerial: ", odds_exec, "\n")
## Odds - Exec-managerial:  0.9425791

Odds được hiểu là tỷ lệ giữa số người xảy ra biến cố (thu nhập >50K) và số người không xảy ra biến cố (thu nhập ≤50K). Từ kết quả trên:

Odds của nhóm “Thủ công và sửa chữa” (Craft-repair) là 0.29, nghĩa là cứ 1 người làm nghề “Thủ công và sửa chữa” có thu nhập >50K thì có khoảng 3.44 người không đạt mức thu nhập này (1 / 0.29 ≈ 3.44). Điều đó cho thấy khả năng có thu nhập cao trong nhóm này thấp hơn rõ rệt so với khả năng không có thu nhập cao.

Odds của nhóm “Quản lý và điều hành” (Exec-managerial) là 0.94, nghĩa là cứ 1 người làm nghề “Quản lý và điều hành” có thu nhập >50K thì có khoảng 1.06 người không đạt mức thu nhập này (1 / 0.94 ≈ 1.06). Nói cách khác, khả năng có thu nhập cao gần như ngang bằng với khả năng không có thu nhập cao trong nhóm nghề nghiệp này.Tóm lại: Nghề “Quản lý và điều hành” có khả năng đạt được mức thu nhập cao cao gấp nhiều lần so với nghề “Thủ công và sửa chữa”.

library(epitools)

# Tạo lại bảng 2x2 với Craft-repair nằm trên
table_2x2 <- matrix(c(
  table_2occ["Craft-repair", "<=50K"], table_2occ["Craft-repair", ">50K"],
  table_2occ["Exec-managerial", "<=50K"], table_2occ["Exec-managerial", ">50K"]
), nrow = 2, byrow = TRUE)

# Đặt tên cho dễ đọc
dimnames(table_2x2) <- list(
  Occupation = c("Craft-repair", "Exec-managerial"),
  Income = c("<=50K", ">50K")
)

# Tính Odds Ratio (method = "wald") với khoảng tin cậy 95%
result_or <- oddsratio(table_2x2, method = "wald", conf.level = 0.95)

# In kết quả
print(result_or)
## $data
##                  Income
## Occupation        <=50K >50K Total
##   Craft-repair     3122  908  4030
##   Exec-managerial  2055 1937  3992
##   Total            5177 2845  8022
## 
## $measure
##                  odds ratio with 95% C.I.
## Occupation        estimate    lower    upper
##   Craft-repair    1.000000       NA       NA
##   Exec-managerial 3.240894 2.942743 3.569254
## 
## $p.value
##                  two-sided
## Occupation        midp.exact  fisher.exact    chi.square
##   Craft-repair            NA            NA            NA
##   Exec-managerial          0 7.477099e-133 9.567349e-131
## 
## $correction
## [1] FALSE
## 
## attr(,"method")
## [1] "Unconditional MLE & normal approximation (Wald) CI"

Tỷ số odds (Odds Ratio) giữa hai nhóm nghề nghiệp là 3.24, nghĩa là tỷ lệ giữa số người đạt thu nhập cao so với số người không đạt thu nhập cao trong nhóm “Quản lý và điều hành” cao hơn khoảng 3.24 lần so với tỷ lệ đó trong nhóm “Thủ công và sửa chữa”.

Điều này phản ánh rằng làm việc trong nhóm nghề nghiệp “Quản lý và điều hành” có lợi thế rõ rệt hơn về khả năng đạt mức thu nhập cao so với nhóm “Thủ công và sửa chữa”.

Phần 3: Các yếu tố tác động đến nghề nghiệp

Bên cạnh vai trò quyết định mức thu nhập, nghề nghiệp của mỗi cá nhân còn chịu sự chi phối bởi nhiều yếu tố khác, trong đó đặc biệt phải kể đến trình độ học vấn và giới tính. Việc nhận diện và đo lường mức độ tác động của các yếu tố này đến cơ hội lựa chọn nghề nghiệp không chỉ góp phần lý giải cơ chế phân tầng nghề nghiệp trong xã hội, mà còn cung cấp căn cứ khoa học để phân tích bất bình đẳng giới và bất bình đẳng cơ hội trên thị trường lao động.

3.1 Ảnh hưởng của trình độ Học vấn đến nghề nghiệp

Học vấn không chỉ quyết định mức thu nhập mà còn ảnh hưởng trực tiếp đến khả năng tiếp cận các ngành nghề khác nhau. Người có trình độ học vấn cao thường có xu hướng lựa chọn các ngành nghề chuyên môn cao, kỹ thuật hoặc quản lý; trong khi những người có học vấn thấp có xu hướng tham gia vào các ngành nghề phổ thông hoặc dịch vụ đơn giản. Việc phân tích mối liên hệ giữa học vấn và nghề nghiệp giúp nhận diện rõ vai trò của học vấn trong sự phân hóa nghề nghiệp trên thị trường lao động.

library(dplyr)
library(knitr)
# Đặt mã hóa UTF-8 để hỗ trợ tiếng Việt
Sys.setlocale("LC_ALL", "en_US.UTF-8")
## [1] "LC_COLLATE=en_US.UTF-8;LC_CTYPE=en_US.UTF-8;LC_MONETARY=en_US.UTF-8;LC_NUMERIC=C;LC_TIME=en_US.UTF-8"
# Bảng tần suất chéo (tính theo hàng): mỗi hàng là 100%
table_edu_occ <- table(b$education, b$occupation)
freq_edu_occ <- prop.table(table_edu_occ, margin = 1) * 100  # Tỷ lệ phần trăm theo từng hàng

# Đưa về dạng dataframe để kable đẹp hơn
df_freq_edu_occ <- as.data.frame.matrix(round(freq_edu_occ, 2))  # Làm tròn 2 chữ số thập phân

# Tạo bảng đẹp với kable
df_freq_edu_occ %>%
  kable(format = "html", caption = "Bảng 3.1: Bảng Tần Suất Chéo (%): Trình độ học vấn và nghề nghiệp") %>%
  kable_styling(full_width = TRUE, bootstrap_options = c("striped", "hover", "condensed", "responsive"), font_size = 12)
Bảng 3.1: Bảng Tần Suất Chéo (%): Trình độ học vấn và nghề nghiệp
Adm-clerical Armed-Forces Craft-repair Exec-managerial Farming-fishing Handlers-cleaners Machine-op-inspct Other-service Priv-house-serv Prof-specialty Protective-serv Sales Tech-support Transport-moving
10th 4.63 0.00 19.88 2.93 5.24 8.66 12.32 23.54 0.73 1.10 0.73 9.76 0.37 10.12
11th 6.39 0.00 16.51 3.15 3.53 11.64 9.45 22.71 1.34 1.91 0.67 13.36 0.57 8.78
12th 9.55 0.27 14.59 3.45 4.24 9.81 9.02 22.02 1.06 2.39 1.33 11.67 0.80 9.81
1st-4th 0.00 0.00 15.23 1.99 11.92 10.60 15.23 25.17 7.28 2.65 0.66 4.64 0.00 4.64
5th-6th 2.08 0.00 14.24 0.35 12.50 13.54 17.36 20.83 4.86 0.35 0.35 3.82 0.35 9.38
7th-8th 1.97 0.00 19.93 3.23 12.39 8.26 16.52 16.34 1.44 1.44 1.62 5.21 0.90 10.77
9th 3.08 0.00 20.44 2.64 6.15 10.77 16.48 21.54 2.20 0.66 0.88 7.03 0.44 7.69
Assoc-acdm 18.95 0.00 11.31 14.19 1.39 2.18 3.27 7.74 0.20 13.49 3.37 14.09 7.14 2.68
Assoc-voc 12.62 0.00 18.97 11.48 3.98 2.14 4.74 8.42 0.31 13.01 3.67 7.96 9.64 3.06
Bachelors 9.87 0.02 4.36 26.43 1.53 0.95 1.17 3.37 0.12 28.97 1.98 15.72 4.42 1.09
Doctorate 1.07 0.00 0.53 14.40 0.27 0.00 0.27 0.00 0.00 80.53 0.00 2.13 0.53 0.27
HS-grad 13.72 0.04 19.34 8.12 4.08 6.13 10.27 12.73 0.48 2.33 2.18 10.69 1.60 8.29
Masters 4.12 0.06 1.29 30.18 0.61 0.31 0.49 1.11 0.00 50.34 0.92 7.93 2.09 0.55
Preschool 4.44 0.00 8.89 0.00 20.00 4.44 24.44 31.11 4.44 2.22 0.00 0.00 0.00 0.00
Prof-school 1.48 0.00 1.29 8.86 0.74 0.00 0.00 0.74 0.00 81.73 0.18 3.14 1.29 0.55
Some-college 18.93 0.03 12.76 13.03 2.61 3.92 4.60 11.44 0.22 6.33 2.96 14.91 4.06 4.19

Dựa trên Bảng 3.1, có thể nhận thấy mối quan hệ khá rõ ràng giữa trình độ học vấn và nghề nghiệp. Các nghề yêu cầu trình độ chuyên môn cao như Exec-managerial và Prof-specialty có tỷ lệ cao những người sở hữu bằng Bachelors, Masters, Doctorate, hoặc Prof-school. Cụ thể, nghề Prof-specialty có đến 80.53% người có trình độ Doctorate và 81.73% có trình độ Prof-school, trong khi Exec-managerial có tỷ lệ cao nhất ở nhóm Masters (30.18%) và Bachelors (26.43%). Ngược lại, các nghề lao động phổ thông hoặc dịch vụ như Handlers-cleaners, Other-service, hay Priv-house-serv lại chủ yếu tập trung ở nhóm có trình độ HS-grad, 12th, hoặc thậm chí thấp hơn, như 7th-8th và 9th. Đặc biệt, các nhóm nghề như Preschool hay 1st-4th cho thấy tỷ lệ cao ở các nghề lao động giản đơn hoặc dịch vụ cá nhân. Bảng số liệu cho thấy rõ ràng rằng trình độ học vấn càng cao thường gắn liền với khả năng tham gia vào các ngành nghề có tính chuyên môn hoặc quản lý cao, từ đó tạo điều kiện cho mức thu nhập cao hơn như đã thể hiện ở Bảng 2.2. Điều này nhấn mạnh vai trò quan trọng của giáo dục đối với cơ hội nghề nghiệp và thu nhập trong xã hội.

Trực quan hóa

# Chuyển thành dataframe để ggplot sử dụng
df_heatmap <- as.data.frame(as.table(freq_edu_occ))
colnames(df_heatmap) <- c("Education", "Occupation", "Percent")

# Vẽ heatmap
ggplot(df_heatmap, aes(x = Occupation, y = Education, fill = Percent)) +
  geom_tile(color = "white") +
  scale_fill_gradient(low = "white", high = "steelblue", name = "%") +
  geom_text(aes(label = sprintf("%.1f", Percent)), size = 3) +
  labs(title = "Tỷ lệ (%) Nghề nghiệp theo Trình độ Học vấn",
       x = "Nghề nghiệp",
       y = "Trình độ học vấn") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, size = 9),
        axis.text.y = element_text(size = 9),
        plot.title = element_text(size = 14, face = "bold", hjust = 0.5))

library(ggplot2)
library(dplyr)

# Tạo bảng tần suất
df_bar <- as.data.frame(table(b$education, b$occupation))
colnames(df_bar) <- c("Education", "Occupation", "Count")

# Vẽ biểu đồ cột nhóm
ggplot(df_bar, aes(x = Education, y = Count, fill = Occupation)) +
  geom_bar(stat = "identity", position = "fill") +  # position = "fill" để thành phần trăm
  scale_y_continuous(labels = scales::percent) +
  labs(title = "Phân bổ Nghề nghiệp theo Trình độ học vấn (%)",
       x = "Trình độ học vấn",
       y = "Tỷ lệ (%)",
       fill = "Nghề nghiệp") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

Kiểm định thống kê

Trong nghiên cứu về thị trường lao động, trình độ học vấn thường được xem là yếu tố nền tảng quyết định đến khả năng lựa chọn và tiếp cận nghề nghiệp. Tuy nhiên, câu hỏi đặt ra là liệu sự khác biệt về nghề nghiệp giữa các nhóm trình độ học vấn có xảy ra một cách ngẫu nhiên hay thực sự phản ánh một mối liên hệ có ý nghĩa thống kê trong tổng thể?

Để kiểm tra vấn đề này, nghiên cứu tiến hành kiểm định Chi-bình phương độc lập (Chi-squared Test of Independence) với giả thuyết như sau:

\[ \left\{ \begin{array}{ll} H_0: & \text{Trình độ học vấn không có mối quan hệ với nghề nghiệp.} \\ H_1: & \text{Trình độ học vấn có mối quan hệ với nghề nghiệp.} \end{array} \right. \]

# Tạo bảng tần suất chéo giữa education và occupation
table_edu_occ <- table(b$education, b$occupation)

# Kiểm định Chi-bình phương
chisq_test <- chisq.test(table_edu_occ)

# Xem kết quả
chisq_test
## 
##  Pearson's Chi-squared test
## 
## data:  table_edu_occ
## X-squared = 15359, df = 195, p-value < 2.2e-16

Với p-value rất nhỏ (< 0.05), chúng ta bác bỏ giả thuyết H₀. Điều này cho thấy: Có mối quan hệ có ý nghĩa thống kê giữa trình độ học vấn và nghề nghiệp của các cá nhân. Nói cách khác, trình độ học vấn có ảnh hưởng rõ rệt đến nghề nghiệp mà một cá nhân đảm nhiệm.

So sánh tỷ lệ làm nghề “Exec-managerial” giữa 2 nhóm trình độ học vấn

Trong thị trường lao động, những vị trí thuộc nhóm “Exec-managerial” (Quản lý và điều hành) thường gắn liền với yêu cầu về trình độ học vấn cao do tính chất công việc đòi hỏi khả năng lãnh đạo, phân tích và ra quyết định. Để kiểm tra xem trình độ học vấn có làm gia tăng cơ hội đảm nhiệm các vị trí quản lý hay không, nghiên cứu tiến hành so sánh tỷ lệ làm nghề “Exec-managerial” giữa hai nhóm trình độ học vấn chính:

  • Nhóm học vấn cao (Higher Education): Bao gồm các cá nhân có trình độ từ Cử nhân trở lên (Bachelors, Masters, Prof-school, Doctorate).

  • Nhóm học vấn thấp (Lower Education): Bao gồm các cá nhân còn lại.

Phép kiểm định thống kê tỷ lệ hai mẫu độc lập được sử dụng nhằm kiểm tra sự khác biệt về tỷ lệ giữa hai nhóm, với giả thuyết:

\[ \left\{ \begin{array}{ll} H_0: & p_{Higher} = p_{Lower} \ (\text{Tỷ lệ làm nghề Exec-managerial giữa 2 trình độ là bằng nhau}) \\ H_1: & p_{Higher} \ne p_{Lower} \ (\text{Tỷ lệ làm nghề Exec-managerial giữa 2 trình độ là khác nhau}) \end{array} \right. \]

b$edu_group <- ifelse(b$education %in% c("Bachelors", "Masters", "Doctorate", "Prof-school"), 
                      "Higher", "Lower")

# Nhóm nghề thành: "Exec-managerial" vs. Others
b$occ_group <- ifelse(b$occupation == "Exec-managerial", "Exec-managerial", "Others")

# Bảng tần suất 2x2
table_edu_exec <- table(b$edu_group, b$occ_group)
print(table_edu_exec)
##         
##          Exec-managerial Others
##   Higher            1926   5662
##   Lower             2066  20508
prop.test(table_edu_exec[, "Exec-managerial"], rowSums(table_edu_exec))
## 
##  2-sample test for equality of proportions with continuity correction
## 
## data:  table_edu_exec[, "Exec-managerial"] out of rowSums(table_edu_exec)
## X-squared = 1301.3, df = 1, p-value < 2.2e-16
## alternative hypothesis: two.sided
## 95 percent confidence interval:
##  0.1517230 0.1728782
## sample estimates:
##     prop 1     prop 2 
## 0.25382182 0.09152122

Kết quả: Vì p-value rất nhỏ < 0.05, → Bác bỏ H0, nghĩa là trình độ học vấn có ảnh hưởng rõ rệt đến khả năng làm nghề Exec-managerial.

Tỷ lệ nghề Exec-managerial ở nhóm “Bachelors” cao hơn nhóm “HS-grad” ≈ 18%.

Relative Risk - RR

Sau khi xác định có sự khác biệt về tỷ lệ đảm nhiệm nghề “Exec-managerial” giữa hai nhóm trình độ học vấn, việc tính toán rủi ro tương đối (Relative Risk - RR) sẽ giúp lượng hóa mức độ khác biệt về khả năng làm nghề quản lý giữa nhóm có trình độ cao và nhóm có trình độ thấp hơn.

Rủi ro tương đối (RR) được hiểu là tỷ lệ xác suất làm nghề quản lý ở nhóm có trình độ ” Lower” so với nhóm “Higher”

library(epitools)

# Giả sử bạn đã có table_edu_exec như sau:
table_edu_exec <- table(b$edu_group, b$occ_group)

# Tạo lại bảng 2x2 với đúng thứ tự mong muốn:
table_2x2 <- matrix(c(
  table_edu_exec["Higher", "Others"], table_edu_exec["Higher", "Exec-managerial"],
  table_edu_exec["Lower", "Others"], table_edu_exec["Lower", "Exec-managerial"]
), nrow = 2, byrow = TRUE)

# Đặt tên hàng, cột cho dễ hiểu:
dimnames(table_2x2) <- list(
  Education = c("Higher", "Lower"),
  Occupation = c("Others", "Exec-managerial")
)

# Tính Relative Risk (RR) với khoảng tin cậy 95%
result_rr <- riskratio(table_2x2, method = "wald", conf.level = 0.95)

# In kết quả
print(result_rr)
## $data
##          Occupation
## Education Others Exec-managerial Total
##    Higher   5662            1926  7588
##    Lower   20508            2066 22574
##    Total   26170            3992 30162
## 
## $measure
##          risk ratio with 95% C.I.
## Education  estimate     lower     upper
##    Higher 1.0000000        NA        NA
##    Lower  0.3605727 0.3408098 0.3814816
## 
## $p.value
##          two-sided
## Education midp.exact  fisher.exact    chi.square
##    Higher         NA            NA            NA
##    Lower           0 9.637841e-255 2.937971e-285
## 
## $correction
## [1] FALSE
## 
## attr(,"method")
## [1] "Unconditional MLE & normal approximation (Wald) CI"

Tỷ số Risk Ratio = 0.8007 cho thấy xác suất để người thuộc nhóm có trình độ học vấn thấp làm nghề “Quản lý và điều hành” chỉ bằng khoảng 36.06% so với nhóm có trình độ học vấn cao.Điều này có nghĩa rằng, người có trình độ học vấn thấp có khả năng tham gia vào nghề “Quản lý và điều hành” thấp hơn khoảng 64% so với người có trình độ học vấn cao.

Nhìn chung, kết quả này cho thấy trình độ học vấn có ảnh hưởng rõ rệt đến khả năng đảm nhận các vị trí nghề nghiệp “Quản lý và điều hành”. Những người có trình độ học vấn cao có lợi thế lớn hơn trong việc tiếp cận các vị trí nghề nghiệp có tính chất quản lý, lãnh đạo, so với những người có trình độ học vấn thấp.

Odds và Oddratio

# Bước 1: Tạo bảng 2x2 với cột "Exec-managerial" nằm ở sau
table_2x2 <- matrix(c(
  table(b$edu_group, b$occ_group)["Higher", "Others"], table(b$edu_group, b$occ_group)["Higher", "Exec-managerial"],
  table(b$edu_group, b$occ_group)["Lower", "Others"], table(b$edu_group, b$occ_group)["Lower", "Exec-managerial"]
), nrow = 2, byrow = TRUE)

# Đặt tên cho hàng và cột
dimnames(table_2x2) <- list(
  Education = c("Higher", "Lower"),
  Occupation = c("Others", "Exec-managerial")
)

# In ra bảng để kiểm tra
print(table_2x2)
##          Occupation
## Education Others Exec-managerial
##    Higher   5662            1926
##    Lower   20508            2066
# Bước 2: Tính Odds từng nhóm
odds_higher <- table_2x2["Higher", "Exec-managerial"] / table_2x2["Higher", "Others"]
odds_lower <- table_2x2["Lower", "Exec-managerial"] / table_2x2["Lower", "Others"]

# In kết quả Odds
cat("Odds nhóm Higher:", round(odds_higher, 4), "\n")
## Odds nhóm Higher: 0.3402
cat("Odds nhóm Lower:", round(odds_lower, 4), "\n")
## Odds nhóm Lower: 0.1007

Odds được hiểu là tỷ lệ giữa số người xảy ra biến cố “làm nghề Quản lý và điều hành” và số người không xảy ra biến cố đó.

Odds của nhóm học vấn cao (Higher): 0.3402.Có nghĩa là trong nhóm người có trình độ học vấn cao, cứ 1 người làm nghề khác thì có khoảng 0.34 người làm nghề “Quản lý và điều hành”. Nói cách khác, trong nhóm người có trình độ học vấn cao, cứ 1 người làm nghề “Quản lý và điều hành” thì có khoảng 2.94(1/0.342) người làm các nghề khác.

Odds của nhóm học vấn thấp (Lower): 0.1007.Tức là trong nhóm người có trình độ học vấn thấp, cứ 1 người làm nghề khác thì chỉ có khoảng 0.10 người làm nghề “Quản lý và điều hành”. Nói cách khác, trong nhóm người có trình độ học vấn thấp, cứ 1 người làm nghề “Quản lý và điều hành” thì có khoảng 9.93 (1/0.1007) người làm các nghề khác.

library(epitools)

# Tạo lại bảng 2x2 với cột "Exec-managerial" nằm sau
table_2x2 <- matrix(c(
  table(b$edu_group, b$occ_group)["Higher", "Others"], table(b$edu_group, b$occ_group)["Higher", "Exec-managerial"],
  table(b$edu_group, b$occ_group)["Lower", "Others"], table(b$edu_group, b$occ_group)["Lower", "Exec-managerial"]
), nrow = 2, byrow = TRUE)

dimnames(table_2x2) <- list(
  Education = c("Higher", "Lower"),
  Occupation = c("Others", "Exec-managerial")
)

# Tính Odds Ratio với khoảng tin cậy 95%
result_or <- oddsratio(table_2x2, method = "wald", conf.level = 0.95)

# In kết quả
print(result_or)
## $data
##          Occupation
## Education Others Exec-managerial Total
##    Higher   5662            1926  7588
##    Lower   20508            2066 22574
##    Total   26170            3992 30162
## 
## $measure
##          odds ratio with 95% C.I.
## Education estimate     lower     upper
##    Higher 1.000000        NA        NA
##    Lower  0.296156 0.2764933 0.3172171
## 
## $p.value
##          two-sided
## Education midp.exact  fisher.exact    chi.square
##    Higher         NA            NA            NA
##    Lower           0 9.637841e-255 2.937971e-285
## 
## $correction
## [1] FALSE
## 
## attr(,"method")
## [1] "Unconditional MLE & normal approximation (Wald) CI"

Tỷ số odds (Odds Ratio) giữa hai nhóm học vấn là 0.2962, nghĩa là tỷ lệ giữa số người làm nghề “Quản lý và điều hành” so với số người làm nghề khác trong nhóm có trình độ học vấn thấp chỉ bằng khoảng 29.62% so với tỷ lệ đó trong nhóm có trình độ học vấn cao.

Nói cách khác, khả năng làm nghề “Quản lý và điều hành” của nhóm có trình độ học vấn cao cao hơn khoảng 3.38 lần so với nhóm có trình độ học vấn thấp.

3.2 Ảnh hưởng của giới tính tới cơ hội nghề nghiệp

Giới tính từ lâu đã ảnh hưởng đáng kể đến sự lựa chọn nghề nghiệp, đặc biệt trong các nền văn hóa có sự phân hóa vai trò giới truyền thống. Một số ngành nghề thường mặc định gắn với giới tính nam (ví dụ: kỹ thuật, cơ khí), trong khi một số nghề khác thường thấy nữ giới chiếm ưu thế (ví dụ: giáo dục, chăm sóc, dịch vụ). Việc kiểm tra mối liên hệ giữa giới tính và cơ hội nghề nghiệp giúp làm rõ vấn đề phân hóa giới trong lựa chọn nghề nghiệp và cơ hội thăng tiến.

library(dplyr)
library(knitr)
library(kableExtra)

# Bảng tần số chéo
freq_table <- table(b$sex, b$occupation)

# Tính bảng tần suất theo hàng (giới tính)
prop_table <- prop.table(freq_table, margin = 1) * 100  # Tính theo hàng (%)

# Chuyển thành dataframe để hiển thị với kable
prop_df <- as.data.frame.matrix(round(prop_table, 2))  # Làm tròn 2 chữ số

# In bảng tần suất chéo đẹp (không có cột Total)
kable(prop_df, caption = "Bảng 3.2: Bảng Tần suất chéo (%): Giới tính và Nghề nghiệp") %>%
  kable_styling(full_width = FALSE, position = "center", bootstrap_options = c("striped", "hover", "condensed"))
Bảng 3.2: Bảng Tần suất chéo (%): Giới tính và Nghề nghiệp
Adm-clerical Armed-Forces Craft-repair Exec-managerial Farming-fishing Handlers-cleaners Machine-op-inspct Other-service Priv-house-serv Prof-specialty Protective-serv Sales Tech-support Transport-moving
Female 25.68 0.00 2.21 11.68 0.66 1.68 5.55 17.97 1.38 15.24 0.78 12.76 3.49 0.92
Male 5.93 0.04 18.71 13.98 4.53 5.82 6.98 7.13 0.04 12.50 2.79 11.46 2.80 7.27

Dựa trên Bảng Tần suất chéo giữa giới tính và nghề nghiệp, có thể nhận thấy sự phân hóa nghề nghiệp theo giới tính khá rõ rệt. Các nghề thuộc khối văn phòng và dịch vụ như Adm-clerical (Nhân viên hành chính - văn thư), Other-service (Dịch vụ khác) và Prof-specialty (Chuyên gia chuyên môn) có tỷ lệ nữ giới tham gia cao hơn đáng kể so với nam giới. Đặc biệt, nghề Adm-clerical (Nhân viên hành chính - văn thư) gần như là lĩnh vực đặc trưng của nữ giới với tỷ lệ cao vượt trội.

Ngược lại, các nghề thiên về lao động tay chân hoặc kỹ thuật như Craft-repair (Thợ thủ công - sửa chữa), Farming-fishing (Nông - ngư nghiệp), Handlers-cleaners (Bốc vác - dọn dẹp), Machine-op-inspct (Vận hành máy móc - kiểm tra), và Transport-moving (Vận tải - di chuyển) lại có tỷ lệ nam giới chiếm ưu thế. Ngoài ra, một số nghề nghiệp có tỷ lệ phân bổ giữa nam và nữ tương đối cân bằng như Exec-managerial (Quản lý - điều hành) và Sales (Bán hàng), tuy nhiên nam giới vẫn nhỉnh hơn ở các vị trí quản lý (Exec-managerial: 13.98% nam so với 11.68% nữ).

Đặc biệt, nghề Priv-house-serv (Giúp việc gia đình) (1.38% nữ so với chỉ 0.04% nam) và Protective-serv (Dịch vụ bảo vệ) (2.79% nam so với 0.78% nữ) cho thấy sự phân hóa giới rất rõ do tính chất công việc.

Nhìn chung, bảng số liệu phản ánh rõ xu hướng phân công lao động theo giới tính trong xã hội, trong đó nữ giới tập trung nhiều hơn vào các nghề nghiệp văn phòng, dịch vụ và chuyên môn, trong khi nam giới chiếm ưu thế ở các ngành nghề nặng nhọc hoặc kỹ thuật.

Trực quan hóa

library(dplyr)
library(tidyr)
library(knitr)
library(kableExtra)
library(ggplot2)

# Tạo bảng tần số chéo giữa sex và occupation
freq_table <- table(b$sex, b$occupation)

# Tạo bảng tần suất theo hàng (giới tính)
prop_table <- prop.table(freq_table, margin = 1) * 100

# Chuyển thành dataframe, thêm cột giới tính (rownames_to_column), pivot_longer để ggplot đọc được
df_prop <- as.data.frame.matrix(round(prop_table, 2)) %>%
  tibble::rownames_to_column(var = "Sex") %>%
  pivot_longer(-c(Sex), names_to = "Occupation", values_to = "Percentage")

# Vẽ biểu đồ
ggplot(df_prop, aes(x = Occupation, y = Percentage, fill = Sex)) +
  geom_bar(stat = "identity", position = position_dodge(width = 0.9)) +
  geom_text(aes(label = paste0(round(Percentage, 1), "%")), 
            position = position_dodge(width = 0.9), 
            vjust = -0.3, size = 3) +
  labs(title = "Tỷ lệ nghề nghiệp theo giới tính",
       x = "Nghề nghiệp", y = "Tỷ lệ (%)", fill = "Giới tính") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        legend.position = "top")

Kiểm định thống kê

Một trong những vấn đề được quan tâm trong phân tích thị trường lao động là liệu giới tính có ảnh hưởng đến sự phân bố nghề nghiệp hay không. Để kiểm tra giả thuyết này, Kiểm định Chi-bình phương về tính độc lập được sử dụng nhằm đánh giá:

\[ \left\{ \begin{array}{ll} H_0: & \text{Giới tính và nghề nghiệp của các cá nhân độc lập với nhau.} \\ H_1: & \text{Giới tính và nghề nghiệp của các cá nhân có mối quan hệ với nhau.} \end{array} \right. \]

# Tạo bảng tần số chéo
freq_table <- table(b$sex, b$occupation)

# Thực hiện kiểm định Chi-bình phương
chisq.test(freq_table)
## 
##  Pearson's Chi-squared test
## 
## data:  freq_table
## X-squared = 5716.8, df = 13, p-value < 2.2e-16

Vì p-value < 0.05, bác bỏ giả thuyết H₀. Điều này có nghĩa là Có mối quan hệ có ý nghĩa thống kê giữa giới tính và nghề nghiệp của các cá nhân

So sánh tỷ lệ nam và nữ làm nghề Craft-repair

Bên cạnh việc kiểm định sự độc lập giữa giới tính và nghề nghiệp, việc so sánh tỷ lệ giữa hai nhóm giới tính đối với một nghề nghiệp cụ thể sẽ giúp làm rõ sự khác biệt cụ thể trong phân bố nghề nghiệp. Trong trường hợp này, chúng ta tiến hành kiểm định hiệu hai tỷ lệ nhằm so sánh tỷ lệ nam và nữ làm nghề Thủ công, sửa chữa (Craft-repair).

\[ \left\{ \begin{array}{ll} H_0: & p_{Female} = p_{Male} \quad (\text{Tỷ lệ nữ làm nghề cơ khí bằng tỷ lệ nam}) \\ H_1: & p_{Female} \neq p_{Male} \quad (\text{Tỷ lệ nữ làm nghề khác tỷ lệ nam}) \end{array} \right. \]

# Bảng tần suất giữa giới tính và nghề Craft-repair
table_sex_craft <- table(b$sex, b$occupation == "Craft-repair")
colnames(table_sex_craft) <- c("Other occupations", "Craft-repair")
print(table_sex_craft)
##         
##          Other occupations Craft-repair
##   Female              9566          216
##   Male               16566         3814
# Tính tỷ lệ làm nghề Craft-repair theo giới tính
prop.table(table_sex_craft, margin = 1)
##         
##          Other occupations Craft-repair
##   Female        0.97791863   0.02208137
##   Male          0.81285574   0.18714426
# prop.test để kiểm định hiệu tỷ lệ
prop.test(x = table_sex_craft[, "Craft-repair"], n = rowSums(table_sex_craft), correct = FALSE)
## 
##  2-sample test for equality of proportions without continuity correction
## 
## data:  table_sex_craft[, "Craft-repair"] out of rowSums(table_sex_craft)
## X-squared = 1555.7, df = 1, p-value < 2.2e-16
## alternative hypothesis: two.sided
## 95 percent confidence interval:
##  -0.1711583 -0.1589675
## sample estimates:
##     prop 1     prop 2 
## 0.02208137 0.18714426

Kết quả kiểm định chi bình phương với hiệu hai tỷ lệ cho giá trị X-squared = 1555,7 với p-value < 2.2e-16, nhỏ hơn mức ý nghĩa 5%, do đó bác bỏ giả thuyết H₀. Điều này cho thấy giới tính có ảnh hưởng đáng kể đến khả năng làm việc trong lĩnh vực cơ khí. Nói cách khác, nam giới có xu hướng làm nghề thủ công, sửa chữa cao hơn đáng kể so với nữ giới trong bộ dữ liệu nghiên cứu.

Relative Risk (RR)

Sau khi tiến hành kiểm định hiệu hai tỷ lệ, để lượng hóa mức độ chênh lệch giữa hai nhóm giới tính, chúng ta tiếp tục sử dụng Relative Risk - RR. Chỉ số này cho phép đánh giá xác suất làm nghề Thủ công, sửa chữa ở nam giới so với nữ giới.

library(epitools)

# 1. Tạo bảng 2x2 với số liệu cụ thể (đảm bảo đúng thứ tự: Other occupations - Craft-repair)
table_sex_craft <- table(b$sex, b$occupation == "Craft-repair")

# 2. Tạo lại bảng 2x2 theo đúng thứ tự mong muốn: hàng là giới tính, cột là nghề nghiệp
table_2x2 <- matrix(c(
  table_sex_craft["Female", "FALSE"], table_sex_craft["Female", "TRUE"],
  table_sex_craft["Male", "FALSE"], table_sex_craft["Male", "TRUE"]
), nrow = 2, byrow = TRUE)

# 3. Đặt tên hàng và cột cho dễ hiểu
dimnames(table_2x2) <- list(
  Gender = c("Female", "Male"),
  Occupation = c("Other occupations", "Craft-repair")
)

# 4. In bảng tần suất để kiểm tra
print(table_2x2)
##         Occupation
## Gender   Other occupations Craft-repair
##   Female              9566          216
##   Male               16566         3814
# 5. Tính Relative Risk (RR) với khoảng tin cậy 95%
result_rr <- riskratio(table_2x2, method = "wald", conf.level = 0.95)

# 6. In kết quả RR
print(result_rr)
## $data
##         Occupation
## Gender   Other occupations Craft-repair Total
##   Female              9566          216  9782
##   Male               16566         3814 20380
##   Total              26132         4030 30162
## 
## $measure
##         risk ratio with 95% C.I.
## Gender   estimate    lower    upper
##   Female 1.000000       NA       NA
##   Male   8.475209 7.405322 9.699669
## 
## $p.value
##         two-sided
## Gender   midp.exact fisher.exact chi.square
##   Female         NA           NA         NA
##   Male            0            0          0
## 
## $correction
## [1] FALSE
## 
## attr(,"method")
## [1] "Unconditional MLE & normal approximation (Wald) CI"

Kết quả phân tích Risk Ratio cho thấy sự khác biệt rất rõ về khả năng tham gia vào nghề “Thủ công và sửa chữa” giữa hai nhóm giới tính.

Giá trị Risk Ratio = 8.48 cho biết khả năng làm nghề “Thủ công và sửa chữa” của nam giới cao hơn khoảng 8.48 lần so với nữ giới.

Nói cách khác, nam giới có xu hướng làm nghề “Thủ công và sửa chữa” cao hơn đáng kể so với nữ giới.

Nhìn chung, giới tính có ảnh hưởng đáng kể đến khả năng tham gia vào nghề “Thủ công và sửa chữa”, với nam giới có tỷ lệ tham gia cao vượt trội so với nữ giới. Điều này phù hợp với đặc điểm của ngành nghề này, vốn đòi hỏi kỹ năng lao động thủ công, cơ khí, thường được nam giới đảm nhiệm nhiều hơn trong thực tế.

Odds và tỷ số Chênh (Odds Ratio - OR)

# Đặt mã hóa UTF-8 để hỗ trợ tiếng Việt
Sys.setlocale("LC_ALL", "en_US.UTF-8")
## [1] "LC_COLLATE=en_US.UTF-8;LC_CTYPE=en_US.UTF-8;LC_MONETARY=en_US.UTF-8;LC_NUMERIC=C;LC_TIME=en_US.UTF-8"
# Tạo bảng tần suất giới tính × nghề Craft-repair (cột "Craft-repair" ở sau)
table_sex_craft <- table(b$sex, b$occupation == "Craft-repair")
colnames(table_sex_craft) <- c("Other occupations", "Craft-repair")

# Thêm cột Total để dễ quan sát
table_sex_craft_with_total <- cbind(table_sex_craft, Total = rowSums(table_sex_craft))
print(table_sex_craft_with_total)
##        Other occupations Craft-repair Total
## Female              9566          216  9782
## Male               16566         3814 20380
# Tính odds cho từng nhóm
odds_female <- table_sex_craft["Female", "Craft-repair"] / table_sex_craft["Female", "Other occupations"]
odds_male <- table_sex_craft["Male", "Craft-repair"] / table_sex_craft["Male", "Other occupations"]

# In kết quả Odds từng nhóm
cat("Odds nhóm Female:", round(odds_female, 4), "\n")
## Odds nhóm Female: 0.0226
cat("Odds nhóm Male:", round(odds_male, 4), "\n")
## Odds nhóm Male: 0.2302

Kết quả phân tích tỷ số odds (Odds) cho thấy sự khác biệt rất lớn giữa nam và nữ trong khả năng làm nghề Craft-repair:

Odds (nữ) = 0,0226, trong nhóm nữ, cứ 1 người làm nghề khác thì có khoảng 0.023 người làm nghề “Craft-repair”. Nói cách khác, Trung bình cứ 1 nữ giới làm nghề “Craft-repair” thì có khoảng 44.25 nữ giới làm nghề khác.

Odds (nam) = 0,2302, nghĩa là trong nhóm nam, cứ 1 người làm nghề khác thì có khoảng 0.23 người làm nghề “Craft-repair”. Nói cách khác, Trung bình cứ 1 nam giới làm nghề “Craft-repair” thì có khoảng 4.34 nam giới làm nghề khác.

Odds Ratio (OR)

library(epitools)

# 1. Tạo bảng tần suất giới tính × nghề Craft-repair (cột "Craft-repair" ở sau cho dễ đọc)
table_sex_craft <- table(b$sex, b$occupation == "Craft-repair")
colnames(table_sex_craft) <- c("Other occupations", "Craft-repair")

# 2. Thêm cột Total để dễ quan sát
table_sex_craft_with_total <- cbind(table_sex_craft, Total = rowSums(table_sex_craft))
print(table_sex_craft_with_total)
##        Other occupations Craft-repair Total
## Female              9566          216  9782
## Male               16566         3814 20380
# 3. Tính Odds Ratio (tỷ số odds) bằng hàm oddsratio(), method = "wald", khoảng tin cậy 95%
result_or <- oddsratio(table_sex_craft, method = "wald", conf.level = 0.95)

# 4. In kết quả Odds Ratio
print(result_or)
## $data
##         
##          Other occupations Craft-repair Total
##   Female              9566          216  9782
##   Male               16566         3814 20380
##   Total              26132         4030 30162
## 
## $measure
##         odds ratio with 95% C.I.
##          estimate    lower    upper
##   Female  1.00000       NA       NA
##   Male   10.19623 8.869725 11.72112
## 
## $p.value
##         two-sided
##          midp.exact fisher.exact chi.square
##   Female         NA           NA         NA
##   Male            0            0          0
## 
## $correction
## [1] FALSE
## 
## attr(,"method")
## [1] "Unconditional MLE & normal approximation (Wald) CI"

Tỷ số odds (Odds Ratio) giữa nam và nữ khi làm nghề “Craft-repair” là 10.20.Điều này có nghĩa là tỷ lệ giữa số nam làm nghề “Craft-repair” so với số nam làm nghề khác cao gấp khoảng 10.2 lần so với tỷ lệ đó ở nữ.

Phần 3: Hồi quy Logistic

1. Giới thiệu

Hồi quy logistic là một phương pháp thống kê dùng để dự đoán xác suất xảy ra của một sự kiện nhị phân (có/không, đúng/sai, thành công/thất bại, v.v.). Mô hình này rất hữu ích khi biến kết quả (biến phụ thuộc) chỉ có 2 giá trị, ví dụ:

  • Một khách hàng có mua hàng không?
  • Một người có mắc bệnh không?
  • Một cá nhân có thu nhập > 50K không?

2. Cấu trúc của mô hình

Ta gọi biến kết quả là \(Y \in \{0,1\}\), trong đó:

  • \(Y = 1\): sự kiện xảy ra (ví dụ: có bệnh)
  • \(Y = 0\): sự kiện không xảy ra

Mục tiêu của mô hình là dự đoán xác suất \(P(Y = 1 \mid X)\) dựa trên các biến độc lập \(X_1, X_2, \dots, X_k\).


3. Hàm logistic

Ta không mô hình hóa trực tiếp xác suất \(P\), mà dùng một hàm chuyển đổi gọi là hàm logistic (hay sigmoid):

\[ P(Y = 1 \mid X) = \frac{1}{1 + e^{-(\beta_0 + \beta_1 X_1 + \cdots + \beta_k X_k)}} \]

Trong đó:

  • \(\beta_0\): hệ số chặn (intercept)
  • \(\beta_j\): hệ số của biến \(X_j\)
  • Hàm logistic luôn cho ra kết quả nằm giữa 0 và 1.

4. Biến đổi logit

Ta cũng có thể viết mô hình dưới dạng logit – log của odds:

\[ \text{logit}(P) = \log\left(\frac{P(Y=1 \mid X)}{1 - P(Y=1 \mid X)}\right) = \beta_0 + \beta_1 X_1 + \cdots + \beta_k X_k \]

Ở đây:

  • Odds = xác suất xảy ra / xác suất không xảy ra.
  • Logit là một phép biến đổi giúp đưa mô hình về dạng tuyến tính.

5. Ý nghĩa hệ số hồi quy

Mỗi hệ số \(\beta_j\) cho biết ảnh hưởng của biến \(X_j\) lên log-odds của việc xảy ra sự kiện.

  • Nếu \(\beta_j > 0\): biến \(X_j\) làm tăng xác suất xảy ra sự kiện.
  • Nếu \(\beta_j < 0\): biến \(X_j\) làm giảm xác suất.
  • Ta thường chuyển sang tỷ số odds (odds ratio) để diễn giải dễ hiểu hơn:

\[ \text{Odds Ratio} = e^{\beta_j} \]

Ví dụ: nếu \(\beta_j = 0.7\), thì odds tăng khoảng \(e^{0.7} \approx 2\), tức là gấp đôi.


6. Ước lượng tham số – Maximum Likelihood

Các hệ số \(\beta\) được ước lượng bằng phương pháp tối đa hóa hàm hợp lý (Maximum Likelihood Estimation – MLE).

  • Ta viết hàm xác suất (likelihood function):

\[ L(\beta) = \prod_{i=1}^{n} \left[ \pi(x_i)^{y_i} (1 - \pi(x_i))^{1 - y_i} \right] \]

  • Với:

\[ \pi(x_i) = \frac{1}{1 + e^{-(\beta_0 + \beta_1 x_{i1} + \cdots + \beta_k x_{ik})}} \]

  • Rồi lấy log để đơn giản hóa:

\[ \ell(\beta) = \sum_{i=1}^{n} \left[ y_i \log(\pi(x_i)) + (1 - y_i) \log(1 - \pi(x_i)) \right] \]

Sau đó, tìm bộ tham số \(\beta\) làm tối đa hóa log-likelihood.

7. Phân tích tác động của giới tính đến xác suất thu nhập cao bằng hồi quy logistic

Trong phần này, em sử dụng hồi quy logistic nhị phân để phân tích tác động của giới tính (sex) đến khả năng có thu nhập cao (>50K) của các cá nhân trong bộ dữ liệu adultcensusincome.

Do biến thu nhập mang tính nhị phân (“>50K” và “<=50K”), mô hình hồi quy logistic là phương pháp thích hợp để ước lượng xác suất một cá nhân có thu nhập cao, và kiểm tra xem giới tính có ảnh hưởng đáng kể đến xác suất này hay không.

Cụ thể:

  • Biến phụ thuộc: income_bin – được tạo từ biến income, với: \[ income\_bin = \begin{cases} 1 & \text{nếu thu nhập } > 50K \\ 0 & \text{nếu thu nhập } \leq 50K \end{cases} \]

  • Biến độc lập: sex – biến định tính gồm hai mức: "Male""Female" (được mã hóa là factor, trong đó "Female" là mức tham chiếu mặc định).

  • Mô hình hồi quy logistic được xây dựng như sau:

\[ \text{logit}(P(Y = 1 \mid \text{sex})) = \beta_0 + \beta_1 \cdot \text{sex} \]

Trong đó:

  • \(Y = 1\): cá nhân có thu nhập >50K
  • \(\beta_0\): log-odds thu nhập cao đối với nhóm "Female" (mức tham chiếu).
  • \(\beta_1\): phần chênh lệch log-odds giữa "Male""Female"

Tỷ số odds được tính bằng:

\[ \text{Odds Ratio} = e^{\beta_1} \]

Cho biết tỷ lệ odds của nam giới có thu nhập >50K so với nữ giới.

# Đọc dữ liệu từ file CSV
df <- read.csv("C:/Users/Admin/Desktop/adultcensusincome.csv", 
               header = TRUE, sep = ",")

# Tải các thư viện cần thiết
library(ggplot2)
library(ggthemes)
library(scales)

# BƯỚC 1: Tạo biến nhị phân từ income (1 nếu >50K, 0 nếu <=50K)
df$income_bin <- ifelse(trimws(df$income) == ">50K", 1, 0)

# BƯỚC 2: Chuyển biến sex thành factor
df$sex <- as.factor(df$sex)

# BƯỚC 3: Xây dựng mô hình hồi quy logistic với biến sex
model_sex <- glm(income_bin ~ sex, data = df, family = binomial)

# BƯỚC 4: Tóm tắt kết quả mô hình
summary(model_sex)
## 
## Call:
## glm(formula = income_bin ~ sex, family = binomial, data = df)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -2.05371    0.03185  -64.47   <2e-16 ***
## sexMale      1.27147    0.03525   36.07   <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: 33851  on 30161  degrees of freedom
## Residual deviance: 32287  on 30160  degrees of freedom
## AIC: 32291
## 
## Number of Fisher Scoring iterations: 4
# BƯỚC 5: Tính Odds Ratio và khoảng tin cậy 95%
odds_ratio <- exp(coef(model_sex))
ci <- exp(confint.default(model_sex))

cat("Kết quả Odds Ratio:\n")
## Kết quả Odds Ratio:
print(round(cbind(Odds_Ratio = odds_ratio, CI_Lower = ci[,1], CI_Upper = ci[,2]), 3))
##             Odds_Ratio CI_Lower CI_Upper
## (Intercept)      0.128    0.120    0.137
## sexMale          3.566    3.328    3.821
# BƯỚC 4: Dự đoán xác suất có thu nhập >50K
df$predicted_prob <- predict(model_sex, type = "response")

# BƯỚC 5: Tính xác suất trung bình theo giới tính
prob_by_sex <- aggregate(predicted_prob ~ sex, data = df, FUN = mean)

Phương trình hồi quy logistic:

\[ \text{logit}(P) = \log\left( \frac{P(Y = 1 \mid \text{sex})}{1 - P(Y = 1 \mid \text{sex})} \right) = -2.0537 + 1.2715 \cdot \text{sex}_{\text{Male}} \] Trong đó:

  • \(Y = 1\) nếu cá nhân có thu nhập >50K, \(Y = 0\) nếu không.
  • \(\text{sex}_{\text{Male}} = 1\) nếu cá nhân là nam, và bằng 0 nếu là nữ.

Kết quả hồi quy cho thấy hệ số của nam giới là 1.271, có ý nghĩa thống kê rất cao (p < 0.001). Hệ số này cho biết nam giới có log-odds thu nhập cao hơn đáng kể so với nữ. Khi chuyển sang odds ratio, kết quả cho thấy nam giới có khả năng đạt được thu nhập >50K cao gấp 3.57 lần nữ giới, trong điều kiện không xét đến các yếu tố khác.

Mức deviance giảm từ 33.851 xuống 32.287 sau khi đưa biến giới tính vào mô hình, cho thấy biến này giúp cải thiện đáng kể độ phù hợp của mô hình.

Xác suất dự đoán:

Từ phương trình trên, ta có thể tính được xác suất dự đoán bằng cách áp dụng hàm logistic:

# Tính xác suất dự đoán từ mô hình đã ước lượng
df$predicted_prob <- predict(model_sex, type = "response")

# Tính xác suất trung bình theo giới tính
library(dplyr)
prob_by_sex <- df %>%
  group_by(sex) %>%
  summarise(mean_prob = mean(predicted_prob))

print(prob_by_sex)
## # A tibble: 2 × 2
##   sex    mean_prob
##   <fct>      <dbl>
## 1 Female     0.114
## 2 Male       0.314

Kết quả cho thấy xác suất trung bình một cá nhân có thu nhập cao (>50K) là khoảng 31.4% đối với nam giới, trong khi chỉ khoảng 11.4% đối với nữ giới. Điều này cho thấy sự khác biệt rõ rệt về khả năng đạt mức thu nhập cao giữa hai giới, với nam giới có xác suất cao hơn đáng kể so với nữ giới.

\[ P(Y = 1 \mid \text{sex}) = \frac{1}{1 + \exp\left[ -(-2.0537 + 1.2715 \cdot \text{sex}_{\text{Male}}) \right]} \]

Cụ thể:

  • Với nữ giới (\(\text{sex}_{\text{Male}} = 0\)):

\[ P(Y = 1 \mid \text{Female}) = \frac{1}{1 + e^{2.0537}} \approx 0.113 \]

  • Với nam giới (\(\text{sex}_{\text{Male}} = 1\)):

\[ P(Y = 1 \mid \text{Male}) = \frac{1}{1 + e^{-(-2.0537 + 1.2715)}} \approx 0.288 \]

# Vẽ biểu đồ xác suất trung bình thu nhập >50K theo giới tính
ggplot(prob_by_sex, aes(x = sex, y = mean_prob, fill = sex)) +
  geom_col(width = 0.6, color = "black", show.legend = FALSE) +
  geom_text(
    aes(label = percent(mean_prob, accuracy = 0.1)),
    vjust = -0.5, size = 6, fontface = "bold"
  ) +
  scale_y_continuous(labels = percent_format(accuracy = 1), limits = c(0, 0.6)) +
  scale_fill_manual(values = c("Female" = "#FF9999", "Male" = "#6699CC")) +
  labs(
    title = "Xác suất trung bình có thu nhập >50K theo giới tính",
    x = "Giới tính",
    y = "Xác suất dự đoán"
  ) +
  theme_minimal() +
  theme(
    plot.title = element_text(hjust = 0.5, size = 18, face = "bold"),
    axis.title = element_text(size = 14),
    axis.text = element_text(size = 13),
    panel.grid.major.y = element_line(color = "gray85")
  )

8. Phân tích tác động của giới tính, học vấn, nghề nghiệp và tình trạng hôn nhân đến xác suất một cá nhân có thu nhập cao (>50K) bằng mô hình hồi quy logistic

Trong phần phân tích này, em xây dựng một mô hình hồi quy logistic nhị phân nhằm ước lượng xác suất một cá nhân có thu nhập >50K dựa trên các đặc điểm nhân khẩu học quan trọng bao gồm:

  • Giới tính (nam hay nữ)

  • Nhóm nghề nghiệp (làm quản lý điều hành hay không)

  • Trình độ học vấn (trình độ cao như đại học, sau đại học hay thấp hơn)

  • Tình trạng hôn nhân (đã kết hôn hay chưa)

  • Nhóm tham chiếu được chọn là:

    • Nữ giới (Female)
    • Trình độ học vấn cao (Higher education: Bachelors, Masters, Prof-school, Doctorate)
    • Đã kết hôn (Married: Married-civ-spouse, Married-spouse-absent, Married-AF-spouse )
    • Nghề quản lý/điều hành (Exec-managerial)

Mô hình hồi quy logistic được thiết lập để kiểm tra xem việc thay đổi mỗi đặc điểm trên có làm thay đổi đáng kể xác suất một người đạt mức thu nhập cao hay không, trong khi giữ các yếu tố còn lại cố định.

Về mặt toán học, mô hình có dạng:

\[ \log\left( \frac{P(Y = 1)}{1 - P(Y = 1)} \right) = \beta_0 + \beta_1 \cdot \text{sex}_{\text{Male}} + \beta_2 \cdot \text{education\_group}_{\text{Lower}} + \beta_3 \cdot \text{marital\_group}_{\text{Single}} + \beta_4 \cdot \text{occupation\_group}_{\text{Other}} \]

  • \(Y = 1\) nghĩa là cá nhân có thu nhập lớn hơn 50K.
  • \(\beta_1\): log-odds chênh lệch giữa nam và nữ (trong đó nữ là nhóm tham chiếu).
  • \(\beta_2\): log-odds chênh lệch giữa nhóm có trình độ học vấn thấp (Lower) so với nhóm học vấn cao (Higher — tham chiếu).
  • \(\beta_3\): log-odds chênh lệch giữa người chưa kết hôn (Single) so với người đã kết hôn (Married — tham chiếu).
  • \(\beta_4\): log-odds chênh lệch giữa nhóm nghề nghiệp không phải quản lý/điều hành (Other) so với nhóm quản lý/điều hành (Exec-managerial — tham chiếu).
library(dplyr)

# BƯỚC 1: Loại bỏ các dòng có giá trị thiếu trong các cột liên quan
df <- df %>%
  filter(!is.na(income) & !is.na(marital.status) & !is.na(occupation) &
         !is.na(sex) & !is.na(education))

# BƯỚC 2: Tạo biến nhóm trình độ học vấn (Higher vs Lower)
df$education_group <- ifelse(df$education %in% c("Bachelors", "Masters", "Prof-school", "Doctorate"), 
                             "Higher", "Lower")

# BƯỚC 3: Tạo biến nhóm tình trạng hôn nhân (Married vs Single)
df$marital_group <- ifelse(df$marital.status %in% c("Married-civ-spouse", "Married-spouse-absent", "Married-AF-spouse"), 
                           "Married", "Single")

# BƯỚC 4: Tạo biến nhóm nghề nghiệp (Exec-managerial vs Other)
df$occupation_group <- ifelse(df$occupation == "Exec-managerial", 
                              "Exec-managerial", "Other")

# BƯỚC 5: Biến nhị phân từ income (1 nếu >50K, 0 nếu <=50K)
df$income_binary <- ifelse(trimws(df$income) == ">50K", 1, 0)

# BƯỚC 6: Chuyển thành factor và đặt nhóm tham chiếu
df$sex <- factor(df$sex)
df$sex <- relevel(df$sex, ref = "Female")  # Đổi nhóm tham chiếu thành Female

df$education_group <- factor(df$education_group)
df$education_group <- relevel(df$education_group, ref = "Higher")

df$marital_group <- factor(df$marital_group)
df$marital_group <- relevel(df$marital_group, ref = "Married")

df$occupation_group <- factor(df$occupation_group)
df$occupation_group <- relevel(df$occupation_group, ref = "Exec-managerial")

# BƯỚC 7: Hồi quy logistic
model <- glm(income_binary ~ sex + occupation_group + education_group + marital_group, 
             family = binomial(link = "logit"), 
             data = df)

# BƯỚC 8: Kết quả hồi quy
summary(model)
## 
## Call:
## glm(formula = income_binary ~ sex + occupation_group + education_group + 
##     marital_group, family = binomial(link = "logit"), data = df)
## 
## Coefficients:
##                       Estimate Std. Error z value Pr(>|z|)    
## (Intercept)            1.39101    0.05925  23.479  < 2e-16 ***
## sexMale                0.33522    0.04361   7.686 1.51e-14 ***
## occupation_groupOther -0.90652    0.04300 -21.082  < 2e-16 ***
## education_groupLower  -1.62602    0.03509 -46.340  < 2e-16 ***
## marital_groupSingle   -2.38010    0.04144 -57.441  < 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: 33851  on 30161  degrees of freedom
## Residual deviance: 24418  on 30157  degrees of freedom
## AIC: 24428
## 
## Number of Fisher Scoring iterations: 5

Phương trình hồi quy logistic:

\[ \log\left( \frac{P(Y = 1)}{1 - P(Y = 1)} \right) = 1.391 + 0.335 \cdot \text{sex}_{\text{Male}} - 0.907 \cdot \text{occupation}_{\text{Other}} - 1.626 \cdot \text{education}_{\text{Lower}} - 2.380 \cdot \text{marital}_{\text{Single}} \]

Kết quả hồi quy cho thấy tất cả các biến đưa vào mô hình đều có ý nghĩa thống kê rất cao (p < 0.001), điều này phản ánh rằng chúng đều có ảnh hưởng rõ rệt đến xác suất đạt được mức thu nhập cao.

Cụ thể, hệ số của biến sexMale là 0.335. Điều này cho thấy nam giới có log-odds thu nhập cao lớn hơn nữ giới, khi giữ các yếu tố khác không đổi. Khi chuyển sang xác suất, kết quả tương đương với việc nam giới có khả năng đạt thu nhập cao cao hơn khoảng 40% so với nữ giới.

Tiếp theo, hệ số của biến occupation_groupOther là -0.907. Điều này phản ánh rằng những người không làm trong nhóm nghề “quản lý điều hành” có log-odds thu nhập cao thấp hơn đáng kể. Cụ thể, họ có khả năng đạt mức thu nhập >50K thấp hơn khoảng 59% so với nhóm nghề quản lý.

Tương tự, biến education_groupLower có hệ số là -1.626, cho thấy những người có trình độ học vấn thấp hơn có khả năng đạt được mức thu nhập cao thấp hơn khoảng 80% so với những người có trình độ học vấn cao.

Cuối cùng, biến marital_groupSingle có hệ số -2.380, là mức ảnh hưởng mạnh nhất trong mô hình. Điều này cho thấy những người chưa kết hôn có khả năng đạt được mức thu nhập >50K thấp hơn khoảng 91% so với những người đã lập gia đình.

Tổng kết lại, mô hình hồi quy logistic này khẳng định rằng các yếu tố xã hội như giới tính, học vấn, nghề nghiệp và tình trạng hôn nhân đều có tác động đáng kể đến xác suất thu nhập cao. Đặc biệt, học vấn và hôn nhân là hai yếu tố có ảnh hưởng mạnh mẽ nhất trong mô hình.

# Tạo bộ dữ liệu mới với các nhóm cụ thể muốn tính xác suất
new_data <- data.frame(
  sex = factor(c("Female", "Male"), levels = levels(df$sex)),
  occupation_group = factor(rep("Exec-managerial", 2), levels = levels(df$occupation_group)),
  education_group = factor(rep("Higher", 2), levels = levels(df$education_group)),
  marital_group = factor(rep("Married", 2), levels = levels(df$marital_group))
)

# Tính xác suất dự đoán từ mô hình logistic đã huấn luyện
predictions <- predict(model, newdata = new_data, type = "response")

# Gộp kết quả
result <- cbind(new_data, Probability = round(predictions, 4))

# Hiển thị bảng đẹp
knitr::kable(result, caption = "Bảng: Xác suất dự đoán thu nhập >50K theo giới tính (giữ các yếu tố khác cố định)", 
             align = "c") |>
  kableExtra::kable_styling(full_width = FALSE, bootstrap_options = c("striped", "hover", "condensed"))
Bảng: Xác suất dự đoán thu nhập >50K theo giới tính (giữ các yếu tố khác cố định)
sex occupation_group education_group marital_group Probability
Female Exec-managerial Higher Married 0.8008
Male Exec-managerial Higher Married 0.8489

\[ P(Y = 1 \mid X) = \frac{1}{1 + \exp\left[ -(1.391 + 0.335 \cdot \text{sex}_{\text{Male}} - 0.907 \cdot \text{occupation}_{\text{Other}} - 1.626 \cdot \text{education}_{\text{Lower}} - 2.380 \cdot \text{marital}_{\text{Single}}) \right]} \]

Cụ thể, với nữ giới, xác suất dự đoán đạt thu nhập cao là khoảng 80.08%, trong khi với nam giới có cùng đặc điểm, xác suất này là khoảng 84.89%. Điều này cho thấy, trong cùng một điều kiện lý tưởng về nghề nghiệp, học vấn và hôn nhân, nam giới vẫn có khả năng đạt thu nhập cao cao hơn khoảng 4.8% so với nữ giới.

library(ggplot2)
library(scales)

# Vẽ biểu đồ xác suất dự đoán theo giới tính
ggplot(result, aes(x = sex, y = Probability, fill = sex)) +
  geom_col(width = 0.6, color = "black", show.legend = FALSE) +
  geom_text(aes(label = percent(Probability, accuracy = 0.1)), 
            vjust = -0.5, size = 6, fontface = "bold") +
  scale_y_continuous(labels = percent_format(accuracy = 1), limits = c(0, 1)) +
  scale_fill_manual(values = c("Female" = "#FF9999", "Male" = "#6699CC")) +
  labs(
    title = "Xác suất dự đoán thu nhập >50K theo giới tính",
    subtitle = "Giữ cố định nghề nghiệp: Exec-managerial, học vấn: cao, hôn nhân: đã kết hôn",
    x = "Giới tính",
    y = "Xác suất dự đoán"
  ) +
  theme_minimal() +
  theme(
    plot.title = element_text(hjust = 0.5, size = 18, face = "bold"),
    plot.subtitle = element_text(hjust = 0.5, size = 13, face = "italic"),
    axis.title = element_text(size = 14),
    axis.text = element_text(size = 12),
    panel.grid.major.y = element_line(color = "gray85")
  )

9.Ảnh hưởng của trình độ học vấn đến xác suất đạt được thu nhập cao

Trong phần phân tích này, em sử dụng hồi quy logistic để kiểm tra xem trình độ học vấn của một cá nhân có ảnh hưởng đáng kể đến xác suất có thu nhập >50K hay không. Đây là một câu hỏi kinh điển trong nghiên cứu về bất bình đẳng thu nhập và kinh tế lao động.

Do là biến định tính có nhiều mức, mô hình hồi quy logistic sẽ so sánh từng mức học vấn với một nhóm tham chiếu là Doctorate.

Mô hình có dạng:

\[ \log\left(\frac{P(Y = 1)}{1 - P(Y = 1)}\right) = \beta_0 + \beta_1 \cdot \text{Edu}_1 + \beta_2 \cdot \text{Edu}_2 + \cdots + \beta_k \cdot \text{Edu}_k \]

Trong đó: - \(Y = 1\) tương ứng với thu nhập >50K. - Mỗi \(\beta_j\) đo lường log-odds tăng thêm khi trình độ học vấn là mức \(j\), so với nhóm tham chiếu. cao nổi bật.

# BƯỚC 1: Tạo biến nhị phân từ income (1 = >50K, 0 = <=50K)
df$income_bin <- ifelse(trimws(df$income) == ">50K", 1, 0)

# BƯỚC 2: Chuyển biến education thành factor
df$education <- as.factor(df$education)

# BƯỚC 3: Đặt "Doctorate" làm nhóm tham chiếu
df$education <- relevel(df$education, ref = "Doctorate")

# BƯỚC 4: Xây dựng mô hình hồi quy logistic với biến education
model_edu <- glm(income_bin ~ education, data = df, family = binomial)

# BƯỚC 5: Xem kết quả mô hình
summary(model_edu)
## 
## Call:
## glm(formula = income_bin ~ education, family = binomial, data = df)
## 
## Coefficients:
##                        Estimate Std. Error z value Pr(>|z|)    
## (Intercept)             1.08091    0.11873   9.104  < 2e-16 ***
## education10th          -3.63801    0.17989 -20.223  < 2e-16 ***
## education11th          -3.90007    0.17905 -21.782  < 2e-16 ***
## education12th          -3.56582    0.22683 -15.720  < 2e-16 ***
## education1st-4th       -4.26589    0.43320  -9.847  < 2e-16 ***
## education5th-6th       -4.21641    0.31789 -13.264  < 2e-16 ***
## education7th-8th       -3.78323    0.21115 -17.917  < 2e-16 ***
## education9th           -3.92582    0.23754 -16.527  < 2e-16 ***
## educationAssoc-acdm    -2.15847    0.13905 -15.523  < 2e-16 ***
## educationAssoc-voc     -2.11032    0.13432 -15.711  < 2e-16 ***
## educationBachelors     -1.39757    0.12211 -11.445  < 2e-16 ***
## educationHS-grad       -2.70728    0.12181 -22.225  < 2e-16 ***
## educationMasters       -0.82257    0.12883  -6.385 1.72e-10 ***
## educationPreschool    -14.64698   79.81447  -0.184    0.854    
## educationProf-school    0.01279    0.15464   0.083    0.934    
## educationSome-college  -2.46683    0.12261 -20.119  < 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: 33851  on 30161  degrees of freedom
## Residual deviance: 29946  on 30146  degrees of freedom
## AIC: 29978
## 
## Number of Fisher Scoring iterations: 12
# BƯỚC 6: Tính Odds Ratio và khoảng tin cậy 95%
odds_ratio <- exp(coef(model_edu))
ci <- exp(confint.default(model_edu))

# In bảng Odds Ratio
result_table <- round(cbind(Odds_Ratio = odds_ratio,
                            CI_lower = ci[,1],
                            CI_upper = ci[,2]), 3)
print(result_table)
##                       Odds_Ratio CI_lower     CI_upper
## (Intercept)                2.947    2.335 3.720000e+00
## education10th              0.026    0.018 3.700000e-02
## education11th              0.020    0.014 2.900000e-02
## education12th              0.028    0.018 4.400000e-02
## education1st-4th           0.014    0.006 3.300000e-02
## education5th-6th           0.015    0.008 2.800000e-02
## education7th-8th           0.023    0.015 3.400000e-02
## education9th               0.020    0.012 3.100000e-02
## educationAssoc-acdm        0.116    0.088 1.520000e-01
## educationAssoc-voc         0.121    0.093 1.580000e-01
## educationBachelors         0.247    0.195 3.140000e-01
## educationHS-grad           0.067    0.053 8.500000e-02
## educationMasters           0.439    0.341 5.650000e-01
## educationPreschool         0.000    0.000 3.776594e+61
## educationProf-school       1.013    0.748 1.371000e+00
## educationSome-college      0.085    0.067 1.080000e-01

Mô hình hồi quy logistic sử dụng để đánh giá mối liên hệ giữa trình độ học vấn và xác suất đạt thu nhập trên 50K cho thấy sự khác biệt rõ rệt giữa các nhóm giáo dục. Trong mô hình này, nhóm tham chiếu là những người có trình độ Doctorate, do đó các hệ số hồi quy phản ánh sự thay đổi log-odds thu nhập cao so với nhóm có học vị tiến sĩ.

Kết quả cho thấy hầu hết các nhóm học vấn thấp hơn đều có hệ số âm, cho thấy log-odds đạt thu nhập >50K giảm đáng kể so với nhóm Doctorate. Cụ thể, những người chỉ tốt nghiệp phổ thông trung học (HS-grad) có hệ số hồi quy là -2.707, tương ứng với odds ratio là 0.067, nghĩa là xác suất họ có thu nhập cao chỉ bằng khoảng 6.7% so với người có bằng tiến sĩ, nếu không xét đến các yếu tố khác. Tương tự, nhóm có trình độ đại học (Bachelors) có odds ratio là 0.247, cho thấy khả năng đạt thu nhập cao chỉ bằng khoảng 24.7% so với nhóm Doctorate.

Mặt khác, một số nhóm như Prof-school (trường chuyên nghiệp sau đại học) có hệ số không đáng kể về mặt thống kê (p = 0.934), cho thấy không có sự khác biệt rõ ràng so với nhóm Doctorate. Các nhóm học vấn rất thấp như “1st-4th grade” hoặc “5th-6th grade” có odds ratio cực kỳ nhỏ (khoảng 0.01–0.03), phản ánh xác suất thu nhập cao gần như không đáng kể.

# BƯỚC 7: Dự đoán xác suất thu nhập >50K cho từng cá nhân
df$predicted_prob <- predict(model_edu, type = "response")

# BƯỚC 9: Tính xác suất trung bình theo từng nhóm học vấn
mean_probs <- aggregate(predicted_prob ~ education, data = df, FUN = mean)
print(mean_probs)
##       education predicted_prob
## 1     Doctorate   7.466667e-01
## 2          10th   7.195122e-02
## 3          11th   5.629771e-02
## 4          12th   7.692308e-02
## 5       1st-4th   3.973510e-02
## 6       5th-6th   4.166667e-02
## 7       7th-8th   6.283662e-02
## 8           9th   5.494505e-02
## 9    Assoc-acdm   2.539683e-01
## 10    Assoc-voc   2.631982e-01
## 11    Bachelors   4.214909e-01
## 12      HS-grad   1.643293e-01
## 13      Masters   5.642286e-01
## 14    Preschool   1.283310e-06
## 15  Prof-school   7.490775e-01
## 16 Some-college   2.000599e-01

Kết quả dự đoán xác suất trung bình từ mô hình hồi quy logistic cho thấy mối quan hệ rất rõ giữa trình độ học vấn và khả năng đạt thu nhập trên 50K. Những người có học vấn cao hơn có xác suất dự đoán cao hơn đáng kể.

Cụ thể, những cá nhân có bằng Doctorate hoặc Prof-school có xác suất thu nhập >50K lần lượt là 74.7% và 74.9%, cao nhất trong các nhóm. Tiếp theo là nhóm có bằng Masters với xác suất khoảng 56.4%, và Bachelors với 42.1%. Hai nhóm có trình độ cao đẳng (Associate academic và vocational) có xác suất lần lượt là 25.4% và 26.3%.

Ngược lại, những nhóm có trình độ học vấn thấp hơn có xác suất rất thấp. Ví dụ, nhóm Preschool gần như bằng 0, trong khi các nhóm từ 1st–4th đến 9th dao động chỉ từ 4% đến 7%. Ngay cả nhóm tốt nghiệp phổ thông (HS-grad) cũng chỉ đạt mức xác suất 16.4% – thấp hơn đáng kể so với nhóm có trình độ đại học trở lên.

# BƯỚC 10: Vẽ biểu đồ xác suất trung bình theo trình độ học vấn
library(ggplot2)

ggplot(mean_probs, aes(x = education, y = predicted_prob)) +
  geom_col(fill = "skyblue", color = "black") +
  geom_text(aes(label = round(predicted_prob, 2)),
            vjust = -0.5, size = 3.5) +
  labs(title = "Xác suất trung bình thu nhập >50K theo trình độ học vấn",
       x = "Trình độ học vấn", y = "Xác suất dự đoán") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        plot.title = element_text(size = 16, hjust = 0.5, face = "bold"))

10. Tác động của giới tính, học vấn, nghề nghiệp và tình trạng hôn nhân đến xác suất một cá nhân có thu nhập cao

Trong phân tích này, em xây dựng mô hình hồi quy logistic để đánh giá tác động đồng thời của bốn đặc điểm nhân khẩu học đến xác suất có thu nhập cao hơn 50K. Cụ thể, các biến được đưa vào mô hình bao gồm:

  • education: trình độ học vấn
  • marital.status: tình trạng hôn nhân.
  • occupation: nhóm nghề nghiệp.
  • sex: giới tính.

Chọn nhóm tham chiếu cho từng biến để mô hình hóa và diễn giải kết quả chính xác hơn. Cụ thể: - Trình độ học vấn: nhóm tham chiếu là "Doctorate" (cao nhất) - Tình trạng hôn nhân: "Never-married" - Nghề nghiệp: "Exec-managerial" (có thể xem là nhóm nghề có vị thế cao) - Giới tính: "Male"

Việc đặt các nhóm này làm tham chiếu giúp ta hiểu rõ tác động tương đố* của các nhóm còn lại so với nhóm chuẩn trong từng biến.

Mô hình có dạng tổng quát:

\[ \log\left(\frac{P(Y = 1)}{1 - P(Y = 1)}\right) = \beta_0 + \sum_j \beta_j \cdot \text{Education}_j + \sum_k \beta_k \cdot \text{Marital}_k + \sum_l \beta_l \cdot \text{Occupation}_l + \beta_s \cdot \text{Sex} \]

Trong đó: - \(Y = 1\): thu nhập >50K - Mỗi nhóm con trong từng biến định tính được biểu diễn bởi một biến giả (dummy variable) - Mỗi hệ số \(\beta_j\) thể hiện sự thay đổi log-odds của khả năng có thu nhập cao khi chuyển từ nhóm tham chiếu sang nhóm đang xét.

# 2. Tạo biến nhị phân từ income (1 = >50K, 0 = <=50K)
df$income_bin <- ifelse(trimws(df$income) == ">50K", 1, 0)

# 3. Chuyển các biến định tính thành factor
df$education <- as.factor(df$education)
df$marital.status <- as.factor(df$marital.status)
df$occupation <- as.factor(df$occupation)
df$sex <- as.factor(df$sex)

# 4. Đặt lại nhóm gốc (reference level) cho từng biến
df$education <- relevel(df$education, ref = "Doctorate")
df$marital.status <- relevel(df$marital.status, ref = "Never-married")
df$occupation <- relevel(df$occupation, ref = "Exec-managerial")
df$sex <- relevel(df$sex, ref = "Male")

# 5. Xây dựng mô hình hồi quy logistic
model_multi_cat <- glm(
  income_bin ~ education + marital.status + occupation + sex,
  data = df, family = binomial(link = "logit")
)

# 6. Xem kết quả chi tiết
summary(model_multi_cat)
## 
## Call:
## glm(formula = income_bin ~ education + marital.status + occupation + 
##     sex, family = binomial(link = "logit"), data = df)
## 
## Coefficients:
##                                      Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                          -0.16796    0.15750  -1.066 0.286236    
## education10th                        -3.22193    0.20802 -15.489  < 2e-16 ***
## education11th                        -3.24894    0.20715 -15.684  < 2e-16 ***
## education12th                        -2.87042    0.25772 -11.138  < 2e-16 ***
## education1st-4th                     -3.89284    0.45450  -8.565  < 2e-16 ***
## education5th-6th                     -3.78180    0.33857 -11.170  < 2e-16 ***
## education7th-8th                     -3.72593    0.23514 -15.846  < 2e-16 ***
## education9th                         -3.59584    0.26157 -13.747  < 2e-16 ***
## educationAssoc-acdm                  -1.97185    0.17028 -11.580  < 2e-16 ***
## educationAssoc-voc                   -2.00638    0.16480 -12.175  < 2e-16 ***
## educationBachelors                   -1.32256    0.14987  -8.825  < 2e-16 ***
## educationHS-grad                     -2.48121    0.15303 -16.214  < 2e-16 ***
## educationMasters                     -0.83200    0.15649  -5.317 1.06e-07 ***
## educationPreschool                  -14.59264  112.33026  -0.130 0.896639    
## educationProf-school                 -0.06052    0.18531  -0.327 0.743995    
## educationSome-college                -2.17718    0.15283 -14.246  < 2e-16 ***
## marital.statusDivorced                0.89382    0.07299  12.245  < 2e-16 ***
## marital.statusMarried-AF-spouse       3.46772    0.48050   7.217 5.32e-13 ***
## marital.statusMarried-civ-spouse      2.82791    0.05592  50.569  < 2e-16 ***
## marital.statusMarried-spouse-absent   0.76568    0.20794   3.682 0.000231 ***
## marital.statusSeparated               0.67493    0.14444   4.673 2.97e-06 ***
## marital.statusWidowed                 1.15221    0.13662   8.434  < 2e-16 ***
## occupationAdm-clerical               -0.91650    0.06891 -13.299  < 2e-16 ***
## occupationArmed-Forces               -1.62413    1.25315  -1.296 0.194964    
## occupationCraft-repair               -0.95889    0.05992 -16.004  < 2e-16 ***
## occupationFarming-fishing            -1.76700    0.11474 -15.400  < 2e-16 ***
## occupationHandlers-cleaners          -1.84776    0.12798 -14.438  < 2e-16 ***
## occupationMachine-op-inspct          -1.33412    0.08576 -15.557  < 2e-16 ***
## occupationOther-service              -1.95708    0.10307 -18.987  < 2e-16 ***
## occupationPriv-house-serv            -3.10010    1.02181  -3.034 0.002414 ** 
## occupationProf-specialty             -0.44179    0.05934  -7.445 9.66e-14 ***
## occupationProtective-serv            -0.53288    0.10450  -5.099 3.41e-07 ***
## occupationSales                      -0.58012    0.05988  -9.689  < 2e-16 ***
## occupationTech-support               -0.37155    0.09485  -3.917 8.96e-05 ***
## occupationTransport-moving           -1.00029    0.08129 -12.305  < 2e-16 ***
## sexFemale                            -0.38820    0.04859  -7.989 1.36e-15 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 33851  on 30161  degrees of freedom
## Residual deviance: 22639  on 30126  degrees of freedom
## AIC: 22711
## 
## Number of Fisher Scoring iterations: 13
# 7. Chuyển hệ số thành odds ratio và khoảng tin cậy
odds_ratios <- exp(coef(model_multi_cat))  # Tỷ lệ odds
conf_intervals <- exp(confint.default(model_multi_cat))  # Khoảng tin cậy 95%

# 8. Tạo bảng kết quả dễ đọc
result_table <- round(cbind(
  Odds_Ratio = odds_ratios,
  CI_lower = conf_intervals[, 1],
  CI_upper = conf_intervals[, 2]
), 3)

# In bảng kết quả
print(result_table)
##                                     Odds_Ratio CI_lower     CI_upper
## (Intercept)                              0.845    0.621 1.151000e+00
## education10th                            0.040    0.027 6.000000e-02
## education11th                            0.039    0.026 5.800000e-02
## education12th                            0.057    0.034 9.400000e-02
## education1st-4th                         0.020    0.008 5.000000e-02
## education5th-6th                         0.023    0.012 4.400000e-02
## education7th-8th                         0.024    0.015 3.800000e-02
## education9th                             0.027    0.016 4.600000e-02
## educationAssoc-acdm                      0.139    0.100 1.940000e-01
## educationAssoc-voc                       0.134    0.097 1.860000e-01
## educationBachelors                       0.266    0.199 3.570000e-01
## educationHS-grad                         0.084    0.062 1.130000e-01
## educationMasters                         0.435    0.320 5.910000e-01
## educationPreschool                       0.000    0.000 1.897519e+89
## educationProf-school                     0.941    0.655 1.353000e+00
## educationSome-college                    0.113    0.084 1.530000e-01
## marital.statusDivorced                   2.444    2.119 2.820000e+00
## marital.statusMarried-AF-spouse         32.064   12.503 8.222700e+01
## marital.statusMarried-civ-spouse        16.910   15.155 1.886900e+01
## marital.statusMarried-spouse-absent      2.150    1.431 3.232000e+00
## marital.statusSeparated                  1.964    1.480 2.607000e+00
## marital.statusWidowed                    3.165    2.422 4.137000e+00
## occupationAdm-clerical                   0.400    0.349 4.580000e-01
## occupationArmed-Forces                   0.197    0.017 2.298000e+00
## occupationCraft-repair                   0.383    0.341 4.310000e-01
## occupationFarming-fishing                0.171    0.136 2.140000e-01
## occupationHandlers-cleaners              0.158    0.123 2.030000e-01
## occupationMachine-op-inspct              0.263    0.223 3.120000e-01
## occupationOther-service                  0.141    0.115 1.730000e-01
## occupationPriv-house-serv                0.045    0.006 3.340000e-01
## occupationProf-specialty                 0.643    0.572 7.220000e-01
## occupationProtective-serv                0.587    0.478 7.200000e-01
## occupationSales                          0.560    0.498 6.300000e-01
## occupationTech-support                   0.690    0.573 8.310000e-01
## occupationTransport-moving               0.368    0.314 4.310000e-01
## sexFemale                                0.678    0.617 7.460000e-01

Mô hình hồi quy logistic được xây dựng để phân tích các yếu tố ảnh hưởng đến khả năng một cá nhân có thu nhập cao (income_bin = 1). Các biến độc lập gồm: trình độ học vấn, tình trạng hôn nhân, nghề nghiệp, và giới tính.

Trình độ học vấn

Kết quả cho thấy học vấn có ảnh hưởng rất rõ rệt. So với nhóm có học vấn cao nhất, những người chỉ học đến phổ thông có xác suất thu nhập cao thấp hơn đáng kể. Ví dụ, người học hết lớp 10 chỉ có khoảng 4% cơ hội có thu nhập cao so với nhóm tham chiếu. Trong khi đó, người có bằng Cử nhân hay Thạc sĩ có xác suất cao hơn, nhưng vẫn chưa bằng nhóm học vấn cao nhất. Điều này cho thấy học vấn càng cao thì khả năng có thu nhập cao càng lớn.

Tình trạng hôn nhân

Những người đã kết hôn hợp pháp có khả năng thu nhập cao gấp gần 17 lần so với người chưa từng kết hôn. Người góa hoặc ly hôn cũng có xác suất cao hơn. Điều này có thể phản ánh sự ổn định hoặc độ tuổi cao hơn dẫn đến thu nhập tốt hơn.

Nghề nghiệp

Nghề nghiệp ảnh hưởng mạnh đến thu nhập. So với các ngành quản lý hoặc chuyên môn cao, các nghề như lao động phổ thông, nông nghiệp hay dịch vụ có odds thấp hơn rõ rệt. Nghề nghiệp càng chuyên môn, càng kỹ thuật thì xác suất thu nhập cao càng lớn.

Giới tính

Phụ nữ có xác suất thu nhập cao thấp hơn 32% so với nam giới, và sự khác biệt này có ý nghĩa thống kê. Điều này phản ánh phần nào khoảng cách thu nhập theo giới.

# 9. Dự đoán xác suất thu nhập >50K theo từng trình độ học vấn
# Giữ các yếu tố khác cố định: Male, Married-civ-spouse, Prof-specialty
edu_levels <- levels(df$education)

# Tạo tập dữ liệu giả định với các trình độ học vấn
new_data <- data.frame(
  education = edu_levels,
  marital.status = factor("Married-civ-spouse", levels = levels(df$marital.status)),
  occupation = factor("Prof-specialty", levels = levels(df$occupation)),
  sex = factor("Male", levels = levels(df$sex))
)

# Dự đoán xác suất
predicted_probs <- predict(model_multi_cat, newdata = new_data, type = "response")

# Tạo bảng kết quả
education_probs <- data.frame(
  Education = edu_levels,
  Probability = round(predicted_probs, 3)
)


# In kết quả
print(education_probs)
##       Education Probability
## 1     Doctorate       0.902
## 2          10th       0.268
## 3          11th       0.263
## 4          12th       0.342
## 5       1st-4th       0.158
## 6       5th-6th       0.173
## 7       7th-8th       0.181
## 8           9th       0.201
## 9    Assoc-acdm       0.561
## 10    Assoc-voc       0.553
## 11    Bachelors       0.710
## 12      HS-grad       0.435
## 13      Masters       0.800
## 14    Preschool       0.000
## 15  Prof-school       0.896
## 16 Some-college       0.510

Kết quả phân tích cho thấy trình độ học vấn có ảnh hưởng rõ rệt đến xác suất đạt thu nhập trên 50.000 USD, ngay cả khi đã kiểm soát các yếu tố như giới tính, nghề nghiệp và tình trạng hôn nhân. Những người có bằng Doctorate, Prof-school hoặc Masters có xác suất rất cao đạt thu nhập cao, lần lượt là 90,2%, 89,6% và 80%. Trong khi đó, nhóm có học vấn thấp như “1st-4th” hay “9th” chỉ có xác suất khoảng 15–20%, thậm chí nhóm “Preschool” gần như bằng 0%.

Những người học cao đẳng học thuật hoặc cao đẳng nghề cũng có xác suất khá cao (trên 55%), cao hơn cả nhóm chỉ học hết phổ thông (43,5%) hay “Some-college” (51%). Điều này cho thấy không chỉ bằng cấp cao mà cả các hình thức đào tạo chuyên môn cũng giúp tăng khả năng có thu nhập cao.

Tóm lại, trình độ học vấn càng cao thì xác suất đạt thu nhập trên 50K càng lớn, thể hiện vai trò quan trọng của giáo dục trong việc nâng cao thu nhập cá nhân

Phần 4: Mô hình hồi quy Probit

Hồi quy Probit là một phương pháp thay thế cho hồi quy logistic khi xử lý mô hình với biến phụ thuộc nhị phân Khác với mô hình logistic sử dụng hàm logistic (sigmoid) để ánh xạ tuyến tính sang xác suất, hồi quy Probit sử dụng hàm phân phối chuẩn tích lũy (CDF) chuẩn.

  • Trong khi mô hình logistic giả định lỗi tuân theo phân phối logistic, mô hình Probit giả định lỗi có phân phối chuẩn chuẩn hóa.
  • Cả hai mô hình đều phù hợp cho biến nhị phân, và thường cho kết quả khá giống nhau.
  • Tuy nhiên, mô hình Probit thường được sử dụng trong các nghiên cứu xã hội học, kinh tế học, và tài chính, khi người phân tích tin rằng hành vi ra quyết định mang tính liên tục tiềm ẩn và chịu ảnh hưởng từ các yếu tố ngẫu nhiên có phân phối chuẩn.

Cấu trúc mô hình

Giả sử ta có biến phụ thuộc nhị phân \(Y \in \{0, 1\}\) và một tập biến độc lập \(X\). Mô hình Probit được định nghĩa như sau:

\[ P(Y = 1 \mid X) = \Phi(\beta_0 + \beta_1 X_1 + \cdots + \beta_k X_k) \]

Trong đó:

  • \(\Phi(\cdot)\): là hàm phân phối chuẩn tích lũy (CDF).
  • \(\beta_0, \beta_1, \dots, \beta_k\): là các hệ số cần ước lượng.
  • Hàm \(\Phi\) đảm bảo đầu ra là một xác suất trong khoảng (0, 1).

Mô hình hóa bằng Probit có thể hiểu như việc giả định rằng đằng sau việc một cá nhân có thu nhập cao hay không là một biến liên tục tiềm ẩn, chỉ được quan sát dưới dạng nhị phân.

  • Các hệ số trong mô hình Probit không thể diễn giải trực tiếp như xác suất hay odds ratio.
  • Thay vào đó, chúng phản ánh ảnh hưởng của biến độc lập lên z-score của phân phối chuẩn.
  • Để có được mức độ thay đổi xác suất tương ứng với sự thay đổi của biến \(X_j\), cần tính hiệu ứng biên (marginal effect).

1. Ảnh hưởng của giới tính đến xác suất thu nhập cao bằng hồi quy Probit

Mô hình Probit có công thức như sau:

\[ P(Y = 1 \mid \text{sex}) = \Phi(\beta_0 + \beta_1 \cdot \text{sex}) \] Trong đó:

  • \(Y = 1\) nếu cá nhân có thu nhập >50K.
  • \(\Phi\)hàm phân phối tích lũy chuẩn chuẩn hóa (Standard Normal CDF).
  • Hệ số \(\beta_1\) cho biết mức độ thay đổi của z-score khi chuyển từ nhóm giới tính tham chiếu (mặc định là "Female") sang nhóm "Male".
# Thư viện cần thiết
library(ggplot2)
library(ggthemes)
library(scales)

# BƯỚC 1: Tạo biến nhị phân từ income (1 = ">50K", 0 = "<=50K")
df$income_bin <- ifelse(trimws(df$income) == ">50K", 1, 0)

# BƯỚC 2: Chuyển biến sex thành factor
df$sex <- as.factor(df$sex)

# BƯỚC 3: Chạy mô hình hồi quy Probit
model_probit <- glm(income_bin ~ sex, data = df, family = binomial(link = "probit"))

# BƯỚC 4: Xem kết quả mô hình
summary(model_probit)
## 
## Call:
## glm(formula = income_bin ~ sex, family = binomial(link = "probit"), 
##     data = df)
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -0.485003   0.009165  -52.92   <2e-16 ***
## sexFemale   -0.722194   0.019024  -37.96   <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: 33851  on 30161  degrees of freedom
## Residual deviance: 32287  on 30160  degrees of freedom
## AIC: 32291
## 
## Number of Fisher Scoring iterations: 4

Kết quả cho thấy hệ số ước lượng cho biến sexFemale là -0.722, với giá trị p < 2e-16, cho thấy mối quan hệ này có ý nghĩa thống kê ở mức rất cao.

Trong mô hình Probit, hệ số hồi quy phản ánh mức độ ảnh hưởng của biến giải thích lên xác suất chuẩn hóa (z-score) của biến phụ thuộc. Cụ thể:

Intercept = -0.485: Đây là z-score trung bình (chưa chuẩn hóa thành xác suất) cho nam giới (vì nam là nhóm tham chiếu).

sexFemale = -0.722: Nữ giới có điểm z thấp hơn nam giới trung bình 0.722 đơn vị, tức là xác suất có thu nhập >50K của họ thấp hơn đáng kể so với nam giới, giữ các yếu tố khác không đổi (ở đây là mô hình đơn biến nên không có yếu tố nào khác).

# Odds ratio không áp dụng trực tiếp cho probit (vì hệ số là theo chuẩn Z), 
# nhưng bạn vẫn có thể xem hệ số hoặc tính xác suất

# BƯỚC 5: Dự đoán xác suất thu nhập >50K
df$predicted_prob <- predict(model_probit, type = "response")

# BƯỚC 6: Tính xác suất trung bình theo giới tính
prob_by_sex <- aggregate(predicted_prob ~ sex, data = df, FUN = mean)
print(prob_by_sex)
##      sex predicted_prob
## 1   Male      0.3138371
## 2 Female      0.1136782

Dựa trên kết quả ước lượng từ mô hình hồi quy Probit, ta có thể thấy sự khác biệt rõ rệt về xác suất có thu nhập cao giữa hai giới.

Cụ thể, đối với nam giới, xác suất dự đoán đạt mức thu nhập trên 50.000 đô mỗi năm là khoảng 31.4%, trong khi với nữ giới, xác suất này chỉ đạt khoảng 11.4%.

# BƯỚC 7: Vẽ biểu đồ xác suất
ggplot(prob_by_sex, aes(x = sex, y = predicted_prob, fill = sex)) +
  geom_col(width = 0.6, color = "black", show.legend = FALSE) +
  geom_text(
    aes(label = percent(predicted_prob, accuracy = 0.1)),
    vjust = -0.5, size = 6, fontface = "bold"
  ) +
  scale_y_continuous(labels = percent_format(accuracy = 1), limits = c(0, 0.6)) +
  scale_fill_manual(values = c("Female" = "#FF9999", "Male" = "#6699CC")) +
  labs(
    title = "Xác suất trung bình có thu nhập >50K theo giới tính (Probit)",
    x = "Giới tính",
    y = "Xác suất dự đoán"
  ) +
  theme_minimal() +
  theme(
    plot.title = element_text(hjust = 0.5, size = 18, face = "bold"),
    axis.title = element_text(size = 14),
    axis.text = element_text(size = 13),
    panel.grid.major.y = element_line(color = "gray85")
  )

2. Đánh giá tác động của giới tính, học vấn, nghề nghiệp và tình trạng Hôn nhân đến xác suất một cá nhân có thu nhập cao (>50K) bằng mô hình hồi quy probit

Để xác định các yếu tố ảnh hưởng đến khả năng có thu nhập cao của một cá nhân, mô hình hồi quy Probit đã được ước lượng với biến phụ thuộc là xác suất thu nhập vượt ngưỡng 50.000 USD mỗi năm. Bốn biến độc lập được đưa vào mô hình dưới dạng nhị phân gồm: giới tính (nam hay nữ), trình độ học vấn (cao hay không), tình trạng hôn nhân (đã kết hôn hay chưa) và nghề nghiệp (quản lý điều hành hay không). Các nhóm tham chiếu lần lượt là: nữ giới, trình độ học vấn cao, đã kết hôn, và làm nghề quản lý điều hành.

# Biến nhị phân: Nam hay Nữ
df$sex_binary <- ifelse(df$sex == "Male", 1, 0)  # Nam = 1, Nữ = 0

# Biến nhị phân: Trình độ học vấn cao (Bachelors, Masters, Prof-school, Doctorate)
high_edu_levels <- c("Bachelors", "Masters", "Prof-school", "Doctorate")
df$edu_high <- ifelse(df$education %in% high_edu_levels, 1, 0)

# Biến nhị phân: Đã kết hôn (Married)
married_status <- c("Married-civ-spouse", "Married-spouse-absent", "Married-AF-spouse")
df$married <- ifelse(df$marital.status %in% married_status, 1, 0)

# Biến nhị phân: Làm nghề quản lý điều hành hay không
df$managerial <- ifelse(df$occupation == "Exec-managerial", 1, 0)
# Mô hình hồi quy Probit
model_probit_binary <- glm(
  income_bin ~ sex_binary + edu_high + married + managerial,
  family = binomial(link = "probit"),
  data = df
)

# Xem kết quả
summary(model_probit_binary)
## 
## Call:
## glm(formula = income_bin ~ sex_binary + edu_high + married + 
##     managerial, family = binomial(link = "probit"), data = df)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -2.01020    0.02316 -86.804   <2e-16 ***
## sex_binary   0.20461    0.02402   8.519   <2e-16 ***
## edu_high     0.94090    0.01992  47.234   <2e-16 ***
## married      1.32742    0.02181  60.874   <2e-16 ***
## managerial   0.52647    0.02491  21.138   <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: 33851  on 30161  degrees of freedom
## Residual deviance: 24390  on 30157  degrees of freedom
## AIC: 24400
## 
## Number of Fisher Scoring iterations: 5

Kết quả mô hình cho thấy tất cả các biến đều có ý nghĩa thống kê cao ở mức 1% (p < 0.001) và có ảnh hưởng rõ rệt đến xác suất thu nhập cao:

Giới tính: Việc là nam giới làm tăng đáng kể xác suất thu nhập cao so với nữ giới, với hệ số ước lượng là 0.205. Điều này phản ánh một phần chênh lệch thu nhập theo giới tính trên thị trường lao động.

Trình độ học vấn: Việc thuộc nhóm trình độ học vấn thấp hơn (so với nhóm có trình độ từ đại học trở lên) làm giảm mạnh xác suất đạt mức thu nhập cao, với hệ số 0.941 cho nhóm học vấn cao (có dấu dương vì biến “edu_high” được mã hóa 1 nếu học vấn cao). Điều này nhấn mạnh vai trò then chốt của học vấn trong phân hóa thu nhập.

Tình trạng hôn nhân: Người chưa kết hôn có xác suất thu nhập cao thấp hơn đáng kể so với người đã kết hôn. Hệ số ước lượng cho biến “married” là 1.327, cho thấy tình trạng hôn nhân có liên hệ chặt chẽ với điều kiện kinh tế.

Nghề nghiệp: Làm việc trong nhóm nghề không phải quản lý điều hành làm giảm khả năng đạt mức thu nhập cao, với hệ số ước lượng 0.526 cho nhóm quản lý. Nghĩa là các vị trí điều hành mang lại lợi thế rõ rệt về thu nhập so với các nghề khác.

Về tổng thể, mô hình hồi quy Probit cho thấy cấu trúc thu nhập trên thị trường lao động có liên quan chặt chẽ đến đặc điểm cá nhân, đặc biệt là trình độ học vấn, tình trạng hôn nhân và ngành nghề. Trong đó, học vấn và hôn nhân là hai yếu tố có hệ số tác động lớn nhất, cho thấy vai trò nổi bật trong việc gia tăng xác suất đạt được mức thu nhập cao

# Tạo dữ liệu giả định (nữ, đã kết hôn, nghề quản lý)
edu_data <- data.frame(
  sex_binary = c(0, 0),         # Nữ
  edu_high = c(0, 1),           # 0 = học vấn thấp, 1 = học vấn cao
  married = c(1, 1),            # Đã kết hôn
  managerial = c(1, 1)          # Nghề quản lý
)

# Đặt nhãn cho dễ đọc
edu_labels <- c("Thấp", "Cao")

# Dự đoán xác suất từ mô hình probit nhị phân
predicted_probs <- predict(model_probit_binary, newdata = edu_data, type = "response")

# Tạo bảng kết quả
result <- data.frame(
  HocVan = edu_labels,
  XacSuat = round(predicted_probs, 3)
)

print(result)
##   HocVan XacSuat
## 1   Thấp   0.438
## 2    Cao   0.784

Khi giữ cố định các yếu tố khác như giới tính (nữ), tình trạng hôn nhân (đã kết hôn) và nghề nghiệp (quản lý/điều hành), mô hình hồi quy Probit cho thấy xác suất dự đoán một cá nhân có thu nhập vượt mức 50.000 đô la mỗi năm phụ thuộc rõ rệt vào trình độ học vấn.

Cụ thể, đối với nhóm có trình độ học vấn thấp (tức là không có bằng đại học trở lên), xác suất để đạt mức thu nhập >50K chỉ vào khoảng 43,8%. Ngược lại, với nhóm trình độ học vấn cao (bao gồm các bằng cấp như Cử nhân, Thạc sĩ, Chuyên ngành nghề và Tiến sĩ), xác suất này tăng đáng kể, đạt mức 78,4%.

library(ggplot2)

# Tạo data frame từ kết quả đã có
prob_data <- data.frame(
  HocVan = c("Thấp", "Cao"),
  XacSuat = c(0.438, 0.784)
)

# Vẽ biểu đồ cột
ggplot(prob_data, aes(x = HocVan, y = XacSuat, fill = HocVan)) +
  geom_col(width = 0.5) +
  geom_text(aes(label = scales::percent(XacSuat, accuracy = 0.1)), 
            vjust = -0.5, size = 5) +
  scale_y_continuous(labels = scales::percent_format(accuracy = 1)) +
  labs(
    title = "Xác suất dự đoán thu nhập >50K theo trình độ học vấn",
    x = "Trình độ học vấn",
    y = "Xác suất thu nhập >50K"
  ) +
  theme_minimal() +
  theme(
    text = element_text(size = 14),
    legend.position = "none"
  )

3. Ảnh hưởng của trình độ học vấn đến xác suất đạt được thu nhập cao

Nhằm đánh giá riêng biệt vai trò của trình độ học vấn đối với khả năng đạt mức thu nhập trên 50.000 USD mỗi năm, một mô hình hồi quy Probit được xây dựng với duy nhất một biến độc lập là education. Biến phụ thuộc là income_bin, được mã hóa nhị phân với giá trị 1 nếu cá nhân có thu nhập vượt ngưỡng 50K và 0 nếu không. Trong mô hình, nhóm trình độ “Doctorate” được chọn làm mốc tham chiếu để thuận tiện cho việc so sánh.

# Biến nhị phân từ income
df$income_bin <- ifelse(trimws(df$income) == ">50K", 1, 0)

# Chuyển education thành factor (nếu chưa)
df$education <- as.factor(df$education)

# Đặt nhóm tham chiếu là Doctorate
df$education <- relevel(df$education, ref = "Doctorate")
model_edu_only <- glm(
  income_bin ~ education,
  data = df,
  family = binomial(link = "probit")
)

# Xem kết quả
summary(model_edu_only)
## 
## Call:
## glm(formula = income_bin ~ education, family = binomial(link = "probit"), 
##     data = df)
## 
## Coefficients:
##                        Estimate Std. Error z value Pr(>|z|)    
## (Intercept)            0.664037   0.070183   9.461  < 2e-16 ***
## education10th         -2.125449   0.096208 -22.092  < 2e-16 ***
## education11th         -2.250672   0.094203 -23.892  < 2e-16 ***
## education12th         -2.090114   0.118193 -17.684  < 2e-16 ***
## education1st-4th      -2.417805   0.198301 -12.193  < 2e-16 ***
## education5th-6th      -2.395701   0.149666 -16.007  < 2e-16 ***
## education7th-8th      -2.195426   0.108892 -20.162  < 2e-16 ***
## education9th          -2.262724   0.119007 -19.013  < 2e-16 ***
## educationAssoc-acdm   -1.326091   0.082197 -16.133  < 2e-16 ***
## educationAssoc-voc    -1.297554   0.079488 -16.324  < 2e-16 ***
## educationBachelors    -0.862118   0.072399 -11.908  < 2e-16 ***
## educationHS-grad      -1.640856   0.071787 -22.857  < 2e-16 ***
## educationMasters      -0.502338   0.076814  -6.540 6.17e-11 ***
## educationPreschool    -5.483231  21.906733  -0.250    0.802    
## educationProf-school   0.007553   0.091359   0.083    0.934    
## educationSome-college -1.505444   0.072328 -20.814  < 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: 33851  on 30161  degrees of freedom
## Residual deviance: 29946  on 30146  degrees of freedom
## AIC: 29978
## 
## Number of Fisher Scoring iterations: 12

Kết quả hồi quy cho thấy hầu hết các hệ số ước lượng của các bậc học thấp hơn đều mang giá trị âm và có ý nghĩa thống kê cao (p < 0.001), điều này cho thấy rằng khi so sánh với nhóm có học vị tiến sĩ, những người có trình độ học vấn thấp hơn có xác suất đạt được mức thu nhập cao thấp hơn một cách đáng kể. Cụ thể, người có bằng thạc sĩ (Masters) có hệ số -0.50, trong khi người chỉ học hết cấp ba (HS-grad) có hệ số lên đến -1.64, cho thấy mức độ bất lợi tăng dần theo chiều giảm của học vấn.

# Tạo dữ liệu giả định
edu_levels <- levels(df$education)

new_data_edu <- data.frame(education = factor(edu_levels, levels = levels(df$education)))

# Dự đoán xác suất
predicted_probs_edu <- predict(model_edu_only, newdata = new_data_edu, type = "response")

# Gộp kết quả vào bảng
edu_prob_table <- data.frame(
  Education = edu_levels,
  Probability = round(predicted_probs_edu, 3)
)

print(edu_prob_table)
##       Education Probability
## 1     Doctorate       0.747
## 2          10th       0.072
## 3          11th       0.056
## 4          12th       0.077
## 5       1st-4th       0.040
## 6       5th-6th       0.042
## 7       7th-8th       0.063
## 8           9th       0.055
## 9    Assoc-acdm       0.254
## 10    Assoc-voc       0.263
## 11    Bachelors       0.421
## 12      HS-grad       0.164
## 13      Masters       0.564
## 14    Preschool       0.000
## 15  Prof-school       0.749
## 16 Some-college       0.200

Kết quả cho thấy, xác suất đạt được thu nhập cao có sự khác biệt rất rõ rệt giữa các nhóm trình độ. Những cá nhân có trình độ cao như “Prof-school” và “Doctorate” đạt xác suất lần lượt là 74.9% và 74.7%, trong khi nhóm “Masters” đạt 56.4% và “Bachelors” là 42.1%. Trong khi đó, những người chỉ tốt nghiệp trung học phổ thông (HS-grad) có xác suất 16.4%, còn những người học hết lớp 10 hoặc 11 chỉ còn khoảng 7% và 5.6%. Đặc biệt, những người học dưới cấp tiểu học như “1st–4th” hay “Preschool” có xác suất gần như bằng 0.

Điều này phản ánh mối quan hệ rất mạnh mẽ và gần như đơn điệu giữa trình độ học vấn và khả năng có thu nhập cao. Trình độ càng cao, xác suất vượt ngưỡng thu nhập càng lớn.

library(ggplot2)

ggplot(edu_prob_table, aes(x = reorder(Education, -Probability), y = Probability)) +
  geom_col(fill = "#4682B4") +
  coord_flip() +
  labs(
    title = "Xác suất đạt thu nhập >50K theo trình độ học vấn",
    x = "Trình độ học vấn",
    y = "Xác suất dự đoán"
  ) +
  theme_minimal(base_size = 13)

4. Đánh giá tác động của giới tính, học vấn, nghề nghiệp và tình trạng Hôn nhân đến xác suất một cá nhân có thu nhập cao (>50K) bằng mô hình hồi quy probit

# BƯỚC 1: Tạo biến nhị phân từ income
df$income_bin <- ifelse(trimws(df$income) == ">50K", 1, 0)

# BƯỚC 2: Chuyển các biến định tính thành factor
df$education <- as.factor(df$education)
df$marital.status <- as.factor(df$marital.status)
df$occupation <- as.factor(df$occupation)
df$sex <- as.factor(df$sex)

# Đặt lại nhóm gốc nếu muốn
df$education <- relevel(df$education, ref = "Doctorate")
df$marital.status <- relevel(df$marital.status, ref = "Never-married")
df$occupation <- relevel(df$occupation, ref = "Exec-managerial")
df$sex <- relevel(df$sex, ref = "Female")

# BƯỚC 3: Chạy mô hình hồi quy Probit
model_probit_multi <- glm(
  income_bin ~ education + marital.status + occupation + sex,
  data = df,
  family = binomial(link = "probit")
)

# BƯỚC 4: Xem kết quả mô hình
summary(model_probit_multi)
## 
## Call:
## glm(formula = income_bin ~ education + marital.status + occupation + 
##     sex, family = binomial(link = "probit"), data = df)
## 
## Coefficients:
##                                     Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                         -0.32905    0.08848  -3.719 0.000200 ***
## education10th                       -1.80341    0.11344 -15.897  < 2e-16 ***
## education11th                       -1.76196    0.11108 -15.862  < 2e-16 ***
## education12th                       -1.61200    0.14181 -11.367  < 2e-16 ***
## education1st-4th                    -2.19260    0.23094  -9.494  < 2e-16 ***
## education5th-6th                    -2.08745    0.17037 -12.252  < 2e-16 ***
## education7th-8th                    -2.05652    0.12367 -16.629  < 2e-16 ***
## education9th                        -1.95844    0.13573 -14.429  < 2e-16 ***
## educationAssoc-acdm                 -1.09216    0.09497 -11.500  < 2e-16 ***
## educationAssoc-voc                  -1.11818    0.09185 -12.174  < 2e-16 ***
## educationBachelors                  -0.71915    0.08270  -8.696  < 2e-16 ***
## educationHS-grad                    -1.39433    0.08452 -16.497  < 2e-16 ***
## educationMasters                    -0.43813    0.08659  -5.060 4.20e-07 ***
## educationPreschool                  -5.28737   18.26413  -0.289 0.772203    
## educationProf-school                -0.01369    0.10213  -0.134 0.893341    
## educationSome-college               -1.21991    0.08444 -14.448  < 2e-16 ***
## marital.statusDivorced               0.47678    0.03799  12.549  < 2e-16 ***
## marital.statusMarried-AF-spouse      1.97091    0.28851   6.831 8.41e-12 ***
## marital.statusMarried-civ-spouse     1.57595    0.02874  54.841  < 2e-16 ***
## marital.statusMarried-spouse-absent  0.42408    0.10931   3.880 0.000105 ***
## marital.statusSeparated              0.35497    0.07486   4.742 2.12e-06 ***
## marital.statusWidowed                0.64941    0.07186   9.037  < 2e-16 ***
## occupationAdm-clerical              -0.54023    0.03895 -13.871  < 2e-16 ***
## occupationArmed-Forces              -1.00115    0.70974  -1.411 0.158367    
## occupationCraft-repair              -0.56385    0.03508 -16.074  < 2e-16 ***
## occupationFarming-fishing           -1.00800    0.06371 -15.822  < 2e-16 ***
## occupationHandlers-cleaners         -1.05651    0.06893 -15.327  < 2e-16 ***
## occupationMachine-op-inspct         -0.79152    0.04885 -16.202  < 2e-16 ***
## occupationOther-service             -1.06968    0.05390 -19.847  < 2e-16 ***
## occupationPriv-house-serv           -1.52480    0.43718  -3.488 0.000487 ***
## occupationProf-specialty            -0.25615    0.03400  -7.534 4.93e-14 ***
## occupationProtective-serv           -0.29889    0.06161  -4.852 1.22e-06 ***
## occupationSales                     -0.34234    0.03469  -9.869  < 2e-16 ***
## occupationTech-support              -0.23147    0.05476  -4.227 2.37e-05 ***
## occupationTransport-moving          -0.58768    0.04728 -12.429  < 2e-16 ***
## sexMale                              0.24590    0.02700   9.108  < 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: 33851  on 30161  degrees of freedom
## Residual deviance: 22614  on 30126  degrees of freedom
## AIC: 22686
## 
## Number of Fisher Scoring iterations: 12

Kết quả hồi quy Probit cho thấy rằng các yếu tố trình độ học vấn, tình trạng hôn nhân, nghề nghiệp và giới tính đều có ảnh hưởng đáng kể đến xác suất một cá nhân đạt thu nhập trên 50.000 đô la mỗi năm.

Trình độ học vấn có tác động rõ rệt và nhất quán. So với nhóm có bằng tiến sĩ (Doctorate – nhóm tham chiếu), hầu hết các nhóm học vấn khác đều có hệ số âm và có ý nghĩa thống kê mạnh (p < 0.001), cho thấy xác suất thu nhập cao giảm mạnh khi trình độ học vấn thấp hơn. Ví dụ, những người chỉ học đến cấp 1 (1st-4th) có hệ số -2.19, còn người có bằng cử nhân là -0.72, phản ánh xác suất thấp hơn tương ứng khi so với nhóm tiến sĩ.

Tình trạng hôn nhân cũng cho thấy ảnh hưởng tích cực đến thu nhập. Người đã kết hôn (đặc biệt là nhóm “Married-civ-spouse”) có hệ số cao nhất (1.58), ngụ ý rằng xác suất thu nhập cao của họ vượt trội so với nhóm “Never-married” (nhóm tham chiếu). Các trạng thái khác như ly hôn, ly thân hoặc goá cũng có hệ số dương, tuy nhiên ở mức độ thấp hơn.

Nghề nghiệp là một yếu tố phân hoá mạnh. Các nhóm nghề phổ thông như lao động dịch vụ, nông nghiệp, thủ công, vận chuyển, làm việc nhà,… đều có hệ số âm đáng kể và có ý nghĩa thống kê, cho thấy xác suất đạt thu nhập cao ở các nhóm này thấp hơn đáng kể so với nhóm quản lý điều hành (Exec-managerial – nhóm tham chiếu). Ngược lại, các ngành kỹ thuật, hỗ trợ kỹ thuật hay chuyên môn có hệ số ít âm hơn, thể hiện sự phân tầng trong cơ hội thu nhập theo loại công việc.

Cuối cùng, giới tính tiếp tục cho thấy sự chênh lệch. Hệ số dương của biến sexMale (0.246) với mức ý nghĩa cao (p < 0.001) chỉ ra rằng nam giới có xác suất đạt thu nhập cao hơn nữ giới, giữ các yếu tố khác không đổi.

# Tạo dữ liệu đầu vào với từng trình độ học vấn
predict_probs <- function(edu) {
  data.frame(
    education = factor(edu, levels = levels(df$education)),
    marital.status = factor("Never-married", levels = levels(df$marital.status)),
    occupation = factor("Exec-managerial", levels = levels(df$occupation)),
    sex = factor("Male", levels = levels(df$sex))
  )
}

# Tạo tập dữ liệu giả định
education_levels <- levels(df$education)
new_data_list <- lapply(education_levels, predict_probs)
new_data <- do.call(rbind, new_data_list)

# 🔍 Dự đoán xác suất (tách riêng phần này)
predicted_probs <- predict(model_probit_multi, newdata = new_data, type = "response")

# Gán tên cho vector xác suất theo trình độ học vấn
names(predicted_probs) <- education_levels

# In ra với tên rõ ràng
print(round(predicted_probs, 3))
##    Doctorate         10th         11th         12th      1st-4th      5th-6th 
##        0.467        0.030        0.033        0.045        0.011        0.015 
##      7th-8th          9th   Assoc-acdm    Assoc-voc    Bachelors      HS-grad 
##        0.016        0.021        0.120        0.115        0.211        0.070 
##      Masters    Preschool  Prof-school Some-college 
##        0.301        0.000        0.461        0.096

Kết quả từ mô hình hồi quy Probit cho thấy rằng trình độ học vấn có ảnh hưởng rõ rệt đến xác suất đạt được mức thu nhập trên 50.000 USD, khi các yếu tố khác như giới tính, tình trạng hôn nhân và nghề nghiệp được giữ cố định. Cụ thể, những cá nhân có bằng tiến sĩ (Doctorate) hoặc chuyên ngành chuyên sâu (Prof-school) có xác suất cao nhất, lần lượt là 46.7% và 46.1%. Mức độ này tiếp tục duy trì ở mức cao đối với những người có bằng thạc sĩ (30.1%) và cử nhân (21.1%).

Ngược lại, những người chỉ hoàn thành trình độ trung học phổ thông (HS-grad) hoặc thấp hơn có xác suất rất thấp, thường dưới 10%. Đặc biệt, các nhóm có trình độ học vấn tiểu học hoặc chưa hoàn thành phổ thông như “1st–4th” hay “Preschool” gần như không có khả năng đạt được mức thu nhập cao, với xác suất gần bằng 0. Trong khi đó, những người học xong các chương trình cao đẳng (associate) có xác suất dao động từ 11.5% đến 12.0%, phản ánh mức độ trung gian giữa học vấn thấp và bậc đại học.

Tóm lại, xác suất có thu nhập cao tăng dần theo trình độ học vấn, cho thấy vai trò thiết yếu của giáo dục, đặc biệt ở các cấp bậc sau đại học, trong việc nâng cao khả năng đạt được mức thu nhập vượt ngưỡng 50.000 USD.

# Cài đặt ggplot2 nếu chưa có
# install.packages("ggplot2")

library(ggplot2)

# Tạo bảng dữ liệu để vẽ
plot_data <- data.frame(
  education = education_levels,
  probability = as.numeric(predicted_probs)
)

# Vẽ biểu đồ
ggplot(plot_data, aes(x = reorder(education, probability), y = probability)) +
  geom_col(fill = "steelblue") +
  coord_flip() +  # Xoay ngang để dễ đọc tên
  labs(
    title = "Xác suất dự đoán thu nhập >50K theo trình độ học vấn",
    x = "Trình độ học vấn",
    y = "Xác suất (giữ các yếu tố khác cố định)"
  ) +
  theme_minimal(base_size = 13)

#Đánh giá ảnh hưởng cận biên của từng biến
library(margins)
marginal_effects <- margins(model_probit_multi)
summary(marginal_effects)

Dựa trên kết quả phân tích hiệu ứng biên (Average Marginal Effects – AME) từ mô hình hồi quy Probit, ta có thể rút ra một số nhận định quan trọng về ảnh hưởng của nghề nghiệp và giới tính đến xác suất có thu nhập trên 50.000 USD, trong điều kiện các yếu tố khác được giữ cố định.

Cụ thể, nam giới có xác suất đạt thu nhập cao cao hơn nữ giới, với AME ước lượng là +5.1%, và giá trị z lớn (z = 9.23), cho thấy ảnh hưởng này có ý nghĩa thống kê cao. Đây là bằng chứng rõ ràng cho thấy tồn tại khoảng cách thu nhập giữa các giới trong mẫu khảo sát.

Ở chiều ngược lại, tất cả các nghề nghiệp được liệt kê đều có ảnh hưởng âm đến xác suất thu nhập cao so với nhóm tham chiếu (có thể là “Exec-managerial” – quản lý điều hành). Ví dụ, nhóm “transport-moving” (vận chuyển) có mức giảm xác suất mạnh nhất, lên tới -14.1%, tiếp theo là “sales” (bán hàng) với -8.6% và “protective-serv” (dịch vụ bảo vệ) với -7.5%. Ngay cả các ngành như “tech-support” (hỗ trợ kỹ thuật), vốn yêu cầu kỹ năng chuyên môn nhất định, cũng có ảnh hưởng âm (-5.9%).

Tóm lại, kết quả này cho thấy nghề nghiệp có ảnh hưởng đáng kể đến khả năng đạt được thu nhập cao, và sự khác biệt theo giới tính vẫn còn hiện hữu một cách rõ rệt trong thị trường lao động.