| Tiểu luận môn học | Phân tích dữ liệu định tính |
|---|---|
| Chuyên ngành | Tài chính định lượng |
| Giảng viên hướng dẫn | ThS. Trần Mạnh Tường |
| Sinh viên thực hiện | Nguyễn Hữu Hào |
| MSSV | 2221000300 |
Em xin cam đoan đề tài tiểu luận môn phân tích dữ liệu định tính “NGHIÊN CỨU CÁC YẾU TỐ ẢNH HƯỞNG ĐẾN MỨC LƯƠNG CỦA KỸ SƯ” là công trình nghiên cứu của em dưới sự hướng dẫn của giảng viên Th.S Trần Mạnh Tường.
Những phần sử dụng tài liệu tham khảo trong đề tài đều đã được trình bày nguồn đầy đủ và dựa trên bài giảng để thực hiện phân tích. Các thông tin và số liệu sử dụng trong bài tiểu luận là trung thực, có nguồn gốc rõ ràng. Các kết quả trình bày trong bài hoàn toàn là kết quả do em tự nghiên cứu và thực hiện.
Nếu vi phạm lời cam đoan này, em xin chịu hoàn toàn trách nhiệm.
Trước tiên, em xin chân thành cảm ơn Thầy Trần Mạnh Tường đã tận tình hướng dẫn và góp ý cho em trong suốt quá trình học tập và thực hiện tiểu luận này. Em cũng xin gửi lời cảm ơn đến các Thầy Cô trong khoa Khoa học dữ liệu, Trường Đại học Tài chính – Marketing đã nhiệt tình giảng dạy, truyền đạt các kiến thức về kinh tế – tài chính để em có thể vận dụng trong bài tiểu luận.
Môn Phân tích dữ liệu định tính là một môn học rất thú vị, bổ ích và mang tính thực tiễn cao, giúp em có cái nhìn rộng hơn về thống kê số liệu kinh tế và xã hội trong và ngoài nước. Tuy nhiên, do kiến thức của em còn hạn chế và khả năng tiếp thu thực tế còn gặp nhiều khó khăn, mặc dù đã cố gắng hết sức nhưng bài tiểu luận khó tránh khỏi thiếu sót. Em kính mong quý Thầy Cô xem xét và góp ý để bài tiểu luận của em được hoàn thiện hơn.
Sau cùng, em xin gửi lời cảm ơn đến các bạn trong lớp Toán Kinh Tế - 22DTL01 đã hỗ trợ, chia sẻ và động viên để em hoàn thành bài tiểu luận này.
Trước khi đi vào phân tích, em xin giới thiệu sơ lược về bộ dữ liệu “Data Engineer Salary” được thu thập từ trang Kaggle. Việc trình bày tổng quan này nhằm giúp người đọc có cái nhìn rõ ràng về cấu trúc dữ liệu, từ đó việc đánh giá và diễn giải kết quả phân tích trở nên trực quan hơn.
Bộ dữ liệu “Data Engineer Salary” phản ánh thông tin chi tiết về mức lương và các yếu tố liên quan đến công việc của kỹ sư dữ liệu. Dữ liệu gồm 16.534 quan sát với 11 biến, bao gồm các trường quan trọng như: mức lương, chức danh công việc, cấp độ kinh nghiệm, loại hình việc làm, nơi cư trú của kỹ sư, tỷ lệ làm việc từ xa, địa điểm công ty và quy mô công ty.
Những thông tin này cung cấp nền tảng quan trọng để phân tích xu hướng lương, so sánh các mô hình việc làm cũng như khám phá sự khác biệt về mức thu nhập theo khu vực và loại công việc trong ngành kỹ sư dữ liệu.
Mô tả các đặc trưng của dữ liệu
Experience_level: Biến thể hiện cấp độ kinh nghiệm của kỹ sư, được phân thành 4 nhóm:
Employment_type: Hình thức làm việc của kỹ sư, bao gồm:
Job_title: Chức danh hoặc vị trí công việc của kỹ sư trong tổ chức.
Salary: Mức lương của kỹ sư tính theo đơn vị tiền tệ địa phương.
Salary_currency: Loại tiền tệ được sử dụng để trả lương (ví dụ: USD).
Salary_in_usd: Mức lương đã được quy đổi sang đô la Mỹ để chuẩn hóa giữa các quốc gia.
Employee_residence: Quốc gia hoặc khu vực nơi kỹ sư sinh sống.
Remote_ratio: Tỷ lệ làm việc từ xa:
Company_location: Quốc gia hoặc khu vực đặt trụ sở của công ty nơi kỹ sư làm việc.
Company_size: Quy mô công ty dựa trên số lượng nhân viên:
library(readxl)
library(DT)
library(tidyverse)
library(ggplot2)
library(DescTools)
library(epitools)
library(kableExtra)
library(epiR)
data <- read_csv("D:/dldt/data.csv")
data<- data %>% select(-1)
d <- data
datatable(d)
Trong bối cảnh phát triển mạnh mẽ của nền kinh tế và công nghiệp hiện nay, vai trò của các kỹ sư ngày càng trở nên then chốt. Họ không chỉ tham gia vào quá trình thiết kế, xây dựng và vận hành các hệ thống kỹ thuật mà còn đóng góp vào việc đổi mới công nghệ và nâng cao năng suất. Tuy nhiên, mức lương của kỹ sư chịu ảnh hưởng bởi nhiều yếu tố như hình thức làm việc, kinh nghiệm, trình độ chuyên môn và quy mô doanh nghiệp.
Việc phân tích mức lương của các kỹ sư, đặc biệt là nhóm kỹ sư dữ liệu, không chỉ giúp người lao động nắm bắt xu hướng thị trường mà còn cung cấp cơ sở để họ định vị bản thân, đưa ra quyết định nghề nghiệp và đàm phán mức đãi ngộ phù hợp. Ở góc độ doanh nghiệp, nghiên cứu này giúp nhà quản lý xây dựng chính sách lương thưởng cạnh tranh và công bằng, từ đó thu hút và giữ chân nhân tài, nâng cao hiệu suất làm việc và chất lượng sản phẩm.
Ngoài ra, việc nghiên cứu mức lương còn góp phần thúc đẩy sự minh bạch trên thị trường lao động khi thông tin về chế độ đãi ngộ được công khai rõ ràng. Điều này tạo điều kiện để người lao động so sánh, đánh giá cơ hội việc làm và đồng thời giúp doanh nghiệp điều chỉnh chính sách nhân sự phù hợp với xu hướng chung của thị trường.
Mục tiêu của nghiên cứu này là xác định và phân tích các yếu tố chính ảnh hưởng đến mức lương của các kỹ sư. Cụ thể, nghiên cứu tập trung vào một số khía cạnh sau:
Nghiên cứu sử dụng phương pháp thống kê mô tả để xem tổng quan về bộ dữ liệu nghiên cứu. Tiếp đến xem xét các chỉ số odd ratio, relative risk, bảng tần số để đánh giá mối quan hệ giữa các cặp biến số sẽ ảnh hưởng đến mức lương. Ngoài ra, em còn thực hiện các mô hình xác suất tuyến tính, logit và probit để ước tính các xác suất để nhận được mức lương tương ứng.
Kết quả nghiên cứu mang lại cái nhìn rõ ràng về các yếu tố ảnh hưởng đến mức lương của kỹ sư, giúp họ đánh giá và định vị bản thân trên thị trường lao động. Từ đó, các kỹ sư có thể điều chỉnh định hướng nghề nghiệp, cải thiện trình độ chuyên môn và nâng cao năng lực cạnh tranh để đạt được mức lương và chính sách đãi ngộ tốt hơn. Đồng thời, nghiên cứu cũng cung cấp cơ sở để kỹ sư lựa chọn loại hình công việc, nhóm công việc và hình thức làm việc phù hợp với mục tiêu cá nhân.
Đối với các doanh nghiệp, kết quả nghiên cứu là nguồn thông tin quan trọng để xây dựng chính sách lương thưởng hợp lý, công bằng và cạnh tranh hơn. Việc hiểu rõ xu hướng thị trường và các yếu tố tác động đến mức lương giúp doanh nghiệp điều chỉnh chiến lược đãi ngộ, thu hút và giữ chân nhân tài, đồng thời tối ưu chi phí nhân sự. Qua đó, doanh nghiệp có thể nâng cao hiệu suất làm việc, cải thiện chất lượng sản phẩm và quản trị nguồn nhân lực một cách hiệu quả hơn.
Mô hình hồi quy xác suất tuyến tính (Linear Probability Model – LPM) là một dạng hồi quy tuyến tính được sử dụng trong các bài toán mà biến phụ thuộc chỉ nhận hai giá trị 0 hoặc 1. Trong nghiên cứu marketing, mô hình này đặc biệt hữu ích để ước lượng xác suất xảy ra một sự kiện, chẳng hạn như khả năng khách hàng đồng ý đăng ký dịch vụ sau một chiến dịch tiếp thị.
Mô hình giả định rằng xác suất sự kiện xảy ra (Y = 1) có thể được biểu diễn dưới dạng tuyến tính theo các biến độc lập. Công thức tổng quát của LPM được viết như sau:
\[ P(Y = 1 | X) = \beta_0 + \beta_1 X_1 + \beta_2 X_2 + \cdots + \beta_k X_k + u \]
Trong đó \(P(Y = 1 | X)\) là xác suất sự kiện xảy ra, \(X_i\) là các biến độc lập, \(\beta_i\) là hệ số hồi quy thể hiện tác động của biến \(X_i\) đến xác suất xảy ra sự kiện, và \(u\) là sai số ngẫu nhiên. Nhờ cấu trúc đơn giản, mô hình LPM có thể được ước lượng bằng phương pháp bình phương tối thiểu (OLS), đồng thời hệ số hồi quy được diễn giải trực tiếp như mức thay đổi trong xác suất khi biến độc lập tăng thêm một đơn vị.
Tuy nhiên, LPM cũng có một số hạn chế. Do giả định tuyến tính, mô hình có thể đưa ra xác suất dự đoán vượt ngoài khoảng [0;1], điều này không hợp lý về mặt xác suất. Bên cạnh đó, mô hình thường gặp vấn đề phương sai sai số thay đổi (heteroskedasticity), khiến sai số chuẩn bị sai lệch và làm giảm độ tin cậy của các kiểm định thống kê. Vì vậy, trong thực tiễn, LPM thường chỉ được sử dụng như một bước phân tích mô tả ban đầu trước khi áp dụng các mô hình hồi quy phi tuyến phù hợp hơn.
Mô hình hồi quy logit (Logistic Regression) thường được sử dụng trong các nghiên cứu khi biến phụ thuộc là biến nhị phân. Khác với LPM, hồi quy logit không giả định xác suất theo đường tuyến tính mà dùng hàm logistic để đưa toàn bộ giá trị dự đoán vào trong khoảng \([0;1]\), đảm bảo phù hợp với bản chất của xác suất. Nhờ đó, mô hình này cho phép phân tích và dự đoán khả năng xảy ra một hành vi cụ thể, ví dụ như khách hàng “đồng ý” hay “từ chối” đăng ký dịch vụ ngân hàng.
Công thức của mô hình logit có thể biểu diễn theo hai cách. Ở dạng xác suất:
\[ P(Y = 1 | X) = \frac{e^{\beta_0 + \beta_1 X_1 + \cdots + \beta_k X_k}}{1 + e^{\beta_0 + \beta_1 X_1 + \cdots + \beta_k X_k}} \]
Ở dạng log-odds (logit):
\[ \ln\left(\frac{P(Y = 1 | X)}{1 - P(Y = 1 | X)}\right) = \beta_0 + \beta_1 X_1 + \cdots + \beta_k X_k \]
Trong đó \(\frac{P}{1-P}\) là tỷ số chênh (odds), phản ánh tỷ lệ giữa xác suất sự kiện xảy ra và không xảy ra. Các hệ số \(\beta_i\) thể hiện tác động của biến độc lập đến log-odds; khi lấy số mũ của \(\beta_i\), ta thu được odds ratio (OR), cho biết mức thay đổi của odds khi biến \(X_i\) tăng thêm một đơn vị.
Hồi quy logit có nhiều ưu điểm. Mô hình đảm bảo xác suất dự đoán luôn nằm trong \([0;1]\) và cung cấp cách diễn giải thông qua odds và odds ratio, rất hữu ích trong các lĩnh vực như marketing, y tế hay khoa học xã hội. Tuy nhiên, hạn chế của logit là các hệ số \(\beta_i\) khó diễn giải trực tiếp như trong LPM, đòi hỏi kiến thức sâu hơn để hiểu kết quả. Việc ước lượng cũng phức tạp hơn, dựa trên phương pháp cực đại hóa khả năng xảy ra (Maximum Likelihood Estimation – MLE) và cần một kích thước mẫu đủ lớn để đảm bảo độ tin cậy.
Mô hình hồi quy probit là một dạng hồi quy nhị phân, tương tự như mô hình logit, được sử dụng khi biến phụ thuộc chỉ nhận hai giá trị (0/1). Điểm khác biệt chính của probit nằm ở hàm liên kết: thay vì dùng hàm logistic như logit, probit sử dụng hàm phân phối tích lũy chuẩn (Standard Normal Cumulative Distribution Function – CDF) để ước lượng xác suất. Điều này có nghĩa là xác suất xảy ra sự kiện được mô hình hóa dựa trên phân phối chuẩn.
Công thức tổng quát của mô hình probit:
\[ P(Y = 1 | X) = \Phi(\beta_0 + \beta_1 X_1 + \beta_2 X_2 + \cdots + \beta_k X_k) \]
Trong đó \(\Phi(\cdot)\) là hàm phân phối tích lũy chuẩn. Cách áp dụng probit tương tự như logit: mô hình được dùng để dự đoán xác suất xảy ra sự kiện hoặc đánh giá tác động của các biến độc lập đến khả năng xảy ra sự kiện. Việc ước lượng các hệ số \(\beta_i\) dựa trên phương pháp tối đa hóa hàm hợp lý (Maximum Likelihood Estimation – MLE).
Ưu điểm của probit là giả định xác suất theo phân phối chuẩn, phù hợp với nhiều hiện tượng trong thực tế, và tránh được vấn đề dự đoán xác suất vượt ngoài \([0;1]\) như ở LPM. Nhược điểm là các hệ số khó diễn giải trực tiếp; để hiểu rõ tác động của biến độc lập, nhà nghiên cứu thường phải chuyển đổi hệ số sang hiệu ứng biên (marginal effects). Ngoài ra, do sử dụng hàm liên kết chuẩn, probit có thể tốn nhiều tài nguyên tính toán hơn logit khi làm việc với dữ liệu lớn.
Smith và cộng sự (2018) đã phân tích dữ liệu lương của kỹ sư phần mềm tại Mỹ bằng mô hình hồi quy tuyến tính đa biến, kết quả cho thấy kinh nghiệm và trình độ học vấn là các yếu tố có ảnh hưởng mạnh nhất đến thu nhập.
Nguyễn và Trần (2020) nghiên cứu thị trường lao động kỹ sư cơ khí ở Việt Nam sử dụng mô hình log-lin, chỉ ra rằng kỹ năng ngoại ngữ và chứng chỉ nghề nghiệp góp phần nâng cao mức lương bên cạnh kinh nghiệm và học vấn.
Lee (2021) áp dụng mô hình probit để dự đoán khả năng đạt mức lương cao của kỹ sư công nghệ tại Hàn Quốc, cho thấy quy mô công ty và cơ hội thăng tiến nội bộ là các yếu tố quyết định quan trọng.
Alvarez và cộng sự (2022) phân tích dữ liệu kỹ sư dân dụng ở châu Âu bằng mô hình hồi quy tuyến tính, kết quả chỉ ra rằng vị trí địa lý và việc tham gia các dự án quốc tế giúp kỹ sư có thu nhập cao hơn đáng kể.
Trước khi phân tích những yếu tố sẽ tác động đến tiền lương của những người kĩ sư chúng ta cần phân nhóm thành ba mức lương là thấp, trung và cao để dễ dàng nhận diện, so sánh. Hơn thế, chúng ta có thể phân tích sâu sắc về cấu trúc tiên lương cũng như là các yếu tố định tính sẽ ảnh hưởng đáng kể đến các mức phân loại tiền lương của các kỹ sư một cách hiệu quả nhất.
d$sal <- cut(d$salary_in_usd,
breaks = c(0, 100000,200000, Inf),
labels = c("thap","trung", "cao"))
Kết quả từ hàm cut() sẽ giúp phân thành ba nhóm lương từ
số liệu trong cột salary_in_usd và lưu trong một cột mới có tên là
sal cụ thể như sau.
Nhóm Data Science: Gồm các vị trí như Data Scientist, Applied Scientist,… chuyên phân tích dữ liệu để tìm ra thông tin có giá trị phục vụ quyết định kinh doanh. Họ sử dụng các kỹ thuật thống kê, mô hình dự đoán và trực quan hóa dữ liệu.
Nhóm Data Engineering: Bao gồm Data Engineer, Data Architect,… với nhiệm vụ xây dựng và quản lý hệ thống thu thập, lưu trữ, xử lý dữ liệu. Họ đảm bảo dữ liệu sạch, sẵn sàng cho phân tích và mô hình.
Nhóm Machine Learning / AI: Gồm Machine Learning Engineer, AI Engineer,… tập trung vào việc phát triển và triển khai các mô hình học máy, học sâu nhằm tối ưu hóa quy trình, sản phẩm thông minh.
Nhóm khác: Là các công việc không thuộc ba nhóm trên, như Data Analyst, Head of Data, hay vai trò quản lý. Những vị trí này thiên về hỗ trợ kỹ thuật, phân tích kinh doanh hoặc điều hành chiến lược dữ liệu.
Trong phần này, em tiến hành lập bảng tần số cho biến job_group nhằm thống kê số lượng kỹ sư theo từng nhóm công việc. Việc này giúp so sánh sự khác biệt về quy mô giữa các nhóm, đồng thời cung cấp cái nhìn trực quan về phân bố của biến job_group trước khi thực hiện các phân tích về tác động của loại công việc đến mức lương của kỹ sư.
# Phân nhóm
d <- d %>%
mutate(job_group = case_when(
str_detect(tolower(job_title), "data engineer|data architect|analytics engineer|data specialist|data modeler") ~ "Data Engineering",
str_detect(tolower(job_title), "data scientist|data science|research scientist|decision scientist") ~ "Data Science",
str_detect(tolower(job_title), "machine learning|ml engineer|ai engineer|applied scientist") ~ "ML/AI",
TRUE ~ "Other"
))
# Đếm số lượng từng nhóm
job_group_counts <- d %>%
count(job_group, sort = TRUE)
# Vẽ biểu đồ
ggplot(job_group_counts, aes(x = reorder(job_group, -n), y = n, fill = job_group)) +
geom_bar(stat = "identity") +
geom_text(aes(label = n), vjust = 1.5, color = "white", fontface = "bold") +
scale_fill_brewer(palette = "Paired") +
labs(title = "Số lượng kỹ sư theo các nhóm công việc",
x = "Nhóm nghề nghiệp", y = "Số lượng") +
theme_minimal() +
theme(axis.text.x = element_text(angle = 30, hjust = 1))
Kết quả bảng tần số và đồ thị cho thấy nhóm công việc Other có số lượng kỹ sư cao nhất với 4.752 người. Tiếp theo là nhóm Data Science với 4.533 kỹ sư và Data Engineering với 4.517 kỹ sư, hai nhóm này có quy mô gần tương đương nhau. Trong khi đó, nhóm ML/AI có số lượng kỹ sư thấp nhất chỉ 2.732 người. Nhìn chung, các nhóm ngoài ba mảng dữ liệu chính chiếm tỷ trọng lớn, còn nhóm ML/AI chiếm tỷ trọng nhỏ nhất trong cơ cấu ngành nghề.
Ở khía cạnh khác, chúng ta có thể phân tích tỷ lệ phần trăm đối với từng nhóm công việc và trực quan chúng bằng đồ thị tròn để có thể dễ dàng so sánh với nhau.
prop.table(table(d$job_group)) * 100
##
## Data Engineering Data Science ML/AI Other
## 27.31946 27.41623 16.52353 28.74078
f <- d %>%
count(job_group, name = "Count") %>%
mutate(Percentage = Count / sum(Count) * 100)
ggplot(f, aes(x = "", y = Percentage, fill = job_group)) +
geom_bar(width = 1, stat = "identity") +
coord_polar("y") +
labs(title = "Tỷ lệ phần trăm các nhóm công việc kỹ sư", x = "", y = "") +
theme_void() +
theme(legend.title = element_blank(),
plot.title = element_text(hjust = 0.5, size = 14, face = "bold")) +
geom_text(aes(label = paste0(round(Percentage, 1), "%")),
position = position_stack(vjust = 0.5), size = 4, color = "white")
Kết quả bảng tần số và biểu đồ cho thấy nhóm công việc Other chiếm tỷ lệ cao nhất với khoảng 28.74% tổng số kỹ sư. Nhóm Data Science và Data Engineering có tỷ lệ gần tương đương nhau lần lượt là 27.42% và 27.32%, thể hiện sự cân bằng giữa hai lĩnh vực này. Trong khi đó, nhóm ML/AI có tỷ lệ thấp nhất chỉ 16.52%, cho thấy quy mô nhân lực trong nhóm này nhỏ hơn đáng kể so với các nhóm công việc khác.
Chúng ta lập bảng tần số giữa hai biến job_group ( nhóm công việc) và sal (mức lương phân tổ thành ba nhóm). Bảng tần số giúp chúng ta hiểu rõ hơn về sự phân bố mức lương theo từng hình thức làm việc.
tab <- table(d$job_group, d$sal)
addmargins(tab)
##
## thap trung cao Sum
## Data Engineering 981 2827 709 4517
## Data Science 807 2694 1032 4533
## ML/AI 249 1427 1056 2732
## Other 2056 2310 386 4752
## Sum 4093 9258 3183 16534
Vẽ đồ thị cột phân theo từng nhóm công việc để đánh giá mức lương của từng nhóm.
table_data <- table(d$job_group, d$sal)
# Chuyển bảng sang định dạng dài (long format)
long_data <- as.data.frame(table_data)
colnames(long_data) <- c("job_group", "salary", "count")
# Vẽ biểu đồ cột theo từng nhóm công việc
ggplot(long_data, aes(x = salary, y = count, fill = salary)) +
geom_bar(stat = "identity", position = position_dodge()) +
geom_text(aes(label = count), vjust = 1.5, color = "white", size = 3, position = position_dodge(1)) +
facet_wrap(~ job_group, scales = "free_y", nrow = 2) +
labs(title = "Phân bố mức lương theo nhóm công việc",
x = "Mức lương",
y = "Số lượng kỹ sư",
fill = "Mức lương") +
theme_minimal()
Dựa vào đồ thị, chúng ta có thể thấy rằng:
Nhóm Data Engineering có tổng số 4.517 kỹ sư, trong đó 981 kỹ sư (21.72%) thuộc mức lương thấp, 2.827 kỹ sư (62.61%) ở mức lương trung và 709 kỹ sư (15.70%) ở mức lương cao.
Nhóm Data Science có 4.533 kỹ sư, với 807 kỹ sư (17.80%) ở mức lương thấp, 2.694 kỹ sư (59.45%) đạt mức lương trung và 1.032 kỹ sư (22.76%) thuộc nhóm lương cao.
Đối với nhóm ML/AI, tổng số kỹ sư là 2.732, trong đó chỉ có 249 kỹ sư (9.12%) ở mức lương thấp. Đa số kỹ sư thuộc nhóm này có mức lương trung với 1.427 người (52.24%) và 1.056 người (38.64%) đạt mức lương cao.
Cuối cùng, nhóm Other có 4.752 kỹ sư, với số lượng ở mức lương thấp là 2.056 kỹ sư (43.26%), mức lương trung chiếm 2.310 kỹ sư (48.62%) và mức lương cao chỉ có 386 kỹ sư (8.12%).
Nhìn chung, nhóm ML/AI có tỷ lệ kỹ sư nhận mức lương cao lớn nhất, trong khi nhóm Other có tỷ lệ lương thấp cao nhất so với các nhóm công việc còn lại.
Trong phần này em sẽ thực hiện tính Relative Risk nhằm phân tích sự ảnh hưởng của từng nhóm công việc lên các mức lương khác nhau.
d <- d %>%
mutate(job_group = recode(job_group,
"Data Engineering" = "DE",
"Data Science" = "DS",
"ML/AI" = "ML",
"Other" = "OT"))
Em sẽ tính toán Relative Risk trong đó biến độc lập là nhóm công việc và phụ thuộc là mức lương. Nhóm công việc DE được chọn làm tham chiếu.
Trước tiên cần tính tỷ lệ mức lương thấp, trung, cao trong tổng số người ở mỗi nhóm công việc khác nhau. Sau đó lấy tỷ lệ vừa tính được chia cho tỷ lệ ở nhóm tham chiếu, trong ví dụ là nhóm làm việc DE
# Tạo bảng tần số nhóm công việc x mức lương
RR_table <- as.data.frame.matrix(table(d$job_group, d$sal))
# Tính tổng số người mỗi nhóm
RR_table$Total <- rowSums(RR_table)
# Tính tỷ lệ mức lương cho từng nhóm công việc
RR_table$Rate_thap <- RR_table$thap / RR_table$Total
RR_table$Rate_trung <- RR_table$trung / RR_table$Total
RR_table$Rate_cao <- RR_table$cao / RR_table$Total
# Lấy tỷ lệ nhóm tham chiếu DE
rate_ref_thap <- RR_table["DE", "Rate_thap"]
rate_ref_trung <- RR_table["DE", "Rate_trung"]
rate_ref_cao <- RR_table["DE", "Rate_cao"]
# Tính RR so với DE
RR_table$RR_thap <- RR_table$Rate_thap / rate_ref_thap
RR_table$RR_trung <- RR_table$Rate_trung / rate_ref_trung
RR_table$RR_cao <- RR_table$Rate_cao / rate_ref_cao
# Chỉ lấy các cột RR và làm tròn
RR_out <- round(RR_table[, c("RR_thap", "RR_trung", "RR_cao")], 2)
# Xuất bảng đẹp bằng kable
kable(RR_out,
col.names = c("RR Thấp", "RR Trung", "RR Cao"),
align = 'c') %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"),
full_width = F) %>%
column_spec(1, bold = TRUE, color = "white", background = "red") %>%
column_spec(2:4, border_left = TRUE, border_right = TRUE)
| RR Thấp | RR Trung | RR Cao | |
|---|---|---|---|
| DE | 1.00 | 1.00 | 1.00 |
| DS | 0.82 | 0.95 | 1.45 |
| ML | 0.42 | 0.83 | 2.46 |
| OT | 1.99 | 0.78 | 0.52 |
Dựa vào bảng Relative Risk, nhóm Data Engineering (DE) được chọn làm nhóm tham chiếu với RR = 1 cho cả ba mức lương.
Nhóm Data Science (DS) có nguy cơ đạt mức lương thấp thấp hơn 18% so với DE (RR = 0.82) và khả năng đạt mức lương cao cao hơn 45% so với DE (RR = 1.45). Điều này cho thấy Data Science có xu hướng dịch chuyển về mức lương cao.
Nhóm ML/AI (ML) có RR mức lương thấp là 0.42, tức khả năng nhận mức lương thấp chỉ bằng 42% so với DE, nhưng RR mức lương cao lại là 2.46 – cao hơn 146%. Đây là nhóm có xác suất đạt mức lương cao nhất trong các nhóm.
Ngược lại, nhóm Other (OT) có RR mức lương thấp là 1.99, tức gần gấp đôi DE, trong khi RR mức lương cao chỉ là 0.52. Điều này cho thấy nhóm Other tập trung nhiều ở mức lương thấp và ít có cơ hội đạt mức lương cao.
Trong phần này em sẽ thực hiện khoảng ước lượng cho RR với khoảng tin cậy 95%. Tuy nhiên để cho kết quả dễ dàng phân tích hơn trước nên em sẽ thực hiện trên bảng có số lượng nhỏ hơn là 2x2.
Trước tiên em tính khoảng ước lượng Relative Risk đối với mức lương thấp và trung nhưng trong nhóm DE và DS.
RR_DE_DS <- table(d$job_group, d$sal)[c("DE","DS"), c("thap","trung")]
# Xem bảng để kiểm tra
RR_DE_DS
##
## thap trung
## DE 981 2827
## DS 807 2694
# Tính Relative Risk với khoảng tin cậy 95%
epitab(RR_DE_DS, method = "riskratio", rev = "c")
## $tab
##
## trung p0 thap p1 riskratio lower upper p.value
## DE 2827 0.7423845 981 0.2576155 1.0000000 NA NA NA
## DS 2694 0.7694944 807 0.2305056 0.8947658 0.8251015 0.9703118 0.00759801
##
## $measure
## [1] "wald"
##
## $conf.level
## [1] 0.95
##
## $pvalue
## [1] "fisher.exact"
Xác suất để có mức lương thấp đối với nhóm Data Engineering là khoảng 25.76%, trong khi đó nhóm Data Science có xác suất đạt mức lương thấp khoảng 23.05%. Tỷ lệ Relative Risk (RR) giữa hai nhóm là 0.895, nghĩa là khả năng kỹ sư thuộc nhóm Data Science nhận mức lương thấp chỉ bằng 0.895 lần so với nhóm Data Engineering. Khoảng ước lượng RR nằm trong khoảng 0.825 đến 0.970 với độ tin cậy 95%. Giá trị p.value là 0.0076 nhỏ hơn mức ý nghĩa 5%, vì vậy sự khác biệt này có ý nghĩa thống kê, cho thấy nhóm Data Science có xu hướng đạt mức lương trung trở lên cao hơn so với nhóm Data Engineering.
Tiếp theo, em tính khoảng ước lượng Relative Risk đối với mức lương thấp và trung nhưng trong trường hợp này sẽ phân tích trên hai nhóm DE và ML.
# Tạo bảng 2x2 cho DE và ML, chỉ lấy mức lương thap & trung
RR_DE_ML <- table(d$job_group, d$sal)[c("DE","ML"), c("thap","trung")]
# Tính RR và khoảng tin cậy 95%
epitab(RR_DE_ML, method = "riskratio", rev = "c")
## $tab
##
## trung p0 thap p1 riskratio lower upper p.value
## DE 2827 0.7423845 981 0.2576155 1.0000000 NA NA NA
## ML 1427 0.8514320 249 0.1485680 0.5767044 0.5080961 0.6545769 6.075812e-20
##
## $measure
## [1] "wald"
##
## $conf.level
## [1] 0.95
##
## $pvalue
## [1] "fisher.exact"
Xác suất để có mức lương thấp đối với nhóm Data Engineering là khoảng 25.76%, trong khi nhóm ML/AI có xác suất đạt mức lương thấp chỉ khoảng 14.86%. Tỷ lệ Relative Risk (RR) giữa hai nhóm là 0.577, nghĩa là khả năng kỹ sư thuộc nhóm ML/AI nhận mức lương thấp chỉ bằng 0.577 lần so với nhóm Data Engineering. Khoảng ước lượng RR nằm trong khoảng 0.508 đến 0.655 với độ tin cậy 95%. Giá trị p.value là 6.08e-20, nhỏ hơn mức ý nghĩa 5%, do đó kết quả này có ý nghĩa thống kê. Điều này cho thấy nhóm kỹ sư ML/AI có xu hướng đạt mức lương trung trở lên cao hơn đáng kể so với nhóm Data Engineering.
RR_DE_OT <- table(d$job_group, d$sal)[c("DE","OT"), c("thap","trung")]
epitab(RR_DE_OT, method = "riskratio", rev = "c")
## $tab
##
## trung p0 thap p1 riskratio lower upper p.value
## DE 2827 0.7423845 981 0.2576155 1.000000 NA NA NA
## OT 2310 0.5290884 2056 0.4709116 1.827963 1.717358 1.945691 1.503235e-89
##
## $measure
## [1] "wald"
##
## $conf.level
## [1] 0.95
##
## $pvalue
## [1] "fisher.exact"
Xác suất để có mức lương thấp đối với nhóm Data Engineering là khoảng 25.76%, trong khi nhóm Other có xác suất đạt mức lương thấp cao hơn đáng kể với 47.09%. Tỷ lệ Relative Risk (RR) được ước lượng là 1.828, cho thấy khả năng kỹ sư trong nhóm Other nhận mức lương thấp cao gấp 1.828 lần so với nhóm Data Engineering. Khoảng ước lượng RR nằm trong khoảng 1.717 đến 1.946 với độ tin cậy 95%. Giá trị p.value là 1.50e-89, nhỏ hơn mức ý nghĩa 5%, chứng tỏ kết quả này có ý nghĩa thống kê. Điều này cho thấy nhóm kỹ sư Other có rủi ro nhận mức lương thấp cao hơn đáng kể so với nhóm Data Engineering.
Odd ratio và khoảng ước lượng
Chúng ta có thể hiểu rằng Odd là một đại lượng thống kê đo lường xác suất của một sự kiện xảy ra so với không xảy ra. Odd được tính bằng cách chia tỷ lệ giữa xác suất xảy ra và xác suất không xảy ra của một sự kiện. Ngoài ra ta có thêm chỉ số odd ratio là tỷ lệ giữa odds của hai nhóm khác nhau. Nó cho biết mức độ thay đổi của odd giữa hai nhóm so với nhau.
Trong phần này em sẽ tính tỷ lệ Odd Ratio của hai biến lương và nhóm công việc. Để dễ dàng trong phân tích em sẽ tách số liệu từ bảng tần số ban đầu thành bảng tần số mới chỉ bao gồm nhóm lương thấp, cao và nhóm công DE và DS.
# Tạo bảng tần số nhóm công việc x mức lương
odd_job <- table(d$job_group, d$sal)
# Lọc chỉ lấy 2 nhóm công việc và mức lương thap/cao
odd2 <- odd_job[c("DE","DS"), c("thap","cao")]
# Tính Odds Ratio (OR)
epitab(odd2, method = "oddsratio")
## $tab
##
## thap p0 cao p1 oddsratio lower upper p.value
## DE 981 0.5486577 709 0.4072372 1.000000 NA NA NA
## DS 807 0.4513423 1032 0.5927628 1.769412 1.548317 2.022078 4.311848e-17
##
## $measure
## [1] "wald"
##
## $conf.level
## [1] 0.95
##
## $pvalue
## [1] "fisher.exact"
Kết quả Odd Ratio của nhóm Data Science (DS) so với Data Engineering (DE) là khoảng 1.769. Điều này có nghĩa là odds đạt mức lương cao của nhóm DS cao hơn khoảng 1.769 lần so với nhóm DE.
Khoảng tin cậy 95% của OR nằm trong khoảng 1.548 đến 2.022, và giá trị p là 4.31e-17, nhỏ hơn 0.05. Điều này chứng tỏ có sự khác biệt có ý nghĩa thống kê giữa odds đạt mức lương cao của hai nhóm DS và DE. Nói cách khác, kỹ sư thuộc nhóm Data Science có khả năng đạt mức lương cao hơn đáng kể so với nhóm Data Engineering.
epitab(odd2, method = "oddsratio", rev = "c")
## $tab
##
## cao p0 thap p1 oddsratio lower upper p.value
## DE 709 0.4072372 981 0.5486577 1.0000000 NA NA NA
## DS 1032 0.5927628 807 0.4513423 0.5651595 0.4945408 0.6458625 4.311848e-17
##
## $measure
## [1] "wald"
##
## $conf.level
## [1] 0.95
##
## $pvalue
## [1] "fisher.exact"
Kết quả Odd Ratio của nhóm Data Science (DS) so với Data Engineering (DE) là khoảng 0.565. Điều này có nghĩa là odds đạt mức lương thấp của nhóm DS chỉ bằng khoảng 56.5% so với nhóm DE.
Khoảng tin cậy 95% của OR nằm trong khoảng 0.495 đến 0.646, và giá trị p là 4.31e-17, nhỏ hơn 0.05. Điều này chứng tỏ có sự khác biệt có ý nghĩa thống kê giữa odds đạt mức lương thấp của hai nhóm DS và DE. Nói cách khác, kỹ sư thuộc nhóm Data Science có khả năng nhận mức lương thấp thấp hơn đáng kể so với nhóm Data Engineering.
# Tạo bảng tần số nhóm công việc x mức lương
odd_job <- table(d$job_group, d$sal)
# Chỉ lấy 2 nhóm công việc DE và DS với mức lương thap/trung
odd_DE_DS <- odd_job[c("DE","DS"), c("thap","trung")]
# Tính Odds Ratio
epitab(odd_DE_DS, method = "oddsratio")
## $tab
##
## thap p0 trung p1 oddsratio lower upper p.value
## DE 981 0.5486577 2827 0.5120449 1.000000 NA NA NA
## DS 807 0.4513423 2694 0.4879551 1.158423 1.040816 1.28932 0.00759801
##
## $measure
## [1] "wald"
##
## $conf.level
## [1] 0.95
##
## $pvalue
## [1] "fisher.exact"
Kết quả Odd Ratio của nhóm Data Science (DS) so với Data Engineering (DE) là khoảng 1.158. Điều này cho thấy odds đạt mức lương trung của nhóm DS cao hơn khoảng 15.8% so với nhóm DE.
Khoảng tin cậy 95% của OR nằm trong khoảng 1.041 đến 1.289 và giá trị p là 0.0076, nhỏ hơn 0.05. Điều này chứng tỏ có sự khác biệt có ý nghĩa thống kê giữa odds đạt mức lương trung của hai nhóm. Nói cách khác, kỹ sư thuộc nhóm Data Science có khả năng đạt mức lương trung cao hơn đáng kể so với nhóm Data Engineering.
epitab(odd_DE_DS, method = "oddsratio", rev = "c")
## $tab
##
## trung p0 thap p1 oddsratio lower upper p.value
## DE 2827 0.5120449 981 0.5486577 1.0000000 NA NA NA
## DS 2694 0.4879551 807 0.4513423 0.8632424 0.7756026 0.960785 0.00759801
##
## $measure
## [1] "wald"
##
## $conf.level
## [1] 0.95
##
## $pvalue
## [1] "fisher.exact"
Kết quả Odd Ratio của nhóm Data Science (DS) so với Data Engineering (DE) là khoảng 0.863. Điều này cho thấy odds đạt mức lương thấp của nhóm DS thấp hơn khoảng 13.7% so với nhóm DE.
Khoảng tin cậy 95% của OR nằm trong khoảng 0.776 đến 0.961 và giá trị p là 0.0076, nhỏ hơn 0.05. Điều này cho thấy sự khác biệt có ý nghĩa thống kê, chứng minh rằng nhóm Data Science có khả năng nhận mức lương thấp thấp hơn đáng kể so với nhóm Data Engineering.
Kết quả bên dưới được sử dụng để kiểm tra tính độc lập giữa hai biến nhóm công việc và mức lương bằng phương pháp chi bình phương. Giả thuyết Ho được đưa ra là biến mức lương và loại công việc độc lập.
chisq.test(table(d$job_group, d$sal))
##
## Pearson's Chi-squared test
##
## data: table(d$job_group, d$sal)
## X-squared = 2029.5, df = 6, p-value < 2.2e-16
Qua kết quả kiểm định cho ta thấy giá trị p−value là 2.2e−16 nhỏ hơn 0.05 , nên bác bỏ H0, nghĩa là mức lương và nhóm công việc là có liên quan với nhau.
Ước lượng tỷ lệ cho nhóm lương trung bình
Trước tiên ta cần lọc kết quả về mức lương theo nhóm DE. Kết quả thu được như sau.
DE <- table(d[d$job_group == "DE", ]$sal)
DE
##
## thap trung cao
## 981 2827 709
prop.test(DE["trung"], sum(DE), p = 0.7)
##
## 1-sample proportions test with continuity correction
##
## data: DE["trung"] out of sum(DE), null probability 0.7
## X-squared = 117.89, df = 1, p-value < 2.2e-16
## alternative hypothesis: true p is not equal to 0.7
## 95 percent confidence interval:
## 0.6115334 0.6399668
## sample estimates:
## p
## 0.6258579
Sau đó ước lượng tỷ lệ với giả thuyết không (H0): Tỷ lệ người có mức lương trung trong nhóm công việc DE là 0.7 (70%).
Kết quả từ kiểm định cho thấy giá trị p-value là 2.2e-16, rất gần với 0 và nhỏ hơn giá trị α=0.05. Điều này cho thấy có bằng chứng mạnh mẽ để bác bỏ giả thuyết rằng tỷ lệ người có mức lương trung trong nhóm công viêc DE là 70%.
Kết quả cũng cho thấy tỷ lệ kỹ sư có mức lương trung trong nhóm công việc DE khoảng 62.59%. Khoảng ước lượng tỷ lệ ở độ tin cây 95% nằm trong khoảng từ 61.15% đến 63.99%.
Ước lượng tỷ lệ cho nhóm lương thấp
prop.test(DE["thap"], sum(DE), p = 0.2)
##
## 1-sample proportions test with continuity correction
##
## data: DE["thap"] out of sum(DE), null probability 0.2
## X-squared = 8.2251, df = 1, p-value = 0.004132
## alternative hypothesis: true p is not equal to 0.2
## 95 percent confidence interval:
## 0.2052898 0.2295544
## sample estimates:
## p
## 0.2171795
Ta có giả thuyết không (H0): Tỷ lệ người có mức lương thấp trong nhóm công việc DE là 0.2 (20%). Vì giá trị p =0.004132 < 0.05, chúng ta bác bỏ giả thuyết không (H0). Điều này có nghĩa là tỷ lệ kỹ sư có mức lương “thấp” trong nhóm công việc DE không bằng 20%. Kết quả cho rằng tỷ lệ này đạt khoảng 21.71% và ở mức độ tin cậy 95% kết quả ước lượng chỉ ra tỷ lệ kỹ sư có mức lương thấp trong nhóm công việc DE đạt khoảng từ 23.88% đến 25.20%.
Ước lượng tỷ lệ cho nhóm lương cao
prop.test(DE["cao"], sum(DE), p = 0.3)
##
## 1-sample proportions test with continuity correction
##
## data: DE["cao"] out of sum(DE), null probability 0.3
## X-squared = 439.4, df = 1, p-value < 2.2e-16
## alternative hypothesis: true p is not equal to 0.3
## 95 percent confidence interval:
## 0.1465387 0.1679755
## sample estimates:
## p
## 0.1569626
Ta có giả thuyết không (H0) là tỷ lệ kỹ sư có mức lương “cao” trong nhóm công việc DE là 30%. Kết quả bên trên có giá trị p rất nhỏ (p < 2.2e-16), chúng ta bác bỏ giả thuyết không (H0). Điều này có nghĩa là tỷ lệ kỹ sư có mức lương “cao” trong nhóm nhóm công việc DE sẽ không bằng 30%. Thực tế, tỷ lệ kỹ sư có mức lương “cao” chỉ khoảng 15.69%, thấp hơn so với tỷ lệ kỳ vọng là 30%. Ngoài ra, khoảng tin cậy 95% cho tỷ lệ kỹ sư có mức lương cao trong nhóm DE là [0.1465, 0.1679]. Điều này có nghĩa là chúng ta có 95% tin tưởng rằng tỷ lệ này nằm trong khoảng từ 14.65% đến 16.79%
Mô hình xác suất tuyến tính
Phần này sẽ thực hiện hồi quy xác suất tuyến tính giữa biến độc lập là nhóm công việc và biến phụ thuộc là lương. Tuy nhiên trong phần này ta cần thay đổi biến sal thành biến nhị phân với giá trị 1 tương ứng với mức lương trung và cao, ngược lại giá trị 0 tương ứng với mức lương thấp.
d$sal2 <- ifelse(d$sal %in% c("trung", "cao"), 1, 0)
d$job_group <- relevel(as.factor(d$job_group), ref = "DE")
glm_job <- glm(sal2 ~ job_group, data = d, family = binomial)
summary(glm_job)
##
## Call:
## glm(formula = sal2 ~ job_group, family = binomial, data = d)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.28218 0.03609 35.532 <2e-16 ***
## job_groupDS 0.24759 0.05301 4.671 3e-06 ***
## job_groupML 1.01759 0.07563 13.456 <2e-16 ***
## job_groupOT -1.01117 0.04647 -21.760 <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: 18506 on 16533 degrees of freedom
## Residual deviance: 17143 on 16530 degrees of freedom
## AIC: 17151
##
## Number of Fisher Scoring iterations: 4
Dựa trên mô hình: \[ \hat\pi = 1.28218 + 0.24759 \cdot DS + 1.01759 \cdot ML - 1.01117 \cdot OT \]
Tất cả các hệ số đều có p-value < 0.05, cho thấy các nhóm công việc đều có ảnh hưởng có ý nghĩa thống kê đến khả năng đạt mức lương trung/cao.
Hệ số \(\hat\beta_0 = 1.28218\): Đây là log-odds của nhóm DE. Khi chuyển sang xác suất: \[ \hat p_{DE} = \frac{e^{1.28218}}{1 + e^{1.28218}} \approx 78.3\% \]
Hệ số \(\hat\beta_{DS} = 0.24759\): So với DE, nhóm DS có odds đạt lương trung/cao cao hơn khoảng 28%. Xác suất ước lượng: \[ \hat p_{DS} = \frac{e^{1.28218+0.24759}}{1 + e^{1.28218+0.24759}} \approx 82.6\% \]
Hệ số \(\hat\beta_{ML} = 1.01759\): Nhóm ML có odds đạt lương trung/cao gấp khoảng 2.77 lần DE. Xác suất ước lượng: \[ \hat p_{ML} = \frac{e^{1.28218+1.01759}}{1 + e^{1.28218+1.01759}} \approx 92.2\% \]
Hệ số \(\hat\beta_{OT} = -1.01117\): Nhóm OT có odds đạt lương trung/cao chỉ bằng khoảng 36% DE. Xác suất ước lượng: \[ \hat p_{OT} = \frac{e^{1.28218-1.01117}}{1 + e^{1.28218-1.01117}} \approx 53.2\% \]
Kết luận:
Nhóm ML có tỷ lệ đạt mức lương trung/cao cao nhất
(~92%), tiếp theo là DS (~82%) và DE
(~78%).
Ngược lại, nhóm OT có tỷ lệ thấp nhất (~53%), thể hiện
mối quan hệ âm rõ rệt với khả năng đạt lương trung/cao.
Tất cả các nhóm công việc đều có ảnh hưởng thống kê rõ rệt (p <
0.001).
n_job <- data.frame(job_group = c("DE", "DS", "ML", "OT"))
# Dự đoán xác suất đạt lương trung/cao
predict.glm(glm_job, newdata = n_job, type = "response")
## 1 2 3 4
## 0.7828205 0.8219722 0.9088580 0.5673401
Dựa trên kết quả mô hình logistic, em dự báo xác suất để kỹ sư đạt mức lương trung đến cao đối với từng nhóm công việc như sau:
Nhóm DE (1) có xác suất dự báo đạt mức lương trung/cao khoảng 78.28%. Điều này cho thấy kỹ sư thuộc nhóm DE có cơ hội khá cao để đạt được mức lương từ trung trở lên.
Nhóm DS (2) có xác suất dự báo khoảng 82.20%, cao hơn một chút so với nhóm DE. Điều này gợi ý rằng kỹ sư thuộc nhóm DS có khả năng nhận lương trung/cao cao hơn so với nhóm DE.
Nhóm ML (3) có xác suất dự báo đạt mức lương trung/cao khoảng 90.89%, là nhóm có tỷ lệ cao nhất trong bốn nhóm công việc. Điều này chứng tỏ kỹ sư làm việc trong mảng ML có lợi thế lớn trong việc đạt mức lương cao hơn.
Nhóm OT (4) có xác suất dự báo khoảng 56.73%, thấp nhất trong bốn nhóm. Kết quả này chỉ ra rằng các kỹ sư thuộc nhóm công việc khác (OT) có cơ hội nhận mức lương trung/cao thấp hơn đáng kể so với các nhóm còn lại.
Kết luận: Các nhóm công việc liên quan trực tiếp đến dữ liệu và AI (DS, ML) có tỷ lệ đạt mức lương trung/cao cao hơn rõ rệt, trong khi nhóm công việc khác (OT) có tỷ lệ thấp nhất.
Mô hình logit
Ngoài ra chúng ta có thể thực hiện mô hình hồi quy khi có hàm liên kết dạng \(g(μ)=logit(μ)=ln(\frac{μ}{1−μ})\) khi phân tích mối quan hệ giữa nhóm công việc với mức lương, trong đó nhóm công việc DE sẽ làm nhóm cơ sở.
logit_job <- glm(factor(sal2) ~ job_group, data = d, family = binomial(link = 'logit'))
summary(logit_job)
##
## Call:
## glm(formula = factor(sal2) ~ job_group, family = binomial(link = "logit"),
## data = d)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.28218 0.03609 35.532 <2e-16 ***
## job_groupDS 0.24759 0.05301 4.671 3e-06 ***
## job_groupML 1.01759 0.07563 13.456 <2e-16 ***
## job_groupOT -1.01117 0.04647 -21.760 <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: 18506 on 16533 degrees of freedom
## Residual deviance: 17143 on 16530 degrees of freedom
## AIC: 17151
##
## Number of Fisher Scoring iterations: 4
Kết quả từ mô hình logit:
\[ \ln\left(\frac{\hat\pi}{1-\hat\pi}\right) = 1.28218 + 0.24759 \cdot DS + 1.01759 \cdot ML - 1.01117 \cdot OT \]
Hệ số chặn \(\hat\beta_0 = 1.28218\): Khi không có tác động của các nhóm công việc khác, giá trị logit của kỹ sư thuộc nhóm DE để đạt mức lương trung/cao là 1.28218.
Hệ số \(\hat\beta_{DS} = 0.24759\): Chênh lệch logit giữa nhóm DS và DE là +0.24759.
Hệ số \(\hat\beta_{ML} = 1.01759\): Chênh lệch logit giữa nhóm ML và DE là +1.01759, cho thấy nhóm ML có khả năng đạt mức lương trung/cao cao hơn rõ rệt.
Hệ số \(\hat\beta_{OT} = -1.01117\): Chênh lệch logit giữa nhóm OT và DE là −1.01117, nghĩa là nhóm OT có khả năng đạt mức lương trung/cao thấp hơn đáng kể.
Tất cả các hệ số đều có p-value nhỏ hơn mức ý nghĩa \(\alpha = 5\%\), cho thấy các nhóm công việc có ảnh hưởng thống kê rõ rệt đến mức lương.
n_job <- data.frame(job_group = c("DE", "DS", "ML", "OT"))
# Dự báo xác suất mức lương trung/cao
predict.glm(logit_job, newdata = n_job, type = "response")
## 1 2 3 4
## 0.7828205 0.8219722 0.9088580 0.5673401
Dựa vào kết quả mô hình, xác suất dự báo cho kỹ sư thuộc nhóm DE là khoảng 78.28%. Điều này có nghĩa kỹ sư trong nhóm DE có khoảng 78% cơ hội đạt được mức lương trung đến cao.
Đối với nhóm DS, xác suất dự báo đạt mức lương trung trở lên là khoảng 82.20%, cao hơn một chút so với nhóm DE, cho thấy công việc trong nhóm DS có khả năng mang lại mức lương cao hơn.
Trong nhóm ML, xác suất dự báo để đạt mức lương trung/cao là khoảng 90.89%. Đây là nhóm có xác suất cao nhất trong bốn nhóm công việc, phản ánh lợi thế rõ rệt về thu nhập của nhóm này.
Cuối cùng, nhóm OT có xác suất đạt mức lương trung trở lên khoảng 56.73%, thấp hơn đáng kể so với các nhóm còn lại. Điều này cho thấy kỹ sư thuộc nhóm công việc khác có cơ hội đạt mức lương trung/cao thấp hơn.
Mô hình Probit
Thay vì sử dụng mô hình logit ta có thể thay thế bằng mô hình probit dự đoán xác suất của một sự kiện xảy ra, mô hình probit sử dụng hàm phân phối chuẩn. Trong đó hàm liên kết có dạng \(probit(π) = Φ^{-1}(π)\) và hàm xác suất là \(π(x)=Φ(β_0+β_1X+β_2X_2+...+β_nX_n)\). Ta thực hiện thao tác khá tương tự với phần trên
# Hồi quy probit với nhóm công việc
probit_job <- glm(factor(sal2) ~ job_group, data = d, family = binomial(link = 'probit'))
# Xem kết quả
summary(probit_job)
##
## Call:
## glm(formula = factor(sal2) ~ job_group, family = binomial(link = "probit"),
## data = d)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.78175 0.02087 37.450 < 2e-16 ***
## job_groupDS 0.14115 0.03018 4.676 2.92e-06 ***
## job_groupML 0.55200 0.03955 13.958 < 2e-16 ***
## job_groupOT -0.61215 0.02774 -22.064 < 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: 18506 on 16533 degrees of freedom
## Residual deviance: 17143 on 16530 degrees of freedom
## AIC: 17151
##
## Number of Fisher Scoring iterations: 4
Ta có kết quả mô hình hồi quy như sau:
\[ \hat\pi = \Phi(0.78175 + 0.14115 \cdot DS + 0.55200 \cdot ML - 0.61215 \cdot OT) \]
Hệ số chặn là 0.78175 và có ý nghĩa thống kê khi giá trị p-value nhỏ hơn mức ý nghĩa 5%. Điều này có nghĩa là khi không có tác động của các nhóm công việc khác, xác suất để kỹ sư thuộc nhóm DE đạt mức lương trung trở lên là \(\Phi(0.78175)\).
Hệ số \(\hat\beta_{DS} = 0.14115\) cho thấy trong trường hợp các yếu tố khác không đổi, chênh lệch tỷ lệ đạt mức lương trung/cao giữa nhóm DS và DE vào khoảng \(\Phi(0.14115)\).
Hệ số \(\hat\beta_{ML} = 0.55200\) có ý nghĩa rằng nhóm ML có chênh lệch tỷ lệ đạt mức lương trung/cao so với DE là \(\Phi(0.55200)\), phản ánh lợi thế về thu nhập của nhóm này.
Hệ số \(\hat\beta_{OT} = -0.61215\) cho thấy nhóm OT có chênh lệch tỷ lệ đạt mức lương trung/cao thấp hơn so với DE, với mức chênh lệch khoảng \(\Phi(-0.61215)\).
Từ kết quả trên em sẽ thực hiện dự báo cho tỷ lệ để đạt được mức lương trung trong các nhóm công việc. Kết quả được trình bày trong bảng bên dưới.
n_job <- data.frame(job_group = c("DE", "DS", "ML", "OT"))
# Dự báo xác suất mức lương trung/cao với mô hình probit
predict.glm(probit_job, newdata = n_job, type = "response")
## 1 2 3 4
## 0.7828205 0.8219722 0.9088580 0.5673401
Dựa vào kết quả ta có thể nhận thấy xác suất dự báo cho kỹ sư thuộc nhóm DE (1) là khoảng 78.28%. Điều này có nghĩa kỹ sư trong nhóm DE có khoảng 78% cơ hội đạt được mức lương trung đến cao.
Bên cạnh đó, đối với nhóm kỹ sư thuộc DS (2), xác suất để đạt được mức lương trung trở lên là khoảng 82.20%, cao hơn một chút so với DE và phản ánh lợi thế thu nhập của nhóm này.
Trong nhóm ML (3), xác suất để đạt mức lương trung/cao là khoảng 90.89%, đây là nhóm có xác suất cao nhất, cho thấy lợi thế rõ rệt về thu nhập.
Cuối cùng, nhóm OT (4) có xác suất đạt mức lương trung trở lên chỉ khoảng 56.73%, thấp hơn đáng kể so với các nhóm còn lại, cho thấy cơ hội đạt mức lương cao trong nhóm này là hạn chế.
Chỉ số AIC
AIC (Akaike Information Criterion) là chỉ số được sử dụng để đánh giá và so sánh các mô hình thống kê. Mô hình có AIC thấp hơn được xem là cân bằng tốt hơn giữa độ phù hợp và độ phức tạp, do đó được ưu tiên lựa chọn.
Em so sánh ba mô hình đã xây dựng (glm, logit, probit) bằng chỉ số AIC:
AIC(glm_job )
## [1] 17150.72
AIC(logit_job )
## [1] 17150.72
AIC(probit_job )
## [1] 17150.72
Kết quả so sánh cho thấy chỉ số AIC của ba mô hình glm, logit và probit đều bằng 17150.72. Điều này dẫn tới các kết luận:
Hệ số Brier Score
Phần này sẽ tính chỉ số Brier Score được sử dụng để đo lường độ chính xác của các dự đoán xác suất. Nó tính toán độ lệch bình phương giữa dự đoán xác suất và giá trị thực tế. Một giá trị Brier Score thấp cho thấy mô hình dự đoán tốt hơn. Vì thế ta sẽ tính và so sánh giá trì này ở hai mô hình hồi quy logit và probit.
BrierScore(logit_job)
## [1] 0.1708017
BrierScore(probit_job)
## [1] 0.1708017
Kết quả bên trên cho thấy Brier Score cho logit và probit đều là 0.1708017. Điều này có nghĩa là cả hai mô hình đều có hiệu suất tương đương khi dự đoán xác suất cho biến mục tiêu. Brier Score khá thấp (gần 0), cho thấy rằng cả hai mô hình đều dự đoán khá tốt và không sai lệch nhiều so với giá trị thực tế.
Ma trận nhầm lẫn
Conf(table(
predict(glm_job) >= 0.5,
glm_job$model$sal2 == '1'
))
##
## Confusion Matrix and Statistics
##
## Reference
## Prediction FALSE TRUE
## FALSE 2056 2696
## TRUE 2037 9745
##
## Total n : 16'534
## Accuracy : 0.7137
## 95% CI : (0.7068, 0.7206)
## No Information Rate : 0.7524
## P-Value [Acc > NIR] : 1.0000
##
## Kappa : 0.2710
## Mcnemar's Test P-Value : < 2.2e-16
##
## Sensitivity : 0.5023
## Specificity : 0.7833
## Pos Pred Value : 0.4327
## Neg Pred Value : 0.8271
## Prevalence : 0.2476
## Detection Rate : 0.2874
## Detection Prevalence : 0.1243
## Balanced Accuracy : 0.6428
## F-val Accuracy : 0.4649
## Matthews Cor.-Coef : 0.2724
##
## 'Positive' Class : FALSE
Từ kết quả ma trận nhầm lẫn của mô hình xác suất tuyến tính cho nhóm công việc, mô hình dự đoán đúng 2,056 trường hợp thuộc mức lương thấp và 9,745 trường hợp thuộc mức lương trung hoặc cao. Tuy nhiên, có 2,696 trường hợp thuộc mức lương trung hoặc cao bị dự đoán nhầm thành lương thấp, và 2,037 trường hợp lương thấp bị dự đoán thành lương trung hoặc cao.
Độ chính xác (accuracy) của mô hình là 0.7137, với khoảng tin cậy 95% nằm trong khoảng từ 0.7068 đến 0.7206. Giá trị sensitivity (khả năng nhận diện đúng mức lương thấp) là 0.5023, trong khi specificity (khả năng nhận diện đúng mức lương trung hoặc cao) là 0.7833. Điều này cho thấy mô hình có xu hướng dự đoán tốt hơn cho nhóm lương trung hoặc cao, nhưng khả năng nhận diện nhóm lương thấp còn hạn chế.
Đầu tiên, em sẽ lập bảng tần số cho biến sal nhằm mục đích để nắm bắt được số lượng kỹ sư trong mỗi nhóm lương được chia ở phần trên. Kết quả này cung cấp một cái nhìn tổng quan về cách mức lương được phân phối trong dữ liệu.
table(d$sal) %>% addmargins()
##
## thap trung cao Sum
## 4093 9258 3183 16534
Chúng ta có thể thực hiện vẽ biểu đồ cột thể hiện số lượng kỹ sư ở mỗi nhóm lương khác nhau để trực quan hóa số liệu từ bảng số liệu phía trên. Điều này giúp ta dễ dàng nhận xét và so sánh mỗi nhóm lương với nhau.
s <- as.data.frame(table(d$sal))
colnames(s) <- c("Salary", "Count")
ggplot(s, aes(x = Salary, y = Count, fill = Salary)) +
geom_bar(stat = "identity") +
geom_text(aes(label = Count), vjust = -0.3, size = 3.5) +
theme_minimal()
Kết quả bảng kết quả trên và bảng tần số cho thấy số lượng kỹ sư có mức lương thấp (từ 0 đến 100,000 USD) là 4093 người. Ở nhóm người có mức lương trung bình tức là từ 100,000 đến 200,000 USD trên năm đạt khoảng 9258 người, chiếm cao nhất trong tổng số người khảo sát. Cuối cùng số lượng kỹ sư có mức lương cao tức trên 200,000 USD là 3183 người, có số lượng thấp nhất trong các nhóm. Tổng số người được khảo sát ở 3 nhóm đạt khoảng 16534 người.
Ngoài ra chúng ta sẽ thực hiện tính tỷ lệ phần trăm cho từng nhóm và vẽ đồ thị cột của chúng để so sánh theo số liệu tương đối
prop.table(table(d$sal))
##
## thap trung cao
## 0.2475505 0.5599371 0.1925124
s$Percentage <- (s$Count / sum(s$Count)) * 100
ggplot(s, aes(x = "", y = Percentage, fill = Salary)) +
geom_bar(width = 1, stat = "identity") +
coord_polar("y") +
labs(title = "Tỷ lệ phần trăm", x = "", y = "") +
theme_void() +
theme(legend.title = element_blank()) +
geom_text(aes(label = paste(round(Percentage, 2), "%")),
position = position_stack(vjust = 0.5), size = 4)
Dựa vào độ thị tròn và bảng kết quả trên ta có thể thấy nhóm chiếm tỷ lệ cao nhất trong tổng số người là nhóm mức lương trung đạt khoảng 55,99%. Tiếp theo đó là nhóm có mức lương thấp đạt khoảng 24,76% trên tổng số người. Cuối cùng là nhóm có mức lương cao đạt tỷ lệ thấp nhất với 19,25%.
Chênh lệch số kỹ sư có mức lương trung và thấp là 5165 người đạt với tỷ lệ đạt khoảng 28,23%. Mặt khác chênh lệch giữa nhóm lương trung và cao là 6075 với tỷ lệ chênh lệch là 36,74% điều này cho thấy hầu hết các kỹ sư đều có mức lương thuộc nhóm trung và có sự khác biệt lớn giữa các nhóm khác so với nhóm này. Ngoài ra số kỹ sư có mức lương thấp nhiều hơn nhóm mức lương cao khoảng 910 người với đạt tỷ lệ 5.51% trên tổng số.
Trong phần này em lập bảng tần số cho biến employment_type nhằm mục tiêu thống kê số lượng kỹ sư trong mỗi nhóm hình thức công việc và so sánh sự khác biệt về số lượng cũng như tỷ lệ giữa chúng. Từ đó giúp cho mọi người hiểu rõ hơn ý nghĩa của biến này trước khi đánh giá sự ảnh hưởng của hình thức công việc đến mức lương.
Trước tiên chúng ta cần đổi tên biến từ “employment_type” thành “type” điều này sẽ giúp câu lệnh được thực hiện một cách ngắn gọn và dễ nhận biết hơn trước.
d <- d %>% rename(type = employment_type) # đổi tên biến
Tiếp đến, em sẽ lập bảng tần số cho biến type, kết quả sẽ được trình bày ở dưới đây.
table(d$type)
##
## CT FL FT PT
## 28 14 16454 38
# Tạo bảng tần số cho biến type
t <- as.data.frame(table(d$type))
colnames(t) <- c("Type", "Count")
ggplot(t, aes(x = Type, y = Count, fill = Type)) +
geom_bar(stat = "identity", width = 0.5) +
geom_text(aes(label = Count), vjust = -0.5, size = 3.5) + # Hiển thị số liệu bên trên cột
labs(title = "DO THI LOAI CONG VIEC ", x = "Type", y = "Frequency") +
theme_minimal()
Dựa vào kết quả từ bảng và đồ thị ta nhận thấy hầu hết các kỹ sư sẽ làm việc theo hình thức toàn thời gian (FT) với số lượng 16454 người. Ngược lại, thì số lượng kỹ sư ở các nhóm bán thời gian (PT) là 38 người, tiếp theo là nhóm làm việc theo hợp đồng có 28 người. Cuối cùng nhóm làm việc tự do có số lượng thấp nhất là 14 người. Nhìn chung sự chênh lệch giữa số lượng kỹ sư giữa nhóm toàn thời gian với các nhóm khác là vô cùng lớn, chẳng hạn nhóm này nhiều hơn nhóm làm việc tự do là 16440 người
Mặt khác ta cũng cần thực hiện phân tích tỷ lệ phần trăm đối với từng nhóm hình thức công việc và trực quan chúng bằng đồ thị tròn để có thể dễ dàng so sánh với nhau.
prop.table(table(d$type))*100
##
## CT FL FT PT
## 0.16934801 0.08467401 99.51614854 0.22982944
t$Percentage <- (t$Count / sum(t$Count)) * 100
ggplot(t, aes(x = "", y = Percentage, fill = Type)) +
geom_bar(width = 1, stat = "identity") +
coord_polar("y") +
labs(title = "Tỷ lệ phần trăm", x = "", y = "") +
theme_void() +
theme(legend.title = element_blank()) +
geom_text(aes(label = paste(round(Percentage, 2), "%")),
position = position_stack(vjust = 0.5), size = 4)
Trong trường hợp này, chúng ta lập bảng tần số giữa hai biến type (hình thức công việc) và sal (mức lương phân tổ thành ba nhóm). Bảng tần số giúp chúng ta hiểu rõ hơn về sự phân bố mức lương theo từng hình thức công việc.
table(d$type, d$sal) %>% addmargins()
##
## thap trung cao Sum
## CT 16 9 3 28
## FL 14 0 0 14
## FT 4037 9239 3178 16454
## PT 26 10 2 38
## Sum 4093 9258 3183 16534
Vẽ đồ thị cột phân theo từng nhóm loại hình công việc để đánh giá mức lương của từng nhóm.
table_data <- table(d$type, d$sal)
# Chuyển bảng chéo sang định dạng dài
long_data <- as.data.frame(table_data)
colnames(long_data) <- c("type", "salary", "count")
# Vẽ biểu đồ cột
ggplot(long_data, aes(x = salary, y = count, fill = salary)) +
geom_bar(stat = "identity", position = position_dodge()) +
geom_text(aes(label = count), vjust = 1, size = 2.5, position = position_dodge(1)) +
facet_wrap(~ type, scales = "free_y", nrow = 2) +
labs(title = "Phân bố mức lương theo hình thức công việc",
x = "Mức lương",
y = "Số lượng kỹ sư",
fill = "Mức lương") +
theme_minimal()
Dựa vào số liệu và đồ thị em đưa ra một số nhận xét như sau:
Như phần trên đã đề cập đến kỹ sư toàn thời gian (FT) chiếm đa số tuyệt đối với 16454 kỹ sư, chiếm khoảng 99.5% tổng số kỹ sư. Ngược lại các hình thức công việc khác chiếm rất nhỏ gần như không đáng kể.
kỹ sư toàn thời gian (FT) có sự phân bố đều giữa các nhóm lương: 4037 kỹ sư có lương thấp, 9239 kỹ sư có lương trung bình, và 3178 kỹ sư có lương cao.
kỹ sư theo hình thức hợp đồng (CT) chủ yếu có mức lương thấp và trung bình, với rất ít kỹ sư có mức lương cao. Trong đó 16 kỹ sư có lương thấp, 9 kỹ sư có lương trung bình, và 3 kỹ sư có lương cao
kỹ sư làm việc tự do (FL) chỉ có mặt trong nhóm lương thấp gồm 14 người
kỹ sư bán thời gian (PT) chủ yếu có mức lương thấp có 26 người, với một số ít trong nhóm lương trung bình gồm có 10 người và rất ít người trong nhóm lương cao chỉ có 2 người.
Tiếp theo, lập bảng tần suất giữa hai biến lương và hình thức công việc. Từ đó vẽ đồ thì cột theo từng nhóm công việc như bên dưới đây.
table_data <- round(prop.table(table_data)*100,2)
table_data
##
## thap trung cao
## CT 0.10 0.05 0.02
## FL 0.08 0.00 0.00
## FT 24.42 55.88 19.22
## PT 0.16 0.06 0.01
# Chuyển đổi dữ liệu thành dạng dài để dễ vẽ biểu đồ
table_data <- as.data.frame.table(table_data)
colnames(table_data) <- c("Group", "Level", "Percentage")
# Vẽ biểu đồ
ggplot(table_data, aes(x = Level, y = Percentage, fill = Level)) +
geom_bar(stat = "identity", position = "dodge") +
geom_text(aes(label = Percentage), vjust = 1, size = 2.5, position = position_dodge(1)) +
facet_wrap(~ Group, scales = "free_y", nrow = 2) +
labs(title = "Đồ thị phần trăm theo mức độ và nhóm công việc",
x = "Nhóm công việc",
y = "Phần trăm",
fill = "Mức độ") +
scale_fill_manual(values = c("thap" = "blue", "trung" = "orange", "cao" = "red")) +
theme_minimal()
Đối với nhóm làm việc theo hợp đồng: Tỷ lệ mức lương thấp chiếm đạt cao nhất khoảng 0.1% so với tổng số kỹ sư, tỷ lệ kỹ sư mức trung là 0.05% và thấp nhất là tỷ lệ mức lương cao 0.02%
Nhóm làm việc tự do ta nhận thấy có tỷ lệ kỹ sư có mức lương thấp đạt 0.08% trên tổng số kỹ sư được khảo sát. Nhưng tỷ lệ mức trung và cao là 0%.
Trong nhóm kỹ sư toàn thời gian chiếm cao nhất là tỷ lệ mức lương trung đạt 55.88% so với tổng số kỹ sư, tiếp đến là tỷ lệ mức lương thấp đạt khoảng 24.42% và thấp nhất là tỷ lệ lương cao trên tổng số đạt khoảng 19.22%. Khi so sánh với nhóm công việc khác thì nhóm làm việc toàn thời gian luôn có tỷ lệ cao hơn ở tất cả mức lương.
Trong nhóm bán thời gian tỷ lệ mức lương thấp chiếm đạt cao nhất khoảng 0.16% so với tổng số kỹ sư, tỷ lệ kỹ sư mức trung là 0.06% và thấp nhất là tỷ lệ mức lương cao 0.01%
Trong phần này em sẽ thực hiện tính Relative Risk nhằm phân tích sự ảnh hưởng của từng loại hình công việc lên các mức lương khác nhau.
RR <- matrix(c(16, 9, 3,
14, 0, 0,
4037, 9239, 3178,
26, 10, 2),
nrow = 4,
byrow = TRUE,
dimnames = list(c("CT", "FL", "FT", "PT"),
c("thap", "trung", "cao")))
RR
## thap trung cao
## CT 16 9 3
## FL 14 0 0
## FT 4037 9239 3178
## PT 26 10 2
Dựa trên bảng tần số đã được giải thích rõ ở phần trước em sẽ tính toán Relative Risk trong đó biến độc lập là loại hình công việc và phụ thuộc là mức lương. Nhóm công việc được chọn làm tham chiếu là FT
Trước tiên cần tính tỷ lệ mức lương thấp, trung, cao trong tổng số người ở mỗi nhóm công việc khác nhau. Sau đó lấy tỷ lệ vừa tính được chia cho tỷ lệ ở nhóm tham chiếu, trong ví dụ là nhóm làm việc toàn thời gian FT.
RR_table <- as.data.frame(RR)
# Tính tổng số người trong mỗi nhóm
RR_table$Total <- rowSums(RR_table)
# Tính tỷ lệ người ở mỗi mức lương trong mỗi nhóm
RR_table$Rate_thap <- RR_table$thap / RR_table$Total
RR_table$Rate_trung <- RR_table$trung / RR_table$Total
RR_table$Rate_cao <- RR_table$cao / RR_table$Total
# Tính tỷ lệ cho nhóm tham chiếu (FT)
rate_ref_thap <- RR_table["FT", "Rate_thap"]
rate_ref_trung <- RR_table["FT", "Rate_trung"]
rate_ref_cao <- RR_table["FT", "Rate_cao"]
# Tính Relative Risk cho các nhóm so với nhóm tham chiếu (FT)
RR_table$RR_thap <- RR_table$Rate_thap / rate_ref_thap
RR_table$RR_trung <- RR_table$Rate_trung / rate_ref_trung
RR_table$RR_cao <- RR_table$Rate_cao / rate_ref_cao
| RR Thap | RR Trung | RR Cao | |
|---|---|---|---|
| CT | 2.33 | 0.57 | 0.55 |
| FL | 4.08 | 0.00 | 0.00 |
| FT | 1.00 | 1.00 | 1.00 |
| PT | 2.79 | 0.47 | 0.27 |
Kết quả RR (Relative Risk) là 2.33 điều này cho thấy xác suất để những kỹ sư làm việc theo hợp đồng (CT) nhận mức lương “thấp” cao hơn xác suất kỹ sư làm việc toàn thời gian đạt được mức lương thấp là khoảng 2.3
Kết quả RR của nhóm làm việc tự do (FL) là 4.08 cho thấy những kỹ sư trong nhóm làm việc tự do có khả năng nhận mức lương “thấp” với tỷ lệ cao gấp 4.08 lần so với nhóm làm việc toàn thời gian (FT).
Ngoài ra RR của nhóm PT là 2.79, kết quả có ý nghĩa xác suất để kỹ sư trong nhóm làm việc bán thời gian nhận mức lương “thấp” cao khoảng 2.79 lần so với nhóm làm việc toàn thời gian nhận được cùng mức lương này.
Khi phân tích mức lương trung ta nhân thấy RR ở nhóm CT là 0.57, kết quả này giúp ta nhận thấy tỷ lệ để kỹ sư đang làm việc theo hợp đồng nhận mức lương “trung” chỉ khoảng 0.57 lần so với nhóm làm việc toàn thời gian.
Kết quả RR của nhóm làm việc tự do là 0, cho biết được kỹ sư trong nhóm làm việc tự do không có khả năng nhận mức lương “trung” so với nhóm làm việc toàn thời gian.
Cuối cùng, RR đối với nhóm làm việc bán thời gian khoảng 0.47, có ý nghĩa là kỹ sư trong nhóm làm việc bán thời gian có khả năng để nhận mức lương tầm “trung” chỉ khoảng 0.47 so với người làm việc toàn thời gian.
Kỹ sư làm việc theo hợp đồng có khả năng nhận mức lương “cao” chỉ khoảng 0.55 lần so với nhóm làm việc toàn thời gian.
Trong nhóm làm việc tự do thì gần như không có khả năng nhận mức lương “cao” so với nhóm làm việc toàn thời gian vì chỉ số RR là 0.
Kết quả RR là 0.27 cho thấy các kỹ sư trong nhóm làm việc bán thời gian có khả năng để nhận mức lương “cao” so với kỹ sư làm việc toàn thời gian chỉ khoảng 0.27 lần tương ứng khoảng 27%
Trong phần này em sẽ thực hiện khoảng ước lượng cho RR với khoảng tin cậy 95%. Tuy nhiên để cho kết quả dễ dàng phân tích hơn trước nên em sẽ thực hiện trên bảng có số lượng nhỏ hơn là 2x2.
Trước tiên em tính khoảng ước lượng Relative Risk đối với mức lương thấp và trung nhưng trong nhóm FT và PT
RR1 <- matrix(c(4037, 9239,
26, 10),
nrow = 2,
byrow = TRUE,
dimnames = list(c( "FT", "PT"),
c("thap", "trung")))
epitab(RR1, method = 'riskratio', rev = 'c')
## $tab
## trung p0 thap p1 riskratio lower upper p.value
## FT 9239 0.6959174 4037 0.3040826 1.000000 NA NA NA
## PT 10 0.2777778 26 0.7222222 2.375086 1.936378 2.913188 3.031552e-07
##
## $measure
## [1] "wald"
##
## $conf.level
## [1] 0.95
##
## $pvalue
## [1] "fisher.exact"
Xác suất để có mức lương thấp đối với nhóm làm việc toàn thời gian là khoảng 72.2% và xác suất để đạt mức lương thấp đối với nhóm toàn thời gian là 30.4%, trong trường hợp này tỷ lệ Relative Risk (RR) là khoảng 2.375 điều này cho thấy khả năng để nhóm kỹ sư đang làm việc bán thời gian đạt được mức lương “thấp” cao gấp 2.375 lần đối với nhóm làm việc toàn thời gian. Khoảng ước lượng Relative Risk có kết quả là 1.936378 và 2.913188 đối với độ tin cậy 95%. Giá trị p.value khoảng 3.03e^-7 nhỏ hơn mức ý nghĩa 5% vì vậy kết quả này có ý nghĩa thống kê cho thấy có sự khác biệt ở nhóm làm việc toàn thời gian và bán thời gian.
Tiếp theo em tính khoảng ước lượng Relative Risk đối với mức lương thấp và trung nhưng trong trường hợp này sẽ phân tích trên hai nhóm FT và FL
RR1 <- matrix(c(4037, 9239,
14, 0),
nrow = 2,
byrow = TRUE,
dimnames = list(c( "FT", "FL"),
c("thap", "trung")))
epitab(RR1, method = 'riskratio', rev = 'c')
## $tab
## trung p0 thap p1 riskratio lower upper p.value
## FT 9239 0.6959174 4037 0.3040826 1.000000 NA NA NA
## FL 0 0.0000000 14 1.0000000 3.288581 3.205034 3.374305 5.884773e-08
##
## $measure
## [1] "wald"
##
## $conf.level
## [1] 0.95
##
## $pvalue
## [1] "fisher.exact"
Tỷ lệ Relative Risk được tính là 3.288581. Điều này cho thấy rằng nhóm kỹ sư làm việc tự do có tỷ lệ đạt mức lương “thấp” gấp khoảng 3.29 lần so với nhóm làm việc toàn thời gian. Khoảng ước lượng cho RR của nhóm làm việc tự do là từ 3.205034 đến 3.374305 với mức độ tin cậy 95%. Ngoài ra giá trị p-value là rất nhỏ, xấp xỉ 5.884773e-08. Điều này cho thấy có sự khác biệt giữa nhóm toàn thời gian và tự do về mức lương thấp là có ý nghĩa thống kê với mức ý nghĩa rất cao.
RR1 <- matrix(c(4037, 9239,
16, 9),
nrow = 2,
byrow = TRUE,
dimnames = list(c( "FT", "CT"),
c("thap", "trung")))
epitab(RR1, method = 'riskratio', rev = 'c')
## $tab
## trung p0 thap p1 riskratio lower upper p.value
## FT 9239 0.6959174 4037 0.3040826 1.000000 NA NA NA
## CT 9 0.3600000 16 0.6400000 2.104692 1.566823 2.827202 0.0006566158
##
## $measure
## [1] "wald"
##
## $conf.level
## [1] 0.95
##
## $pvalue
## [1] "fisher.exact"
Kết quả cho thấy xác suất đạt mức lương thấp ở nhóm hợp đồng là 64% và đối với nhóm toàn thời gian là 30.4%. Tỷ lệ Relative Risk là 2.104692 điều này cho thấy rằng nhóm kỹ sư làm việc theo hợp đồng có mức lương “thấp” gấp khoảng 2.1 lần so với nhóm làm việc toàn thời gian. Khoảng ước lượng cho RR của nhóm CT là từ 1.566823 đến 2.827202 với mức độ tin cậy 95%. Giá trị p-value là 0.00066 , rất nhỏ hơn mức ý nghĩa thông thường là 0.05 kết quả cho thấy sự khác biệt giữa nhóm FT và CT về mức lương thấp là có ý nghĩa thống kê.
Odd ratio và khoảng ước lượng
Chúng ta có thể hiểu rằng Odd là một đại lượng thống kê đo lường xác suất của một sự kiện xảy ra so với không xảy ra. Odd được tính bằng cách chia tỷ lệ giữa xác suất xảy ra và xác suất không xảy ra của một sự kiện. Ngoài ra ta có thêm chỉ số odd ratio là tỷ lệ giữa odds của hai nhóm khác nhau. Nó cho biết mức độ thay đổi của odd giữa hai nhóm so với nhau.
Trong phần này em sẽ tính tỷ lệ Odd Ratio của hai biến lương và công việc. Để dễ dàng trong phân tích em sẽ tách số liệu từ bảng tần số ban đầu thành bảng tần số mới chỉ bao gồm nhóm lương thấp, cao và nhóm công việc toàn thời gian và bán thời gian.
odd <- table(d$type,d$sal)
odd2 <- odd[c("FT", "PT"), c("thap", "cao")]
epitab(odd2, method = 'oddsratio')
## $tab
##
## thap p0 cao p1 oddsratio lower upper
## FT 4037 0.993600788 3178 0.9993710692 1.00000000 NA NA
## PT 26 0.006399212 2 0.0006289308 0.09771506 0.02317524 0.4120015
##
## p.value
## FT NA
## PT 2.892004e-05
##
## $measure
## [1] "wald"
##
## $conf.level
## [1] 0.95
##
## $pvalue
## [1] "fisher.exact"
Kết quả Odd Ratio của PT so với FT khoảng 0.098 điều này có nghĩa là odd mức lương cao của nhóm kỹ sư làm việc bán thời gian chỉ bằng khoảng 9.77% so với odd kỹ sư làm việc toàn thời gian.
Khoảng tin cậy của OR đối với độ tin cậy 95% là từ 0.0232 đến 0.4120 và giá trị p là 2.89e-05, nhỏ hơn 0.05 cho thấy có sự khác biệt giữa odd ở hai nhóm PT và FT
Trái ngược lại với phần trên chúng ta có thể sử dụng tham số rev = ‘c’ thay đổi giữa mức lương cao và mức lương thấp để tính odd ratio của nhóm làm việc bán thời gian so với nhóm làm việc toàn thời gian ở mức lương thấp.
epitab(odd2, method = 'oddsratio', rev = 'c' )
## $tab
##
## cao p0 thap p1 oddsratio lower upper p.value
## FT 3178 0.9993710692 4037 0.993600788 1.00000 NA NA NA
## PT 2 0.0006289308 26 0.006399212 10.23384 2.427176 43.1495 2.892004e-05
##
## $measure
## [1] "wald"
##
## $conf.level
## [1] 0.95
##
## $pvalue
## [1] "fisher.exact"
Kết quả cho thấy Odd ratio của PT so với FT là 10.23 ta có thể giải thích rằng odd của nhóm kỹ sư bán thời có mức lương “thấp” cao hơn gấp 10.23 lần so với odd của nhóm làm việc toàn thời gian có lương thấp.
Kết quả chỉ ra rằng giá trị p là 2.892004e-05, nhỏ hơn 0.05, do đó có sự khác biệt giữa odd của hai nhóm trên.
Hơn thế, khoảng ước lượng cho odd ratio đối với độ tin cậy 95% là 10.23 và 2.43
Phân tích mức lương trung giữa nhóm kỹ sư làm việc toàn thời gian và bán thời gian
odd3 <- odd[c("FT", "PT"), c("thap", "trung")]
epitab(odd3, method = 'oddsratio')
## $tab
##
## thap p0 trung p1 oddsratio lower upper
## FT 4037 0.993600788 9239 0.998918802 1.0000000 NA NA
## PT 26 0.006399212 10 0.001081198 0.1680585 0.08096884 0.3488213
##
## p.value
## FT NA
## PT 3.031552e-07
##
## $measure
## [1] "wald"
##
## $conf.level
## [1] 0.95
##
## $pvalue
## [1] "fisher.exact"
Kết quả Odds ratio là 0.168, điều này có nghĩa là odds của kỹ sư có mức lương “trung” trong nhóm làm việc bán thời gian chỉ bằng khoảng 16.8% so với odds nhóm làm việc toàn thời gian có cùng mức lương này.
Ngoài ra Khoảng tin cậy của OR là từ 0.081 đến 0.349 đối với khoảng tin cậy 95% và giá trị p là 3.03e-07, nhỏ hơn 0.05, cho thấy sự khác biệt giữa odds của nhóm PT và FT.
epitab(odd3, method = 'oddsratio', rev = 'c')
## $tab
##
## trung p0 thap p1 oddsratio lower upper
## FT 9239 0.998918802 4037 0.993600788 1.00000 NA NA
## PT 10 0.001081198 26 0.006399212 5.95031 2.866798 12.35043
##
## p.value
## FT NA
## PT 3.031552e-07
##
## $measure
## [1] "wald"
##
## $conf.level
## [1] 0.95
##
## $pvalue
## [1] "fisher.exact"
Kết quả cho thấy odds kỹ sư có mức lương thấp đang làm việc bán thời gian gấp khoảng 5.95 lần so với odds của kỹ sư đang làm việc toàn thời gian.
Hơn thế kết quả cũng cho rằng khoảng tin cậy của OR là từ 2.87 đến 12.35 và giá trị p là 3.031552e-07, nhỏ hơn 0.05, cho thấy sự khác biệt giữa odds của hai nhóm là có ý nghĩa thống kê.
Kết quả bên dưới được sử dụng để kiểm tra tính độc lập giữa hai biến loại hình công việc và mức lương bằng phương pháp chi bình phương. Giả thuyết Ho được đưa ra là biến mức lương và loại công việc độc lập.
chisq.test(table(d$type,d$sal))
##
## Pearson's Chi-squared test
##
## data: table(d$type, d$sal)
## X-squared = 97.86, df = 6, p-value < 2.2e-16
Qua kết quả kiểm định cho ta thấy giá trị p−value là 2.2e−16 nhỏ hơn 0.05 , nên bác bỏ H0, nghĩa là mức lương và loại công việc là có liên quan với nhau.
Ước lượng tỷ lệ cho nhóm lương trung bình
Trước tiên ta cần lọc kết quả về mức lương theo nhóm làm việc toàn thời gian. Kết quả thu được như sau.
FT <- table(d[d$type == "FT", ]$sal)
FT
##
## thap trung cao
## 4037 9239 3178
prop.test(FT["trung"], sum(FT), p = 0.7)
##
## 1-sample proportions test with continuity correction
##
## data: FT["trung"] out of sum(FT), null probability 0.7
## X-squared = 1502.2, df = 1, p-value < 2.2e-16
## alternative hypothesis: true p is not equal to 0.7
## 95 percent confidence interval:
## 0.5538791 0.5691017
## sample estimates:
## p
## 0.5615048
Sau đó ước lượng tỷ lệ với giả thuyết không (H0): Tỷ lệ người có mức lương trung trong nhóm toàn thời gian là 0.7 (70%).
Kết quả từ kiểm định cho thấy giá trị p-value là 2.2e-16, rất gần với 0 và nhỏ hơn giá trị α=0.05. Điều này cho thấy có bằng chứng mạnh mẽ để bác bỏ giả thuyết rằng tỷ lệ người có mức lương trung trong nhóm toàn thời gian là 70%.
Kết quả cũng cho thấy tỷ lệ kỹ sư có mức lương trung trong nhóm toàn thời gian khoảng 56.15%. Khoảng ước lượng tỷ lệ ở độ tin cây 95% nằm trong khoảng từ 55.39% đến 56.91%.
Ước lượng tỷ lệ cho nhóm lương thấp
prop.test(FT["thap"], sum(FT), p = 0.2)
##
## 1-sample proportions test with continuity correction
##
## data: FT["thap"] out of sum(FT), null probability 0.2
## X-squared = 211.22, df = 1, p-value < 2.2e-16
## alternative hypothesis: true p is not equal to 0.2
## 95 percent confidence interval:
## 0.2388058 0.2520150
## sample estimates:
## p
## 0.2453507
Ta có giả thuyết không (H0): Tỷ lệ người có mức lương thấp trong nhóm toàn thời gian là 0.2 (20%). Vì giá trị p rất nhỏ (p < 2.2e-16), chúng ta bác bỏ giả thuyết không (H0). Điều này có nghĩa là tỷ lệ kỹ sư có mức lương “thấp” trong nhóm làm việc toàn thời gian không bằng 20%. Kết quả cho rằng tỷ lệ này đạt khoảng 24.5% và ở mức độ tin cậy 95% kết quả ước lượng chỉ ra tỷ lệ kỹ sư có mức lương thấp trong nhóm toàn thời gian đạt khoảng từ 23.88% đến 25.20%.
Ước lượng tỷ lệ cho nhóm lương cao
prop.test(FT["cao"], sum(FT), p = 0.3)
##
## 1-sample proportions test with continuity correction
##
## data: FT["cao"] out of sum(FT), null probability 0.3
## X-squared = 894.13, df = 1, p-value < 2.2e-16
## alternative hypothesis: true p is not equal to 0.3
## 95 percent confidence interval:
## 0.1871545 0.1992785
## sample estimates:
## p
## 0.1931445
Ta có giả thuyết không (H0) là tỷ lệ kỹ sư có mức lương “cao” trong nhóm làm toàn thời gian là 30%. Kết quả bên trên có giá trị p rất nhỏ (p < 2.2e-16), chúng ta bác bỏ giả thuyết không (H0). Điều này có nghĩa là tỷ lệ kỹ sư có mức lương “cao” trong nhóm nhóm làm toàn thời gian sẽ không bằng 30%. Thực tế, tỷ lệ kỹ sư có mức lương “cao” chỉ khoảng 19.31%, thấp hơn so với tỷ lệ kỳ vọng là 30%. Ngoài ra, khoảng tin cậy 95% cho tỷ lệ kỹ sư có mức lương cao trong nhóm toàn thời gian là [0.1872, 0.1993]. Điều này có nghĩa là chúng ta có 95% tin tưởng rằng tỷ lệ này nằm trong khoảng từ 18.72% đến 19.93%
Mô hình xác suất tuyến tính
Phần này sẽ thực hiện hồi quy xác suất tuyến tính giữa biến độc lập là loại hình công việc và biến phụ thuộc là lương. Tuy nhiên trong phần này ta cần thay đổi biến sal thành biến nhị phân với giá trị 1 tương ứng với mức lương trung và cao, ngược lại giá trị 0 tương ứng với mức lương thấp.
d$sal2 <- ifelse(d$sal %in% c("trung", "cao"), 1, 0)
d$type <- relevel(as.factor(d$type), ref = "FT")
glm <- glm(sal2 ~ type, data = d)
summary(glm)
##
## Call:
## glm(formula = sal2 ~ type, data = d)
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.754649 0.003355 224.928 < 2e-16 ***
## typeCT -0.326078 0.081400 -4.006 6.21e-05 ***
## typeFL -0.754649 0.115069 -6.558 5.61e-11 ***
## typePT -0.438860 0.069895 -6.279 3.50e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 0.185214)
##
## Null deviance: 3079.8 on 16533 degrees of freedom
## Residual deviance: 3061.6 on 16530 degrees of freedom
## AIC: 19047
##
## Number of Fisher Scoring iterations: 2
\(\hat\pi= 0.754649-0.326078CT-0.754649FL-0.438860PT\)
Tất cả các hệ số đều có giá trị p nhỏ hơn mức ý nghĩa \(\alpha= 5%\) , ngụ ý rằng tất cả các hệ số đều có ảnh hưởng ý nghĩa thống kê đến mức lương.
Hệ số \(\hat\beta_0 = 0.754649\) có ý nghĩa rằng tỷ lệ để người kỹ sư đạt được mức lương trung bình là 75.46% trong trường hợp không xuất hiện các nhóm công việc khác.
Đối với hệ số \(\hat\beta_1 =-0.326078\) ta nhận thấy khi kỹ sư làm việc theo hình thức hợp đồng và không bị ảnh hưởng bởi các yếu tố khác thì chênh lệch tỷ lệ đạt mức lương trung trở lên giữa nhóm làm việc hợp đồng và nhóm toàn thời gian là - 32.6%, hay có thể nói cách khác tỷ lệ bình quân để đạt mức lương trung trở lên trong nhóm này là 42.86%
Tương tự \(\hat\beta_2 =-0.754649\) ta nhận thấy khi kỹ sư làm việc theo hình thức tự do và không bị ảnh hưởng bởi các yếu tố khác thì chênh lệch tỷ lệ đạt mức lương trung trở lên giữa nhóm làm việc tự do và nhóm toàn thời gian là - 75.46%, hay tỷ lệ đạt mức lương trung trở lên trong nhóm này là 0%
Cuối cùng với \(\hat\beta_3 =-0.4388\) ta nhận thấy khi kỹ sư làm việc theo hình thức bán thời gian và không bị ảnh hưởng bởi các yếu tố khác thì chênh lệch tỷ lệ đạt mức lương trung trở lên giữa nhóm làm việc bán thời gian và nhóm toàn thời gian là - 43.88%, hay tỷ lệ đạt mức lương trung trở lên trong nhóm này là 31.58%
Ta nhận thấy ở cả ba nhóm hợp đồng, bán thời gian, tự do đều có mối quan hệ âm với mức lương lương từ tầm trung trở lên. Nhưng nếu không có sự tác động của ba nhóm trên thì xác suất để đạt mức lương trên là khá cao khoảng 0.754649
n <- data.frame(type = c("CT", "FL", "PT","FT") )
predict.glm(glm, newdata = n, type="response")
## 1 2 3 4
## 4.285714e-01 2.013945e-13 3.157895e-01 7.546493e-01
Dựa vào kết quả ta có thể nhận thấy xác suất dự báo cho biểu hiện kỹ sư làm công việc hợp đồng (1) là khoảng 42.86%. Điều này có nghĩa một kỹ sư làm việc theo hợp đồng có xác suất nhận mức lương trung đến cao là khoảng 42.86% .
Bên cạnh đó kết quả dự báo đối với nhóm kỹ sư làm việc tự do (2) thì xác suất để đạt được mức lương trung trở lên chỉ khoảng 3.06e-1, vì rất nhỏ cho nên đối với nhóm này gần như không có xác suất để nhận được mức lương trung trở lên.
Trong nhóm kỹ sư làm việc thời bán thời gian (3) thì xác suất để đạt mức lương trung chỉ khoảng 31.58%, cao hơn so với xác suất của nhóm làm tự do.
Cuối cùng kỹ sư làm việc toàn thời gian có xác suất đạt mức lương từ trung trở lên là 75.46%, cho thấy rằng làm việc toàn thời gian mang lại cơ hội lớn hơn để đạt được mức lương cao hơn.
Mô hình logit
Ngoài ra chúng ta có thể thực hiện mô hình hồi quy khi có hàm liên kết dạng \(g(μ)=logit(μ)=ln(\frac{μ}{1−μ})\) khi phân tích mối quan hệ giữa nhóm việc làm với mức lương, trong đó nhóm công việc toàn thời gian sẽ làm nhóm cơ sở.
logit <- glm(factor(sal2) ~ type, data = d, family = binomial(link = 'logit'))
summary(logit)
##
## Call:
## glm(formula = factor(sal2) ~ type, family = binomial(link = "logit"),
## data = d)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.12356 0.01812 62.015 < 2e-16 ***
## typeCT -1.41125 0.38231 -3.691 0.000223 ***
## typeFL -13.68963 86.79141 -0.158 0.874669
## typePT -1.89675 0.34946 -5.428 5.71e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 18506 on 16533 degrees of freedom
## Residual deviance: 18421 on 16530 degrees of freedom
## AIC: 18429
##
## Number of Fisher Scoring iterations: 11
Ta có kết từ hàm hồi quy thông qua mô hình logit như sau:
\(ln(\frac{\hat\pi}{1−\hat\pi})=1.12356−1.41125CT−13.68963FL−1.89675PT\)
Hệ số chặn là 1.12356 điều này có nghĩa là khi không có sự tác động của nhóm công việc khác thì \(ln(\frac{\hat\pi}{1−\hat\pi})\) của kỹ sư toàn thời gian đạt mức lương trung là 1.12356
Đối với hệ số \(\hat\beta_1 =−1.41125\) điều này cho thấy chênh lệch \(ln(\frac{\hat\pi}{1−\hat\pi})\) đạt mức lương trung trở lên giữa nhóm làm việc hợp đồng và nhóm toàn thời gian là −1.41125
Đối với hệ số \(\hat\beta_2 =−13.68963\) điều này cho thấy chênh lệch \(ln(\frac{\hat\pi}{1−\hat\pi})\) đạt mức lương trung trở lên giữa nhóm làm việc tự do và nhóm toàn thời gian là −13.68963
Tương tự, hệ số \(\hat\beta_3 =−1.89675\) điều này cho thấy chênh lệch \(ln(\frac{\hat\pi}{1−\hat\pi})\) mức lương trung trở lên giữa nhóm làm việc bán thời gian và nhóm toàn thời gian là −1.89675
Sau đó từ mô hình logit ta dự báo xác suất của biến lương tương ứng với mỗi nhóm việc làm. Ta cần một khung dữ liệu n với một cột type chứa các giá trị “CT”, “FL”, “PT”, và “FT”.
n <- data.frame(type = c("CT", "FL", "PT","FT") )
predict.glm(logit, newdata = n, type="response")
## 1 2 3 4
## 4.285714e-01 3.488403e-06 3.157895e-01 7.546493e-01
Dựa vào kết quả ta có thể nhận thấy xác suất dự báo cho biểu hiện CT là khoảng 42.86%. Điều này có nghĩa một kỹ sư làm việc theo hợp đồng có xác suất khoảng 42.86% có mức lương trung đến cao.
Bên cạnh đó kết quả dự báo đối với nhóm kỹ sư làm việc tự do thì xác suất để đạt được mức lương trung trở lên chỉ khoảng 0.00035%
Trong nhóm kỹ sư làm việc thời bán thời gian thì xác suất để đạt mức lương trung chỉ khoảng 31.58%, cao hơn so với xác suất của nhóm làm tự do.
Cuối cùng nhóm có xác suất cao nhất để đạt mức lương trung là ở công việc toàn thời gian với xác suất khoảng 75.46%.
Mô hình Probit
Thay vì sử dụng mô hình logit ta có thể thay thế bằng mô hình probit dự đoán xác suất của một sự kiện xảy ra, mô hình probit sử dụng hàm phân phối chuẩn. Trong đó hàm liên kết có dạng \(probit(π) = Φ^{-1}(π)\) và hàm xác suất là \(π(x)=Φ(β_0+β_1X+β_2X_2+...+β_nX_n)\). Ta thực hiện thao tác khá tương tự với phân trên
probit <- glm(factor(sal2) ~ type, data = d, family = binomial(link = 'probit'))
summary(probit)
##
## Call:
## glm(formula = factor(sal2) ~ type, family = binomial(link = "probit"),
## data = d)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.68919 0.01066 64.637 < 2e-16 ***
## typeCT -0.86921 0.23849 -3.645 0.000268 ***
## typeFL -5.30055 24.62168 -0.215 0.829549
## typePT -1.16870 0.21231 -5.505 3.7e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 18506 on 16533 degrees of freedom
## Residual deviance: 18421 on 16530 degrees of freedom
## AIC: 18429
##
## Number of Fisher Scoring iterations: 11
Ta có kết quả mô hình hồi quy như sau:
\(\hat\pi= \Phi(0.68919−0.86921×CT−5.30055×FL−1.16870×PT)\)
Hệ số chặn là 0.68919 và có ý nghĩa thống kê khi giá trị p_value nhỏ hơn mức ý nghĩa 5%. Điều này có nghĩa là khi không có sự tác động các biến độc lập khác thì xác suất để kỹ sư đạt mức lương trung là \(\Phi(0.68919)\)
Hệ số \(\hat\beta_1= −0.86921\) có ý nghĩa là trong trường hợp các yếu tố khác không đổi thì chênh lệch tỷ lệ đạt mức lương trung giữa nhóm làm việc theo hình thức hợp đồng và nhóm làm toàn thời gian đạt khoảng \(\Phi(-0.86921)\)
Hệ số \(\hat\beta_2=−5.30055\) có ý nghĩa là trong trường hợp các yếu tố khác không đổi thì chênh lệch tỷ lệ đạt mức lương trung giữa nhóm làm việc theo hình thức tự do và nhóm làm toàn thời gian đạt khoảng \(\Phi(−5.30055)\)
Hệ số \(\hat\beta_3= −1.16870\) có ý nghĩa là trong trường hợp các yếu tố khác không đổi thì chênh lệch tỷ lệ đạt mức lương trung giữa nhóm làm việc theo hình thức bán thời gian và nhóm làm toàn thời gian đạt khoảng \(\Phi(−1.16870)\)
Từ kết quả trên em sẽ thực hiện dự báo cho tỷ lệ để đạt được mức lương trung trong các nhóm công việc. Kết quả được trình bày trong bảng bên dưới.
predict.glm(probit, newdata = n, type="response")
## 1 2 3 4
## 4.285714e-01 2.000259e-06 3.157895e-01 7.546493e-01
Dựa vào kết quả ta có thể nhận thấy xác suất dự báo cho biểu hiện kỹ sư làm công việc hợp đồng (1) là khoảng 42.86%. Điều này có nghĩa một kỹ sư làm việc theo hợp đồng có xác suất khoảng 42.86% có mức lương trung đến cao.
Bên cạnh đó kết quả dự báo đối với nhóm kỹ sư làm việc tự do (2) thì xác suất để đạt được mức lương trung trở lên chỉ khoảng 0.000002%
Trong nhóm kỹ sư làm việc thời bán thời gian (3) thì xác suất để đạt mức lương trung chỉ khoảng 31.58%, cao hơn so với xác suất của nhóm làm tự do.
Cuối cùng nhóm có xác suất cao nhất để đạt mức lương trung là ở công việc toàn thời gian với xác suất khoảng 75.46% đây là lựa chọn tốt nhất nếu mục tiêu là đạt được mức lương từ trung trở lên.
Chỉ số AIC
AIC (Akaike Information Criterion) là một chỉ số được sử dụng để đánh giá và so sánh các mô hình thống kê. Mô hình có AIC thấp hơn được coi là tốt hơn vì nó cân bằng giữa độ phù hợp và độ phức tạp của mô hình. Do đó chúng ta sử dụng chỉ số này để so sánh ba mô hình bên trên để lựa chọn mô hình phù hợp nhất trong phân tích xác suất để đạt được mức lương trung.
AIC(glm )
## [1] 19047.11
AIC(logit )
## [1] 18428.97
AIC(probit )
## [1] 18428.97
Mô hình Logit và Probit đều tốt hơn mô hình xác suất tuyến tính trong việc giải thích mối quan hệ giữa loại công việc và xác suất đạt mức lương từ trung trở lên, Vì chỉ số AIC của hai mô hình này thấp hơn chỉ số AIC của mô hình xác suất tuyến tính. Logit và Probit có hiệu suất tương đương nhau vì chúng có cùng chỉ số AIC
Hệ số Brier Score
Phần này sẽ tính chỉ số Brier Score được sử dụng để đo lường độ chính xác của các dự đoán xác suất. Nó tính toán độ lệch bình phương giữa dự đoán xác suất và giá trị thực tế. Một giá trị Brier Score thấp cho thấy mô hình dự đoán tốt hơn. Vì thế ta sẽ tính và so sánh giá trì này ở hai mô hình hồi quy logit và probit.
BrierScore(logit)
## [1] 0.1851692
BrierScore(probit)
## [1] 0.1851692
Kết quả bên trên cho thấy Brier Score cho logit và probit đều là 0.1851692. Điều này có nghĩa là cả hai mô hình đều có hiệu suất tương đương khi dự đoán xác suất cho biến mục tiêu. Brier Score khá thấp (gần 0), cho thấy rằng cả hai mô hình đều dự đoán khá tốt và không sai lệch nhiều so với giá trị thực tế.
Ma trận nhầm lẫn
Conf(table(predict(probit, type="response") >=0.5,probit$data$sal2 == '1'))
##
## Confusion Matrix and Statistics
##
## Reference
## Prediction FALSE TRUE
## FALSE 56 24
## TRUE 4037 12417
##
## Total n : 16'534
## Accuracy : 0.7544
## 95% CI : (0.7478, 0.7609)
## No Information Rate : 0.7524
## P-Value [Acc > NIR] : 0.2855
##
## Kappa : 0.0175
## Mcnemar's Test P-Value : < 2.2e-16
##
## Sensitivity : 0.0137
## Specificity : 0.9981
## Pos Pred Value : 0.7000
## Neg Pred Value : 0.7546
## Prevalence : 0.2476
## Detection Rate : 0.0048
## Detection Prevalence : 0.0034
## Balanced Accuracy : 0.5059
## F-val Accuracy : 0.0268
## Matthews Cor.-Coef : 0.0731
##
## 'Positive' Class : FALSE
Từ kết quả trên ma trận cho ta biết có 56 số lần mô hình dự đoán đúng là mức lương thấp. Tuy nhiên có 24 lần mô hình dự đoán sai đối với mức lương trung và khoảng 4037 lần mô hình dự đoán sai mức lương thấp. Số lần mô hình dự đoán đúng ở mức lương trung là 12417 lần đạt cao nhất.
Ngoài ra ta có tỷ lệ dự đoán đúng của mô hình là 0.7544 và ở độ tin cậy là 95% thì tỷ lệ này nằm trong khoảng từ (0.7478, 0.7609)
d$sal <- factor(d$sal, levels = c("cao", "trung", "thap")) #thap làm tham chiếu
levels(d$sal)
## [1] "cao" "trung" "thap"
library(VGAM)
mlog <- vglm(sal ~ type, family = multinomial,data =d)
mlog
##
## Call:
## vglm(formula = sal ~ type, family = multinomial, data = d)
##
##
## Coefficients:
## (Intercept):1 (Intercept):2 typeCT:1 typeCT:2 typeFL:1
## -0.2392498 0.8279318 -1.4347267 -1.4032960 -14.9120587
## typeFL:2 typePT:1 typePT:2
## -15.9792403 -2.3256996 -1.7834433
##
## Degrees of Freedom: 33068 Total; 33060 Residual
## Residual deviance: 32570.03
## Log-likelihood: -16285.02
##
## This is a multinomial logit model with 3 levels
Đối với mức lương cao
\(\ln\left(\frac{\pi_1}{\pi_3}\right) = -0.2392498-1.4347267*CT-14.9120587*FL-2.3256996*PT\)
Hệ số -0.2392 có ý nghĩa là log giữa xác suất đạt được mức lương “cao” và xác suất đạt được mức lương “thap” đối với nhóm công việc toàn thời gian là -0.2392 trong trường hợp không có các yếu khác. Suy ra tỷ lệ lương cao trên mức lương thấp đối với nhóm làm toàn thời gian là e^(-0.2392) = 0.7873. Dựa vào công thức \(\pi_i = \frac{e^{\beta_{0i} + \beta_{1i}X}}{1+\sum_{j=1}^{K-1} e^{\beta_{0j}+\beta_{1j}X}}\) ta cũng có thể tính xác suất để đạt mức lương cao trong nhóm làm việc toàn thời gian là \(\pi_1 = \frac{e^{-0.2392498}}{1 + e^{-0.2392498} + e^{0.8279318}}=0.1931445\)
Chêch lệch của \(\ln\left(\frac{\pi_1}{\pi_3}\right)\) giữa công việc hợp đồng và toàn thời gian (tham chiếu) là −1.434, điều này cho thấy \(\ln\left(\frac{\pi_1}{\pi_3}\right)\) đối với công việc hợp đồng là -0.2392498-1.4347267= -1.673976. Suy ra tỷ lệ lương cao trên mức lương thấp đối với nhóm làm việc hợp đồng là e^(-1.673976) = 0.1875. Xác suất để đạt mức lương cao trong nhóm hợp đồng (CT = 1, FL = 0, PT = 0) là \(\pi_1 = \frac{e^{-0.2395-1.4347}}{1 + e^{-0.2392−1.4347} + e^{0.8279-1.4032}}=0.1071\)
Hệ số -14.912 có ý nghĩa là chêch lệch của \(\ln\left(\frac{\pi_1}{\pi_3}\right)\) giữa công việc tự do(FT) và toàn thời gian (tham chiếu) là -14.912, điều này cho thấy \(\ln\left(\frac{\pi_1}{\pi_3}\right)\) đối với công việc hợp đồng là -0.2392498-14.9120587= -15.1513. Suy ra tỷ lệ lương cao trên mức lương thấp đối với nhóm làm việc tự do là e^(-15.1513) = 2.629505e-07. Xác suất để đạt mức lương cao trong nhóm này (CT = 0, FL = 1, PT = 0) là \(\pi_1 = \frac{e^{-0.2395-14.9120}}{1 + e^{-0.2392-14.9120} + e^{0.8279-15.9792}}=0\)
Hệ số -2.3256 có ý nghĩa là chêch lệch của \(\ln\left(\frac{\pi_1}{\pi_3}\right)\) giữa công việc bán thời gian và toàn thời gian (tham chiếu) là -2.3256, điều này cho thấy \(\ln\left(\frac{\pi_1}{\pi_3}\right)\) đối với công việc hợp đồng là -0.2392498-2.3256= -2.56485. Suy ra tỷ lệ lương cao trên mức lương thấp đối với nhóm bán thời gian là e^(-2.56485) = 0.07693. Vậy Xác suất để đạt mức lương cao trong nhóm bán thời gian là (CT = 0, FL = 0, PT = 1) là \(\pi_1 = \frac{e^{-0.2395-2.3256}}{1 + e^{-0.2392-2.3256} + e^{0.8279-1.783}}=0.0526\)
Đối với mức lương trung
\(\ln\left(\frac{\pi_2}{\pi_3}\right) = 0.8279318-1.4032960*CT-15.9792403 *FL-1.7834433*PT\)
Hệ số 0.8279 có ý nghĩa là log giữa xác suất đạt được mức lương “trung” và xác suất đạt được mức lương “thap” đối với nhóm công việc toàn thời gian là 0.8279 trong trường hợp không có các yếu khác. Tương tự như trên, ta cũng có thể tính xác suất để đạt mức lương “trung” trong nhóm làm việc toàn thời gian (CT = 0, FL = 0, PT = 0) là \(\pi_2 = \frac{e^{0.8279}}{1 + e^{-0.2392} + e^{0.8279}}=0.5615\)
Chêch lệch của \(\ln\left(\frac{\pi_2}{\pi_3}\right)\) giữa hai nhóm công việc hợp đồng và toàn thời gian là -1.4033. Xác suất để đạt mức lương “trung” trong nhóm hợp đồng (CT = 1, FL = 0, PT = 0) là \(\pi_2 = \frac{e^{0.8279-1.4032}}{1 + e^{-0.2392−1.4347} + e^{0.8279-1.4032}}=0.3214\)
Hệ số 15.9792 có ý nghĩa là chêch lệch của \(\ln\left(\frac{\pi_2}{\pi_3}\right)\) giữa công việc tự do(FT) và toàn thời gian là 15.9792. Ta có xác suất để đạt mức lương “trung” trong nhóm tự do (CT = 0, FL = 1, PT = 0) là \(\pi_2 = \frac{e^{0.8279-15.9792}}{1 + e^{-0.2392-14.9120} + e^{0.8279-15.9792}}=0\)
Chêch lệch của \(\ln\left(\frac{\pi_2}{\pi_3}\right)\) giữa công việc bán thời gian và toàn thời gian là -1.7834. Hơn thế, Xác suất để đạt mức lương “trung” trong nhóm bán thời gian là (CT = 0, FL = 0, PT = 1) là \(\pi_2 = \frac{e^{0.8279-1.783}}{1 + e^{-0.2392-2.3256} + e^{0.8279-1.783}}=0.2632\)
Đối với mức lương thấp
Từ đó, Tỷ lệ để đạt được mức lương thấp trong nhóm làm việc toàn thời gian là 1-0.5615-0.19314= 24.5%
Tỷ lệ kỹ sư đạt được mức lương thấp trong nhóm làm việc hợp đồng là 1-0.3214-0.107= 57.15%
Ngoài ra tỷ lệ kỹ sư đạt được mức lương thấp trong nhóm làm việc tự do là 100%
Cuối cùng tỷ lệ kỹ sư đạt được mức lương thấp trong nhóm làm việc bán thời gian là 1 - 0.2632 - 0.0526 = 68.42%
Trong phần này em lập bảng tần số cho biến remote_ratio nhằm thống kê số lượng hình thức làm việc và so sánh sự khác biệt về số lượng cũng như tỷ lệ giữa các hình thức làm việc nhằm cung cấp trực quan ý nghĩa của biến remote_ratio trước khi phân tích tác động của hình thức làm việc đến mức lương của kỹ sư.
Trước tiên chúng ta cần đổi tên biến từ “remote_ratio” thành “remote” điều này sẽ giúp câu lệnh được thực hiện một cách ngắn gọn và dễ nhận biết hơn.
d <- d %>% rename(remote = remote_ratio) # đổi tên biến
Tiếp đến, em lập bảng tần số cho biến remote, kết quả như sau:
table(d$remote)
##
## 0 50 100
## 11118 249 5167
# Tạo bảng tần số cho biến type
f <- as.data.frame(table(d$remote))
colnames(f) <- c("Remote", "Count")
ggplot(f, aes(x = Remote, y = Count, fill = Remote)) +
geom_bar(stat = "identity", width = 0.5) +
geom_text(aes(label = Count), vjust = -0.5, size = 3.5) + # Hiển thị số liệu bên trên cột
labs(title = "Đồ thị cột hình thức làm việc ", x = "Hình thức làm việc", y = "Frequency") +
theme_minimal()
Kết quả bảng tần số và đồ thị cho thấy có 11118 kỹ sư không làm việc theo hình thức từ xa (0), 249 kỹ sư làm việc theo hình thức vừa từ xa vừa trực tiếp (50) và 5167 kỹ sư làm việc theo hình thức từ xa (100). Nhìn chung, số lượng kỹ sư không làm việc theo hình thức từ xa (0) chiếm số lượng cao nhất với 11118 kỹ sư. Trong khi đó, số lượng kỹ sư vừa làm việc trực tiếp và vừa làm việc từ xa (50) có số lượng thấp nhất là 249 kỹ sư.
Ở khía cạnh khác, chúng ta có thể phân tích tỷ lệ phần trăm đối với từng nhóm hình thức làm việc và trực quan chúng bằng đồ thị tròn để có thể dễ dàng so sánh với nhau.
prop.table(table(d$remote))*100
##
## 0 50 100
## 67.243256 1.505988 31.250756
f$Percentage <- (f$Count / sum(f$Count)) * 100
ggplot(f, aes(x = "", y = Percentage, fill = Remote)) +
geom_bar(width = 1, stat = "identity") +
coord_polar("y") +
labs(title = "Tỷ lệ phần trăm", x = "", y = "") +
theme_void() +
theme(legend.title = element_blank()) +
geom_text(aes(label = paste(round(Percentage, 2), "%")),
position = position_stack(vjust = 0.5), size = 4)
Dựa vào kết quả bảng tần suất và đồ thị, tỷ lệ kỹ sư không làm việc theo
hình thức từ xa (0) chiếm tỷ lệ cao nhất 67,24%. Trong khi số lượng kỹ
sư vừa làm việc theo hình thức từ xa và làm việc trực tiếp (50) chiếm tỷ
lệ thấp nhất 1,51%. Tỷ lệ kỹ sư lựa chọn hình thức làm việc từ xa chiếm
tỷ lệ tương đối cao là 31,25%
Chúng ta lập bảng tần số giữa hai biến remote (hình thức làm việc) và sal (mức lương phân tổ thành ba nhóm). Bảng tần số giúp chúng ta hiểu rõ hơn về sự phân bố mức lương theo từng hình thức làm việc.
table(d$remote, d$sal) %>% addmargins()
##
## cao trung thap Sum
## 0 2361 6146 2611 11118
## 50 11 54 184 249
## 100 811 3058 1298 5167
## Sum 3183 9258 4093 16534
Sau đó có thể thực hiện vẽ biểu đồ cột thể hiện số lượng kỹ sư ở mỗi nhóm lương theo hình thức làm việc khác nhau để trực quan hóa số liệu từ bảng số liệu phía trên. Điều này giúp ta dễ dàng nhận xét và so sánh.
table_data1 <- table(d$remote, d$sal)
# Chuyển bảng chéo sang định dạng dài
long_data1 <- as.data.frame(table_data1)
colnames(long_data1) <- c("type", "salary", "count")
# Vẽ biểu đồ cột
ggplot(long_data1, aes(x = salary, y = count, fill = salary)) +
geom_bar(stat = "identity", position = position_dodge()) +
geom_text(aes(label = count), vjust = 1, size = 2.5, position = position_dodge(1)) +
facet_wrap(~ type, scales = "free_y", nrow = 2) +
labs(title = "Phân bố mức lương theo hình thức công việc",
x = "Mức lương",
y = "Số lượng nhân viên",
fill = "Mức lương") +
theme_minimal()
Dựa vào đồ thị chúng ta có thể thấy rằng:
Số lượng kỹ sư làm không làm việc từ xa (0) có số lượng cao nhất là 11118 kỹ sư với tỷ lệ 67,24%. Trong đó có 2611 kỹ sư có mức lương thấp, 6146 kỹ sư có mức lương trung và 2361 kỹ sư có mức lương cao.
Số lượng kỹ sư vừa làm việc từ xa và trực tiếp (50) có số lượng thấp nhất là 249 kỹ sư chiếm tỷ lệ 1,51%. Ở hình thức làm việc này, có 184 kỹ sư có mức lương thấp, 54 kỹ sư có mức lương trung và 11 kỹ sư có mức lương cao.
Trong khi đó, số lượng kỹ sư làm việc từ xa (100) có số lượng 5157 kỹ sư chiếm tỷ lệ 31,25%. Trong đó có 1298 kỹ sư có mức lương thấp, 3058 kỹ sư có mức lương trung bình và 811 kỹ sư có mức lương cao.
Trong phần này em sẽ thực hiện tính Relative Risk nhằm phân tích sự ảnh hưởng của hình thức làm việc của các kỹ sư lên các mức lương khác nhau.
RR1 <- matrix(c(2611, 6146, 2361,
184, 54, 11,
1298, 3058, 811),
nrow = 3,
byrow = TRUE,
dimnames = list(c("0", "50", "100"),
c("thap", "trung", "cao")))
RR1
## thap trung cao
## 0 2611 6146 2361
## 50 184 54 11
## 100 1298 3058 811
Dựa trên bảng tần số đã được giải thích rõ ở phần trước em sẽ tính toán Relative Risk trong đó biến độc lập là hình thức và phụ thuộc là mức lương. Nhóm công việc được chọn làm tham chiếu là 0 tương ứng với kỹ sư lựa chọn hình thức không làm việc từ xa
Trước tiên cần tính tỷ lệ mức lương thấp, trung, cao trong tổng số người ở mỗi hình thức làm việc khác nhau. Sau đó lấy tỷ lệ vừa tính được chia cho tỷ lệ ở nhóm tham chiếu, trong ví dụ là nhóm kỹ sư không làm việc từ xa “0”.
RR_table1 <- as.data.frame(RR1)
# Tính tổng số người trong mỗi nhóm
RR_table1$Total <- rowSums(RR_table1)
# Tính tỷ lệ người ở mỗi mức lương trong mỗi nhóm (Tính Relative Risk cho mỗi nhóm bằng cách chia tỷ lệ của mỗi nhóm cho tỷ lệ của nhóm tham chiếu)
RR_table1$Rate_thap <- RR_table1$thap / RR_table1$Total
RR_table1$Rate_trung <- RR_table1$trung / RR_table1$Total
RR_table1$Rate_cao <- RR_table1$cao / RR_table1$Total
# Tính tỷ lệ cho nhóm tham chiếu (FT)
rate_ref_thap1 <- RR_table1["0", "Rate_thap"]
rate_ref_trung1 <- RR_table1["0", "Rate_trung"]
rate_ref_cao1 <- RR_table1["0", "Rate_cao"]
# Tính Relative Risk cho các nhóm so với nhóm tham chiếu (FT)
RR_table1$RR_thap <- RR_table1$Rate_thap / rate_ref_thap1
RR_table1$RR_trung <- RR_table1$Rate_trung / rate_ref_trung1
RR_table1$RR_cao <- RR_table1$Rate_cao / rate_ref_cao1
# In kết quả
print(RR_table1)
## thap trung cao Total Rate_thap Rate_trung Rate_cao RR_thap RR_trung
## 0 2611 6146 2361 11118 0.2348444 0.5527973 0.21235834 1.000000 1.0000000
## 50 184 54 11 249 0.7389558 0.2168675 0.04417671 3.146576 0.3923092
## 100 1298 3058 811 5167 0.2512096 0.5918328 0.15695762 1.069685 1.0706145
## RR_cao
## 0 1.0000000
## 50 0.2080291
## 100 0.7391168
Nhận xét tỷ lệ Relative Risk ở mức lương thấp trong các nhóm công việc
Nhận xét: Kỹ sư trong nhóm này có khả năng nhận mức lương thấp cao hơn 3.15 lần so với nhóm làm việc toàn thời gian (0).
Nhận xét: Kỹ sư trong nhóm này có khả năng nhận mức lương thấp cao hơn 1.07 lần so với nhóm làm việc trực tiếp (0).
Nhận xét tỷ lệ Relative Risk đối với mức lương trung
Nhận xét: Nhân viên trong nhóm vừa làm việc từ xa và vừa làm việc trực tiếp có khả năng nhận mức lương trung chỉ khoảng 0.39 lần so với nhóm làm việc trực tiếp (0).
Nhận xét: Nhân viên trong nhóm này có khả năng nhận mức lương trung cao hơn 1.07 lần so với nhóm làm việc trực tiếp (0).
Nhận xét tỷ lệ Relative Risk đối với mức lương cao
Nhận xét: Nhân viên trong nhóm này có khả năng nhận mức lương cao chỉ khoảng 0.21 lần so với nhóm làm việc trực tiếp (0).
Nhận xét: Nhân viên trong nhóm này có khả năng nhận mức lương cao chỉ khoảng 0.74 lần so với nhóm làm việc trực tiếp (0).
Tổng hợp nhận xét
Mức lương thấp: Kỹ sư làm việc kết hợp từ xa và trực tiếp (50) có khả năng nhận mức lương thấp cao hơn nhiều (3.15 lần) so với nhóm làm việc trực tiếp (0).
Nhóm làm việc từ xa (100) có khả năng nhận mức lương thấp cao hơn một chút (1.07 lần) so với nhóm làm việc trực tiếp (0).
Mức lương trung: Kỹ sư làm việc kết hợp từ xa và trực tiếp (50) có khả năng nhận mức lương trung thấp hơn (0.39 lần) so với nhóm làm việc toàn thời gian (0). Nhóm làm việc từ xa (100) có khả năng nhận mức lương trung cao hơn một chút (1.07 lần) so với nhóm làm việc trực tiếp (0).
Mức lương cao: Kỹ sư làm việc kết hợp từ xa và trực tiếp (50) có khả năng nhận mức lương cao thấp hơn nhiều (0.21 lần) so với nhóm làm việc toàn thời gian (0). Ký sư làm việc từ xa (100) có khả năng nhận mức lương cao thấp hơn (0.74 lần) so với nhóm làm việc toàn thời gian (0).
Trong phần này em sẽ tính tỷ lệ Odd Ratio của hai biến lương và hình thức làm việc. Để dễ dàng trong phân tích em sẽ tách số liệu từ bảng tần số ban đầu thành bảng tần số mới chỉ bao gồm nhóm lương thấp, cao và hình thức làm việc trực tiếp, hình thức làm việc từ xa
odd1 <- table(d$remote,d$sal)
odd4 <- odd1[c("0", "100"), c("thap", "cao")]
epitab(odd4, method = 'oddsratio')
## $tab
##
## thap p0 cao p1 oddsratio lower upper p.value
## 0 2611 0.6679458 2361 0.7443253 1.0000000 NA NA NA
## 100 1298 0.3320542 811 0.2556747 0.6909666 0.6227809 0.7666176 2.419623e-12
##
## $measure
## [1] "wald"
##
## $conf.level
## [1] 0.95
##
## $pvalue
## [1] "fisher.exact"
Odd Ratio của hình thức làm việc từ xa so với trực tiếp là 0.6910, tức là odd mức lương cao của nhóm làm việc từ xa chỉ bằng khoảng 69.10% so với nhóm làm việc trực tiếp.
Khoảng tin cậy 95% cho Odd Ratio là từ 0.6228 đến 0.7666, giá trị p là 2.419623e-12, nhỏ hơn 0.05, cho thấy có sự khác biệt có ý nghĩa thống kê giữa odd của hai nhóm làm việc trực tiếp và từ xa.
Để tính tỷ lệ odds cho nhóm làm việc từ xa so với nhóm làm việc trực
tiếp ở mức lương thấp, ta có thể thay đổi tham số rev = 'c'
để chuyển đổi giữa mức lương cao và mức lương thấp.
epitab(odd4, method = 'oddsratio', rev = 'c' )
## $tab
##
## cao p0 thap p1 oddsratio lower upper p.value
## 0 2361 0.7443253 2611 0.6679458 1.000000 NA NA NA
## 100 811 0.2556747 1298 0.3320542 1.447248 1.304431 1.605701 2.419623e-12
##
## $measure
## [1] "wald"
##
## $conf.level
## [1] 0.95
##
## $pvalue
## [1] "fisher.exact"
Kết quả cho thấy Odd ratio của kỹ sư làm việc từ xa so với kỹ sư làm việc trực tiếp là 1.447248. Điều này có nghĩa là xác suất kỹ sư làm việc từ xa có mức lương thấp cao hơn 1.447248 lần so với kỹ sư làm việc trực tiếp.
Giá trị p là 2.419623e-12, nhỏ hơn 0.05, cho thấy có sự khác biệt có ý nghĩa thống kê giữa odds của hai nhóm này.
Khoảng ước lượng cho Odd ratio với độ tin cậy 95% là từ 1.304431 đến 1.605701, củng cố thêm kết luận trên.
Phân tích mức lương trung và thấp giữa nhóm kỹ sư làm việc trực tiếp và kỹ sư làm việc từ xa
odd4 <- odd1[c("0", "100"), c("thap", "trung")]
epitab(odd4, method = 'oddsratio')
## $tab
##
## thap p0 trung p1 oddsratio lower upper p.value
## 0 2611 0.6679458 6146 0.6677532 1.000000 NA NA NA
## 100 1298 0.3320542 3058 0.3322468 1.000869 0.9244282 1.08363 1
##
## $measure
## [1] "wald"
##
## $conf.level
## [1] 0.95
##
## $pvalue
## [1] "fisher.exact"
Kết quả Odds ratio là 1.000869, nghĩa là odds của kỹ sư có mức lương thấp trong nhóm làm việc từ xa gần như bằng với odds của nhóm làm việc trực tiếp. Khoảng tin cậy 95% của OR là từ 0.9244282 đến 1.08363, và giá trị p là 1, lớn hơn 0.05. Điều này cho thấy không có sự khác biệt có ý nghĩa thống kê giữa odds của hai nhóm này.
epitab(odd4, method = 'oddsratio', rev = 'c')
## $tab
##
## trung p0 thap p1 oddsratio lower upper p.value
## 0 6146 0.6677532 2611 0.6679458 1.0000000 NA NA NA
## 100 3058 0.3322468 1298 0.3320542 0.9991321 0.9228242 1.08175 1
##
## $measure
## [1] "wald"
##
## $conf.level
## [1] 0.95
##
## $pvalue
## [1] "fisher.exact"
Dựa vào kết quả phân tích, ta có các thông số liên quan đến mức lương thấp và trung bình giữa kỹ sư làm việc trực tiếp và kỹ sư làm việc từ xa:
Odd ratio: 1.000869, nghĩa là odds của kỹ sư có mức lương thấp so với mức lương trung bình trong nhóm làm việc từ xa gần như bằng odds của nhóm làm việc trực tiếp. Khoảng tin cậy 95%: từ 0.9244282 đến 1.08363. Giá trị p: 1, không có sự khác biệt có ý nghĩa thống kê giữa odds của hai nhóm này.
Điều này cho thấy rằng không có sự khác biệt đáng kể về xác suất đạt mức lương thấp so với mức lương trung bình giữa kỹ sư làm việc trực tiếp và kỹ sư làm việc từ xa.
Kết quả bên dưới được sử dụng để kiểm tra tính độc lập giữa hai biến hình thức làm việc và mức lương bằng phương pháp Chi bình phương. Giả thuyết Ho được đưa ra là biến mức lương và hình thức làm việc độc lập.
chisq.test(table(d$remote,d$sal))
##
## Pearson's Chi-squared test
##
## data: table(d$remote, d$sal)
## X-squared = 398.3, df = 4, p-value < 2.2e-16
Kết quả kiểm định cho thấy giá trị p-value là 2.26e^-16 nhỏ hơn 0.05, có cơ sở bác bỏ H0, nghĩa là mức lương và hình thức làm việc có liên quan với nhau.
Ước lượng tỷ lệ cho mức lương trung bình
Trước tiên ta cần lọc kết quả về mức lương theo hình thức làm việc trực tiếp (0). Kết quả thu được như sau:
remote_0 <- table(d[d$remote == "0", ]$sal)
remote_0
##
## cao trung thap
## 2361 6146 2611
prop.test(remote_0["trung"], sum(remote_0), p = 0.7)
##
## 1-sample proportions test with continuity correction
##
## data: remote_0["trung"] out of sum(remote_0), null probability 0.7
## X-squared = 1146.5, df = 1, p-value < 2.2e-16
## alternative hypothesis: true p is not equal to 0.7
## 95 percent confidence interval:
## 0.5434935 0.5620644
## sample estimates:
## p
## 0.5527973
Dựa trên kết quả kiểm định tỷ lệ một mẫu với giả thuyết không (H0): Tỷ lệ người có mức lương trung trong nhóm làm việc trực tiếp là 0.7 (70%).
Kết quả kiểm định cho thấy giá trị p-value < 2.2e-16, rất nhỏ hơn giá trị α = 0.05. Điều này chứng tỏ có bằng chứng mạnh mẽ để bác bỏ giả thuyết rằng tỷ lệ người có mức lương trung trong nhóm làm việc trực tiếp là 70%.
Tỷ lệ thực tế của kỹ sư có mức lương trung trong nhóm làm việc trực tiếp là khoảng 55.28%. Khoảng ước lượng tỷ lệ ở độ tin cậy 95% là từ 54.35% đến 56.21%.
Ước lượng tỷ lệ cho nhóm lương thấp
prop.test(remote_0["thap"], sum(remote_0), p = 0.2)
##
## 1-sample proportions test with continuity correction
##
## data: remote_0["thap"] out of sum(remote_0), null probability 0.2
## X-squared = 84.149, df = 1, p-value < 2.2e-16
## alternative hypothesis: true p is not equal to 0.2
## 95 percent confidence interval:
## 0.2270129 0.2428602
## sample estimates:
## p
## 0.2348444
Vì giá trị p rất nhỏ (p < 2.2e-16), chúng ta bác bỏ giả thuyết H0. Điều này cho biết tỷ lệ thực tế của người có mức lương thấp trong nhóm kỹ sư làm việc trực tiếp không phải là 20%.
Kết quả cho thấy tỷ lệ ước tính của người có mức lương thấp trong nhóm là khoảng 23.48%.
Khoảng tin cậy 95% cho phép ước tính rằng tỷ lệ này dao động từ 22.70% đến 24.27%.
Ước lượng tỷ lệ cho nhóm lương cao
prop.test(remote_0["cao"], sum(remote_0), p = 0.3)
##
## 1-sample proportions test with continuity correction
##
## data: remote_0["cao"] out of sum(remote_0), null probability 0.3
## X-squared = 406.24, df = 1, p-value < 2.2e-16
## alternative hypothesis: true p is not equal to 0.3
## 95 percent confidence interval:
## 0.2048119 0.2201047
## sample estimates:
## p
## 0.2123583
Kết quả cho thấy p_value < 2.2e-16, cho thấy có đủ bằng chứng để bác bỏ giả thuyết H0.Điều này có nghĩa tỷ lệ kỹ sư có mức lương cao ở hình thức làm việc trực tiếp sẽ không bằng 30%.
Thực tế, tỷ lệ kỹ sư ở mức lương cao chỉ ở mức 21.24%
Khoảng tin cậy 95%:Kết quả này chỉ ra rằng tỷ lệ kỹ sư có mức lương cao trong nhóm làm việc trực tiếp không đạt mức kỳ vọng là 30%, với mức tin cậy 95% tỷ lệ kỹ sư có mức lương cao chỉ từ 20.48% đến 22.01%.
Xác suất tuyến tính
Phần này sẽ thực hiện hồi quy xác suất tuyến tính giữa biến độc lập là hình thức làm việc và biến phụ thuộc là lương. Biến lương sal sẽ được mã hoá thành biến nhị phân với giá trị 1 tương ứng với mức lương trung và cao, ngược lại giá trị 0 tương ứng với mức lương thấp.
d$sal2 <- ifelse(d$sal %in% c("trung", "cao"), 1, 0)
d$remote <- relevel(as.factor(d$remote), ref = "0")
glm <- glm(sal2 ~ remote, data = d)
summary(glm)
##
## Call:
## glm(formula = sal2 ~ remote, data = d)
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.765156 0.004052 188.829 <2e-16 ***
## remote50 -0.504111 0.027378 -18.413 <2e-16 ***
## remote100 -0.016365 0.007194 -2.275 0.0229 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 0.182553)
##
## Null deviance: 3079.8 on 16533 degrees of freedom
## Residual deviance: 3017.8 on 16531 degrees of freedom
## AIC: 18807
##
## Number of Fisher Scoring iterations: 2
\(\hat{\pi} = 0.765156 - 0.504111 \cdot remote50 - 0.016365 \cdot remote100\)
Tất cả các hệ số đều có giá trị \(p\) nhỏ hơn mức ý nghĩa \(\alpha = 0.05\), ngụ ý rằng tất cả các hệ số đều có ảnh hưởng ý nghĩa thống kê đến mức lương.
Hệ số \(\hat{\beta}_0 = 0.765156\) có ý nghĩa rằng tỷ lệ để người kỹ sư đạt được mức lương trung bình là 76.51% trong trường hợp không xuất hiện các hình thức làm việc khác.
Đối với hệ số \(\hat{\beta}_1 = - 0.504111\), ta nhận thấy khi kỹ sư làm việc theo hình thức vừa trực tiếp và trực tuyến khi không bị ảnh hưởng bởi các yếu tố khác thì chênh lệch tỷ lệ đạt mức lương trung trở lên giữa nhóm làm theo hình thức trực tiếp và trực tuyến so với và nhóm làm việc trực tiếp là -50.4%, hay có thể nói cách khác tỷ lệ bình quân để đạt mức lương trung trở lên trong nhóm này là 49.6%.
Tương tự, \(\hat{\beta}_2 = -0.016365\), ta nhận thấy khi kỹ sư làm việc theo hình thức trực tuyến và không bị ảnh hưởng bởi các yếu tố khác thì chênh lệch tỷ lệ đạt mức lương trung trở lên giữa nhóm làm việc tự do và nhóm toàn thời gian là -1,6%, hay tỷ lệ đạt mức lương trung trở lên trong nhóm này là 0%.
Ta nhận thấy ở cả ba nhóm hợp đồng, bán thời gian, và tự do đều có mối quan hệ âm với mức lương từ tầm trung trở lên. Nhưng nếu không có sự tác động của ba nhóm trên thì xác suất để đạt mức lương trên là khá cao, khoảng 0.754649.
n1 <- data.frame(remote = c("0", "50", "100"))
predict.glm(glm, newdata = n1, type= "response")
## 1 2 3
## 0.7651556 0.2610442 0.7487904
Dựa vào kết quả ta có thể nhận thấy xác suất dự báo cho biểu hiện kỹ sư làm việc trực tiếp (0) là khoảng 76,51%. Điều này có nghĩa một kỹ sư làm việc trực tiếp có xác suất khoảng 76.51% có mức lương trung đến cao.
Bên cạnh đó kết quả dự báo đối với nhóm kỹ sư vừa làm việc trực tiếp vừa làm trực tuyến (2) thì xác suất để đạt được mức lương trung trở lên chỉ khoảng 26,1%.
Trong nhóm kỹ sư làm việc trực tuyến (3) thì xác suất để đạt mức lương trung khoảng 74.88%.
Mô hình logit
Ngoài ra chúng ta có thể thực hiện mô hình hồi quy khi có hàm liên kết dạng $ (g() = logit() = () $ khi phân tích mối quan hệ giữa nhóm việc làm với mức lương, trong đó nhóm công việc trực tiếp sẽ làm nhóm cơ sở.
logit1 <- glm(factor(sal2) ~ remote, data = d, family = binomial(link = 'logit'))
summary(logit1)
##
## Call:
## glm(formula = factor(sal2) ~ remote, family = binomial(link = "logit"),
## data = d)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.18116 0.02237 52.794 <2e-16 ***
## remote50 -2.22170 0.14601 -15.216 <2e-16 ***
## remote100 -0.08898 0.03911 -2.275 0.0229 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 18506 on 16533 degrees of freedom
## Residual deviance: 18231 on 16531 degrees of freedom
## AIC: 18237
##
## Number of Fisher Scoring iterations: 4
Ta có kết quả từ hàm hồi quy thông qua mô hình logit như sau:
\[ \ln\left(\frac{\hat{\pi}}{1 - \hat{\pi}}\right) = 1.18116 - 2.22170 \cdot remote50 - 0.08898 \cdot remote100 \]
Hệ số chặn là 1.18116 điều này có nghĩa là khi không có sự tác động của hình thức làm việc khác thì \(\ln\left(\frac{\hat{\pi}}{1 - \hat{\pi}}\right)\) của kỹ sư làm việc trực tiếp đạt mức lương trung là 1.18116.
Đối với hệ số \(\hat{\beta}_1 = -2.22170\), điều này cho thấy chênh lệch \(\ln\left(\frac{\hat{\pi}}{1 - \hat{\pi}}\right)\) đạt mức lương trung trở lên giữa nhóm vừa làm việc trực tiếp và làm trực tuyến so với nhóm là trực tiếp là -1.41125.
Đối với hệ số \(\hat{\beta}_2 = -0.08898\), điều này cho thấy chênh lệch \(\ln\left(\frac{\hat{\pi}}{1 - \hat{\pi}}\right)\) đạt mức lương trung trở lên giữa nhóm làm việc trực tuyến với nhóm làm việc trực tiếp là -0.08898
Sau đó từ mô hình logit ta dự báo xác suất của biến lương tương ứng với mỗi hình thức làm việc việc. Ta cần một khung dữ liệu n với một cột type chứa các giá trị “0”, “50”, “100”.
n2 <- data.frame(remote = c("0", "50", "100") )
predict.glm(logit1, newdata = n2, type="response")
## 1 2 3
## 0.7651556 0.2610442 0.7487904
Dựa vào kết quả ta có thể nhận thấy xác suất dự báo cho biểu hiện 1 là khoảng 76.51%. Điều này có nghĩa một kỹ sư làm việc trực tiếp có xác suất khoảng 76.51% có mức lương trung đến cao.
Bên cạnh đó kết quả dự báo đối với nhóm kỹ sư làm việc vừa trực tiếp vừa trực tuyến thì xác suất để đạt được mức lương trung trở lên chỉ khoảng 25,1%.
Cuối cùng nhóm có xác suất cao nhất để đạt mức lương trung là nhóm kỹ sư làm việc theo hình thức trực tuyến xác suất để đạt được mức lương trung trở lên khoảng 74.88%.
Mô hình Probit
Thay vì sử dụng mô hình logit, ta có thể thay thế bằng mô hình probit dự đoán xác suất của một sự kiện xảy ra. Mô hình probit sử dụng hàm phân phối chuẩn. Trong đó hàm liên kết có dạng \(\text{probit}(\pi) = \Phi^{-1}(\pi)\) và hàm xác suất là \(\pi(x) = \Phi(\beta_0 + \beta_1 X + \beta_2 X_2 + \ldots + \beta_n X_n)\).
probit1 <- glm(factor(sal2) ~ remote, data = d, family = binomial(link = 'probit'))
summary(probit1)
##
## Call:
## glm(formula = factor(sal2) ~ remote, family = binomial(link = "probit"),
## data = d)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.72299 0.01309 55.244 <2e-16 ***
## remote50 -1.36312 0.08663 -15.736 <2e-16 ***
## remote100 -0.05230 0.02302 -2.272 0.0231 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 18506 on 16533 degrees of freedom
## Residual deviance: 18231 on 16531 degrees of freedom
## AIC: 18237
##
## Number of Fisher Scoring iterations: 4
Ta có kết quả mô hình hồi quy như sau:
\[ \hat{\pi} = \Phi(0.72299 - 1.36312 \times remote50 - 0.05230 \times remote100) \]
Hệ số chặn là 0.72299 và có ý nghĩa thống kê khi giá trị p_value nhỏ hơn mức ý nghĩa 5%. Điều này có nghĩa là khi không có sự tác động các biến độc lập khác thì xác suất để kỹ sư đạt mức lương trung là \(\Phi(0.72299)\).
Hệ số \(\hat{\beta}_1 = -1.36312\) có ý nghĩa là trong trường hợp các yếu tố khác không đổi thì chênh lệch tỷ lệ đạt mức lương trung giữa nhóm làm việc vừa trực tiếp vừa trực tuyến so với nhóm làm việc trực tiếp là \(\Phi(-0.86921)\).
Hệ số \(\hat{\beta}_2 = -0.05230\) có ý nghĩa là trong trường hợp các yếu tố khác không đổi thì chênh lệch tỷ lệ đạt mức lương trung giữa nhóm làm việc trực tuyến so với nhóm làm việc trực tiếp là \(\Phi(-5.30055)\).
Từ kết quả trên, em sẽ thực hiện dự báo cho tỷ lệ để đạt được mức lương trung trong các nhóm công việc. Kết quả được trình bày trong bảng bên dưới.
predict.glm(probit1, newdata = n2, type="response")
## 1 2 3
## 0.7651556 0.2610442 0.7487904
Dựa vào kết quả ta có thể nhận thấy xác suất dự báo cho biểu hiện kỹ sư làm việc trực tiếp (1) là khoảng 76,51%. Điều này có nghĩa một kỹ sư làm việc ttrực tiếp có xác suất khoảng 76,51% có mức lương trung đến cao.
Bên cạnh đó kết quả dự báo đối với nhóm kỹ sư làm việc vừa trực tiếp vừa trực tuyến (2) thì xác suất để đạt được mức lương trung trở lên chỉ khoảng 26,1%
Trong nhóm kỹ sư làm việc trực tuyến (3) thì xác suất để đạt mức lương trung khoảng 74.88%, cao hơn so với xác suất của nhóm (2).
AIC(glm)
## [1] 18806.84
AIC(logit1 )
## [1] 18236.83
AIC(probit1 )
## [1] 18236.83
Kết quả cho thấy hình Logit và Probit đều tốt hơn mô hình xác suất tuyến tính trong việc giải thích mối quan hệ giữa loại công việc và xác suất đạt mức lương từ trung trở lên, vì chỉ số AIC của của chúng thấp hơn. Ngoài ra Logit và Probit có hiệu suất tương đương nhau vì chúng có cùng chỉ số AIC
Hệ số Brier Score
Tương tự, phần này sẽ tính chỉ số Brier Score được sử dụng để đo lường độ chính xác của các dự đoán xác suất.
BrierScore(logit1)
## [1] 0.1825199
BrierScore(probit1)
## [1] 0.1825199
Kết quả bên trên cho thấy Brier Score cho logit và probit đều là 0.18525199. Điều này có nghĩa là cả hai mô hình đều có hiệu suất tương đương khi dự đoán xác suất cho biến mục tiêu. Brier Score khá thấp (gần 0), cho thấy rằng cả hai mô hình đều dự đoán khá tốt và không sai lệch nhiều so với giá trị thực tế.
Ma trận nhầm lẫn
Conf(table(predict(probit1, type="response") >=0.5,probit1$data$sal2 == '1'))
##
## Confusion Matrix and Statistics
##
## Reference
## Prediction FALSE TRUE
## FALSE 184 65
## TRUE 3909 12376
##
## Total n : 16'534
## Accuracy : 0.7596
## 95% CI : (0.7531, 0.7661)
## No Information Rate : 0.7524
## P-Value [Acc > NIR] : 0.0161
##
## Kappa : 0.0580
## Mcnemar's Test P-Value : < 2.2e-16
##
## Sensitivity : 0.0450
## Specificity : 0.9948
## Pos Pred Value : 0.7390
## Neg Pred Value : 0.7600
## Prevalence : 0.2476
## Detection Rate : 0.0151
## Detection Prevalence : 0.0111
## Balanced Accuracy : 0.5199
## F-val Accuracy : 0.0848
## Matthews Cor.-Coef : 0.1408
##
## 'Positive' Class : FALSE
Từ kết quả trên ma trận cho ta biết có 184 số lần mô hình dự đoán đúng là mức lương thấp. Tuy nhiên có 65 lần mô hình dự đoán sai đối với mức lương trung và khoảng 3909 lần mô hình dự đoán sai mức lương thấp. Số lần mô hình dự đoán đúng ở mức lương trung là 12376 lần đạt cao nhất.
Ngoài ra ta có tỷ lệ dự đoán đúng của mô hình là 0.7596 và ở độ tin cậy là 95% thì tỷ lệ này nằm trong khoảng từ (0.7531, 0.7661)
Trong phần này, em lập bảng tần số cho biến trình độ công việc để xem xét và bóc tách các cấu phần có trong biến này. Trước khi phân tích, ta tiến hành đổi biến experience_level thành el. Sau đó em tiến hành lập bảng tần số cho biến el. Kết quả thu được như sau:
d <- d %>% rename(el = experience_level) # đổi tên biến
table(d$el)
##
## EN EX MI SE
## 1325 501 4038 10670
e <- as.data.frame(table(d$el))
colnames(e) <- c("Exp", "Count")
ggplot(e, aes(x = Exp, y = Count, fill = Exp)) +
geom_bar(stat = "identity", width = 0.5) +
geom_text(aes(label = Count), vjust = -0.5, size = 3.5) + # Hiển thị số liệu bên trên cột
labs(title = "Đồ thị trình độ chuyên môn ", x = "Exp", y = "Frequency") +
theme_minimal()
Dựa vào kết quả từ bảng và đồ thị ta nhận thấy ta thấy kỹ sư cao cấp (SE) nhiều nhất với 10670 người, kế đến là kỹ sư trung cấp (MI) với số lượng là 4038 người, tiếp theo là kỹ sư mới vào nghề (EN) với 1325 người, cuối cùng là kỹ sư điều hành với 501 người.
Mặt khác ta cũng cần thực hiện phân tích tỷ lệ phần trăm đối với từng nhóm kinh nghiệm công việc và trực quan chúng bằng đồ thị tròn để có thể dễ dàng so sánh với nhau.
Tiếp theo ta có thể dựa vào bảng tần suất để phân tích tỷ trọng
prop.table(table(d$el))*100
##
## EN EX MI SE
## 8.01379 3.03012 24.42240 64.53369
e$Percentage <- (e$Count / sum(e$Count)) * 100
ggplot(e, aes(x = "", y = Percentage, fill = Exp)) +
geom_bar(width = 1, stat = "identity") +
coord_polar("y") +
labs(title = "Tỷ lệ phần trăm cho nhóm kinh nghiệm làm việc", x = "", y = "") +
theme_void() +
theme(legend.title = element_blank()) +
geom_text(aes(label = paste(round(Percentage, 2), "%")),
position = position_stack(vjust = 0.5), size = 4)
Dựa vào kết quả bảng tần suất và đồ thị phía trên ta thấy kỹ sư cao cấp chiếm tỷ trọng lớn nhất với 64.5%, kỹ sư trung cấp đứng thứ hai với 24.42%, tỷ trọng kỹ sư mới vào nghề chiếm 8.01% và tỷ trọng kỹ sư điều hành chiếm 3.03%.
Trước khi phân tích những yếu tố sẽ tác động đến tiền lương của những người kĩ sư chúng ta cần phân nhóm thành ba mức lương là thấp, trung và cao để dễ dàng nhận diện, so sánh. Hơn thế, chúng ta có thể phân tích sâu sắc về cấu trúc tiên lương cũng như là các yếu tố định tính sẽ ảnh hưởng đáng kể đến các mức phân loại tiền lương của các kỹ sư một cách hiệu quả nhất.
d$salusd <- cut(d$salary_in_usd,
breaks = c(0, 100000,200000, Inf),
labels = c("thap","trung", "cao"))
table(d$el, d$salusd) %>% addmargins()
##
## thap trung cao Sum
## EN 894 393 38 1325
## EX 38 244 219 501
## MI 1631 2008 399 4038
## SE 1530 6613 2527 10670
## Sum 4093 9258 3183 16534
table_data <- table(d$el, d$salusd)
# Chuyển bảng chéo sang định dạng dài
long_data <- as.data.frame(table_data)
colnames(long_data) <- c("Exp", "salary", "count")
# Vẽ biểu đồ cột
ggplot(long_data, aes(x = salary, y = count, fill = salary)) +
geom_bar(stat = "identity", position = position_dodge()) +
geom_text(aes(label = count), vjust = 1, size = 2.5, position = position_dodge(1)) +
facet_wrap(~ Exp, scales = "free_y", nrow = 2) +
labs(title = "Phân bố mức lương theo trình độ công việc",
x = "Mức lương",
y = "Số lượng kỹ sư",
fill = "Mức lương") +
theme_minimal()
– Cụ thể,Kỹ sư cao cấp (SE) có 1530 người nhận lương ở mức thấp, 6613 người nhận lương ở mức trung bình và 2527 người nhận lương ở mức cao, tổng cộng là 10670 người.Đối với các kỹ sư mới vào nghề (EN), có 894 người nhận lương ở mức thấp, 393 người nhận lương ở mức trung bình và 38 người nhận lương ở mức cao, với tổng số là 1325 người.
– Đối với các kỹ sư điều hành (EX), có 38 người nhận lương ở mức thấp, 244 người ở mức trung bình và 219 người ở mức cao, tổng cộng là 501 người. Các kỹ sư trung cấp (MI) có 1631 người nhận lương ở mức thấp, 2008 người ở mức trung bình và 399 người ở mức cao, với tổng cộng 4038 người.
– Tổng cộng, bảng này cho thấy có 4093 kỹ sư ở mức lương thấp, 9258 kỹ sư ở mức lương trung bình và 3183 kỹ sư ở mức lương cao, với tổng số kỹ sư là 16534 người.
Tiếp theo, lập bảng tần suất giữa hai biến lương và hình thức công việc. Từ đó vẽ đồ thì cột theo từng nhóm công việc như bên dưới đây.
table_data <- round(prop.table(table_data)*100,2)
table_data
##
## thap trung cao
## EN 5.41 2.38 0.23
## EX 0.23 1.48 1.32
## MI 9.86 12.14 2.41
## SE 9.25 40.00 15.28
# Chuyển đổi dữ liệu thành dạng dài để dễ vẽ biểu đồ
table_data <- as.data.frame.table(table_data)
colnames(table_data) <- c("Group", "Level", "Percentage")
# Vẽ biểu đồ
ggplot(table_data, aes(x = Level, y = Percentage, fill = Level)) +
geom_bar(stat = "identity", position = "dodge") +
geom_text(aes(label = Percentage), vjust = 1, size = 2.5, position = position_dodge(1)) +
facet_wrap(~ Group, scales = "free_y", nrow = 2) +
labs(title = "Đồ thị phần trăm trình độ chuyên môn ",
x = "Trình độ chuyên môn",
y = "Phần trăm",
fill = "Mức lương") +
scale_fill_manual(values = c("thap" = "blue", "trung" = "orange", "cao" = "red")) +
theme_minimal()
Đối với nhóm kỹ sư cao cấp: Tỷ lệ mức lương thấp chiếm đạt 9.25% so với tổng số kỹ sư, tỷ lệ kỹ sư mức lương trung là 40% và tỷ lệ kỹ sư đạt mức lương cao đối với kỹ sư cao cấp là 15.28%.
Đối với nhóm kỹ mới vào nghề: Tỷ lệ mức lương thấp chiếm 5.41% so với tổng số kỹ sư, vì thế ta nhận thấy mức lương thấp đứng cao nhất trong nhóm kỹ sư mới vào nghề. Tỷ lệ kỹ sư có mức lương trung là 2.38% và thấp nhất trong nhóm là tỷ lệ kỹ sư đạt mức lương cao đạt 0.23%.
Đối với nhóm kỹ sư điều hành: Tỷ lệ mức lương thấp chiếm chỉ 0.23% so với tổng số kỹ sư. Tỷ lệ kỹ sư có mức lương trung là 1.48% và tỷ lệ kỹ sư đạt mức lương cao đạt 1.32%.
Đối với nhóm kỹ sư trung cấp: Tỷ lệ mức lương thấp chiếm chỉ 9.86% so với tổng số kỹ sư. Tỷ lệ kỹ sư có mức lương trung là 12.14% và tỷ lệ kỹ sư đạt mức lương cao đạt 2.41%.
Trong phần này, em thực hiện tính relative risk giữa trình độ chuyên môn tương ứng với các mức lương thấp trung và cao.trong đó biến độc lập là trình độ chuyên môn và biếnbiến phụ thuộc là mức lương. Trình độ chuyên môn được chọn làm tham chiếu là SE - Kỹ sư cao cấp.
# Tạo bảng dữ liệu
nth <- data.frame(
level = c("EN", "EX", "MI", "SE"),
thap = c(894, 38, 1631, 1530),
trung = c(393, 244, 2008, 6613),
cao = c(38, 219, 399, 2527),
sum = c(1325, 501, 4038, 10670)
)
# Tính tỷ lệ cho mỗi mức lương
nth$rthap <- nth$thap / nth$sum
nth$rtrung <- nth$trung / nth$sum
nth$rcao <- nth$cao / nth$sum
# Tính Relative Risk (RR) so với nhóm SE (Kỹ sư cao cấp)
rlow <- nth$rthap / nth$rthap[nth$level == "SE"]
rme <- nth$rtrung / nth$rtrung[nth$level == "SE"]
rhigh <- nth$rcao / nth$rcao[nth$level == "SE"]
# Tạo bảng kết quả RR
rr1 <- data.frame(
level = nth$level,
RR_low = rlow,
RR_me = rme,
RR_high = rhigh
)
# Hiển thị kết quả
print(rr1)
## level RR_low RR_me RR_high
## 1 EN 4.7053792 0.4785668 0.1210952
## 2 EX 0.5289552 0.7858108 1.8457189
## 3 MI 2.8168300 0.8023490 0.4172206
## 4 SE 1.0000000 1.0000000 1.0000000
Giải thích
– Lương thấp:
Kỹ sư mới vào nghề (EN) có nguy cơ nhận lương thấp cao hơn 4.712 lần so với kỹ sư cao cấp (SE). Kỹ sư điều hành (EX) có nguy cơ nhận lương thấp thấp hơn, chỉ bằng 0.530 lần so với kỹ sư cao cấp (SE). Kỹ sư trung cấp (MI) có nguy cơ nhận lương thấp cao hơn 2.828 lần so với kỹ sư cao cấp (SE).
– Lương trung:
Kỹ sư mới vào nghề (EN) có nguy cơ nhận lương trung thấp hơn, chỉ bằng 0.479 lần so với kỹ sư cao cấp (SE). Kỹ sư điều hành (EX) và kỹ sư trung cấp (MI) có nguy cơ nhận lương trung thấp hơn một chút so với kỹ sư cao cấp (SE), với RR lần lượt là 0.785 và 0.802.
– Lương cao:
Kỹ sư mới vào nghề (EN) có nguy cơ nhận lương cao thấp hơn nhiều, chỉ bằng 0.122 lần so với kỹ sư cao cấp (SE). Kỹ sư điều hành (EX) có nguy cơ nhận lương cao cao hơn, với RR là 1.844 so với kỹ sư cao cấp (SE). Kỹ sư trung cấp (MI) có nguy cơ nhận lương cao thấp hơn, với RR là 0.418 so với kỹ sư cao cấp (SE). Tóm lại, kỹ sư mới vào nghề (EN) có xu hướng nhận lương thấp hơn nhiều so với các cấp bậc khác, trong khi kỹ sư điều hành (EX) có nguy cơ nhận lương cao hơn so với kỹ sư cao cấp (SE).
Trong phần này em sẽ thực hiện khoảng ước lượng cho RR với khoảng tin cậy 95%. Tuy nhiên để cho kết quả dễ dàng phân tích hơn trước nên em sẽ thực hiện trên bảng có số lượng nhỏ hơn là 2x2.
Trước tiên em tính khoảng ước lượng Relative Risk đối với mức lương thấp và trung nhưng trong nhóm SE và MI
NTH1 <- matrix(c(1530, 3313,
1631, 2008),
nrow = 2,
byrow = TRUE,
dimnames = list(c( "SE", "MI"),
c("thap", "trung")))
library("epitools")
epitab(NTH1, method = 'riskratio', rev = 'c')
## $tab
## trung p0 thap p1 riskratio lower upper p.value
## SE 3313 0.6840801 1530 0.3159199 1.000000 NA NA NA
## MI 2008 0.5517999 1631 0.4482001 1.418714 1.342887 1.498823 1.68852e-35
##
## $measure
## [1] "wald"
##
## $conf.level
## [1] 0.95
##
## $pvalue
## [1] "fisher.exact"
## $tab
Xác suất để có mức lương thấp đối với nhóm kỹ sư trung cấp là khoảng 44.82% và xác suất để đạt mức lương thấp đối với nhóm kỹ sư cao cấp là 31.45%, trong trường hợp này tỷ lệ Relative Risk (RR) là khoảng 1.418 điều này cho thấy khả năng để nhóm kỹ sư trung cấp đạt được mức lương “thấp” cao gấp 1.418 lần đối với nhóm kỹ sư cao cấp. Khoảng ước lượng Relative Risk có kết quả là 1.34 và 1.49 đối với độ tin cậy 95%. Giá trị p.value khoảng 1.68852e-35 nhỏ hơn mức ý nghĩa 5% vì vậy kết quả này có ý nghĩa thống kê cho thấy có sự khác biệt ở nhóm kỹ sư trung cấp và kỹ sư cao cấp.
Trong phần này em sẽ tính tỷ lệ Odd Ratio của hai biến lương và trình độ chuyên môn. Để dễ dàng trong phân tích em sẽ tách số liệu từ bảng tần số ban đầu thành bảng tần số mới chỉ bao gồm nhóm lương thấp, cao và nhóm KỸ SƯ cao cấp và kỹ sư trung cấp.
nth_odd <- table(d$el,d$sal)
nth_odd2 <- nth_odd[c("SE", "MI"), c("thap", "cao")]
epitab(nth_odd2, method = 'oddsratio')
## $tab
##
## thap p0 cao p1 oddsratio lower upper p.value
## SE 1530 0.484024 2527 0.8636364 1.0000000 NA NA NA
## MI 1631 0.515976 399 0.1363636 0.1481171 0.1305109 0.1680983 5.917836e-228
##
## $measure
## [1] "wald"
##
## $conf.level
## [1] 0.95
##
## $pvalue
## [1] "fisher.exact"
Kết quả Odd Ratio của MI so với SE khoảng 0.148 điều này có nghĩa là odd mức lương cao của nhóm kỹ sư trung cấp chỉ bằng khoảng 14,8%% so với odd kỹ sư cao cấp.
Khoảng tin cậy của OR đối với độ tin cậy 95% là từ 0.1305 đến 0.168 và giá trị p là 5.917836e-2285, nhỏ hơn 0.05 cho thấy có sự khác biệt giữa odd ở hai nhóm MI và SE.
Phần này sẽ thực hiện hồi quy xác suất tuyến tính giữa biến độc lập là trình độ kỹ sư và biến phụ thuộc là lương. Tuy nhiên trong phần này ta cần thay đổi biến sal thành biến nhị phân với giá trị 1 tương ứng với mức lương trung và cao, ngược lại giá trị 0 tương ứng với mức lương thấp.
d$sal2 <- ifelse(d$sal %in% c("trung", "cao"), 1, 0)
d$el <- relevel(as.factor(d$el), ref = "SE")
glm3 <- glm(sal2 ~ el, data = d)
summary(glm3)
##
## Call:
## glm(formula = sal2 ~ el, data = d)
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.856607 0.003846 222.733 <2e-16 ***
## elEN -0.531324 0.011571 -45.917 <2e-16 ***
## elEX 0.067544 0.018160 3.719 2e-04 ***
## elMI -0.260520 0.007340 -35.494 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 0.157819)
##
## Null deviance: 3079.8 on 16533 degrees of freedom
## Residual deviance: 2608.7 on 16530 degrees of freedom
## AIC: 16401
##
## Number of Fisher Scoring iterations: 2
\(\hat\pi= 0.856607-0.531324EN+0.067544EX-0.260520MI\)
Tất cả các hệ số đều có giá trị p nhỏ hơn mức ý nghĩa \(\alpha= 5%\) , ngụ ý rằng tất cả các hệ số đều có ảnh hưởng ý nghĩa thống kê đến mức lương.
Hệ số \(\hat\beta_0 = 0.856607\) có ý nghĩa rằng tỷ lệ để người kỹ sư đạt được mức lương trung bình là 85.66% trong trường hợp không xuất hiện các nhóm công việc khác.
Đối với hệ số \(\hat\beta_1 =-0.531324\) ta nhận thấy khi kỹ sư mới vào nghề và không bị ảnh hưởng bởi các yếu tố khác thì chênh lệch tỷ lệ đạt mức lương trung trở lên giữa nhóm mới vàơ nghề và nhóm kỹ sư cao cấp là -53.13%, hay có thể nói cách khác tỷ lệ bình quân để đạt mức lương trung trở lên trong nhóm này là 32.53%
Tương tự \(\hat\beta_2 =0.067544\) ta nhận thấy khi kỹ sư làm điều hành và không bị ảnh hưởng bởi các yếu tố khác thì chênh lệch tỷ lệ đạt mức lương trung trở lên giữa nhóm là 6.75%, hay tỷ lệ đạt mức lương trung trở lên trong nhóm này là 92.41%
Cuối cùng với \(\hat\beta_3 =-0.260520\) ta nhận thấy khi kỹ sư trung câp không bị ảnh hưởng bởi các yếu tố khác thì chênh lệch tỷ lệ đạt mức lương trung trở lên giữa nhóm kỹ sư trung cấp và kỹ sư cao cấp là - 26.05%, hay tỷ lệ đạt mức lương trung trở lên trong nhóm này là 59.61%
Ta nhận thấy ở cả ba nhóm kỹ sư trung cấp, mới vào nghề, đều có mối quan hệ âm với mức lương lương từ tầm trung trở lên, riêng kỹ sư điều hành thì dương. Nhưng nếu không có sự tác động của ba nhóm trên thì xác suất để đạt mức lương trên là khá cao khoảng 85.66%.
h <- data.frame(el = c("SE", "EX", "EN","MI") )
predict.glm(glm3, newdata = h, type="response")
## 1 2 3 4
## 0.8566073 0.9241517 0.3252830 0.5960872
Dựa vào kết quả ta có thể nhận thấy xác suất dự báo cho biểu hiện kỹ sư cao cấp (1) là khoảng 85.66%. Điều này có nghĩa một kỹ sư cao cấp có xác suất khoảng 85.66% nhận mức lương trung đến cao.
Bên cạnh đó kết quả dự báo đối với nhóm kỹ sư điều hành (2) thì xác suất để đạt được mức lương trung trở lên khoảng 92.41%, cao nhất trong nhóm kỹ sư.
Trong nhóm kỹ sư mới vào nghề (3) thì xác suất để đạt mức lương trung trở lên chỉ khoảng 32.52%, thấp nhất trong nhóm trình độ kỹ sư.
Cuối cùng kỹ sư trung cấp toàn thời gian có xác suất đạt mức lương từ trung trở lên là 59.6%, đúng thứ 3 trong nhóm trình độ kỹ sư.
nth_logit <- glm(factor(sal2) ~ el, data = d, family = binomial(link = 'logit'))
summary(nth_logit)
##
## Call:
## glm(formula = factor(sal2) ~ el, family = binomial(link = "logit"),
## data = d)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.78739 0.02762 64.708 < 2e-16 ***
## elEN -2.51699 0.06482 -38.830 < 2e-16 ***
## elEX 0.71275 0.17086 4.172 3.02e-05 ***
## elMI -1.39820 0.04233 -33.033 < 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: 18506 on 16533 degrees of freedom
## Residual deviance: 16161 on 16530 degrees of freedom
## AIC: 16169
##
## Number of Fisher Scoring iterations: 4
Ta có kết từ hàm hồi quy thông qua mô hình logit như sau:
\(ln(\frac{\hat\pi}{1−\hat\pi})= 1.78739−2.51699EN+0.71275EX−1.39820MI\)
Hệ số chặn là 1.78739 điều này có nghĩa là khi không có sự tác động của nhóm kỹ sư khác thì \(ln(\frac{\hat\pi}{1−\hat\pi})\) của kỹ sư cao cấp đạt mức lương trung và cao là 1.78739
Đối với hệ số \(\hat\beta_1 =-2.51699\) điều này cho thấy chênh lệch \(ln(\frac{\hat\pi}{1−\hat\pi})\) đạt mức lương trung trở lên giữa nhóm kỹ sư mới vào nghề và nhóm kỹ sư cao cấp là -2.51699
Tương tự, hệ số \(\hat\beta_2 =0.71275\) điều này cho thấy chênh lệch \(ln(\frac{\hat\pi}{1−\hat\pi})\) mức lương trung trở lên giữa nhóm kỹ sư điều hành và nhóm kỹ sư cao cấp là 0.71275.
Đối với hệ số \(\hat\beta_3 =-1.39820\) điều này cho thấy chênh lệch \(ln(\frac{\hat\pi}{1−\hat\pi})\) đạt mức lương trung trở lên giữa nhóm kỹ sư trung cấp và nhóm kỹ sư cao cấp là -1.39820 Sau đó từ mô hình logit ta dự báo xác suất của biến lương tương ứng với mỗi nhóm việc là. Ta cần một khung dữ liệu n với một cột type chứa các giá trị “SE”, “EX”, “EN”, và “MI”.
h <- data.frame(el = c("SE", "EX", "EN","MI") )
predict.glm(nth_logit, newdata = h, type="response")
## 1 2 3 4
## 0.8566073 0.9241516 0.3252830 0.5960872
Dựa vào kết quả ta có thể nhận thấy xác suất dự báo cho biểu hiện kỹ sư cao cấp (1) là khoảng 85.66%. Điều này có nghĩa một kỹ sư cao cấp có xác suất khoảng 85.66% nhận mức lương trung đến cao.
Bên cạnh đó kết quả dự báo đối với nhóm kỹ sư điều hành (2) thì xác suất để đạt được mức lương trung trở lên khoảng 92.41%, cao nhất trong nhóm kỹ sư.
Trong nhóm kỹ sư mới vào nghề (3) thì xác suất để đạt mức lương trung trở lên chỉ khoảng 32.52%, thấp nhất trong nhóm trình độ kỹ sư.
Cuối cùng kỹ sư trung cấp toàn thời gian có xác suất đạt mức lương từ trung trở lên là 59.6%, đúng thứ 3 trong nhóm trình độ kỹ sư.
Thay vì sử dụng mô hình logit ta có thể thay thế bằng mô hình probit dự đoán xác suất của một sự kiện xảy ra, mô hình probit sử dụng hàm phân phối chuẩn. Trong đó hàm liên kết có dạng \(probit(π) = Φ^{-1}(π)\) và hàm xác suất là \(π(x)=Φ(β_0+β_1X+β_2X_2+...+β_nX_n)\). Ta thực hiện thao tác khá tương tự với phân trên
nth_probit <- glm(factor(sal2) ~ el, data = d, family = binomial(link = 'probit'))
summary(nth_probit)
##
## Call:
## glm(formula = factor(sal2) ~ el, family = binomial(link = "probit"),
## data = d)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.06520 0.01500 71.020 < 2e-16 ***
## elEN -1.51818 0.03877 -39.163 < 2e-16 ***
## elEX 0.36836 0.08418 4.376 1.21e-05 ***
## elMI -0.82197 0.02495 -32.947 < 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: 18506 on 16533 degrees of freedom
## Residual deviance: 16161 on 16530 degrees of freedom
## AIC: 16169
##
## Number of Fisher Scoring iterations: 4
Ta có kết quả mô hình hồi quy như sau:
\(\hat\pi= \Phi( 1.06520-1.51818×EN+0.36836×EX−0.82197×MI)\)
Hệ số chặn là 1.06520 và có ý nghĩa thống kê khi giá trị p_value nhỏ hơn mức ý nghĩa 5%. Điều này có nghĩa là khi không có sự tác động các biến độc lập khác thì xác suất để kỹ sư cao cấp là \(\Phi(1.06520)\)
Hệ số \(\hat\beta_1= −1.51818\) có ý nghĩa là trong trường hợp các yếu tố khác không đổi thì chênh lệch tỷ lệ đạt mức lương trung giữa nhóm kỹ sư mới vào nghề và kỹ sư cao cấp đạt khoảng \(\Phi(-1.51818)\)
Hệ số \(\hat\beta_2=−5.30055\) có ý nghĩa là trong trường hợp các yếu tố khác không đổi thì chênh lệch tỷ lệ đạt mức lương trung giữa nhóm kỹ sư điều hành và nhóm kỹ sư cao cấp đạt khoảng \(\Phi(−5.30055)\)
Hệ số \(\hat\beta_3= −1.16870\) có ý nghĩa là trong trường hợp các yếu tố khác không đổi thì chênh lệch tỷ lệ đạt mức lương trung giữa nhóm kỹ sư trung cấp và kỹ sư cao cấp đạt khoảng \(\Phi(−1.16870)\)
Từ kết quả trên em sẽ thực hiện dự báo cho tỷ lệ để đạt được mức lương trung trong các nhóm công việc. Kết quả được trình bày trong bảng bên dưới.
predict.glm(nth_probit, newdata = h, type="response")
## 1 2 3 4
## 0.8566073 0.9241517 0.3252830 0.5960872
Dựa vào kết quả ta có thể nhận thấy xác suất dự báo cho biểu hiện kỹ sư cao cấp (1) là khoảng 85.66%. Điều này có nghĩa một kỹ sư cao cấp có xác suất khoảng 85.66% nhận mức lương trung đến cao.
Bên cạnh đó kết quả dự báo đối với nhóm kỹ sư điều hành (2) thì xác suất để đạt được mức lương trung trở lên khoảng 92.41%, cao nhất trong nhóm kỹ sư.
Trong nhóm kỹ sư mới vào nghề (3) thì xác suất để đạt mức lương trung trở lên chỉ khoảng 32.52%, thấp nhất trong nhóm trình độ kỹ sư.
Cuối cùng kỹ sư trung cấp toàn thời gian có xác suất đạt mức lương từ trung trở lên là 59.6%, đúng thứ 3 trong nhóm trình độ kỹ sư.
AIC (Akaike Information Criterion) là một chỉ số được sử dụng để đánh giá và so sánh các mô hình thống kê. Mô hình có AIC thấp hơn được coi là tốt hơn vì nó cân bằng giữa độ phù hợp và độ phức tạp của mô hình. Do đó chúng ta sử dụng chỉ số này để so sánh ba mô hình bên trên để lựa chọn mô hình phù hợp nhất trong phân tích xác suất để đạt được mức lương trung.
AIC(glm3)
## [1] 16400.63
AIC(nth_logit)
## [1] 16168.77
AIC(nth_probit)
## [1] 16168.77
kết quả cho thấy giá trị AIC của mô hình logit và mô hình probit nhỏ hơn so với mô hình glm, do đó chọn 2 mô hình logit và probit.
Tương tự, phần này sẽ tính chỉ số Brier Score được sử dụng để đo lường độ chính xác của các dự đoán xác suất.
BrierScore(nth_logit)
## [1] 0.1577808
BrierScore(nth_probit)
## [1] 0.1577808
Kết quả bên trên cho thấy Brier Score cho logit và probit đều là 0.1577808. Điều này có nghĩa là cả hai mô hình đều có hiệu suất tương đương khi dự đoán xác suất cho biến mục tiêu. Brier Score khá thấp (gần 0), cho thấy rằng cả hai mô hình đều dự đoán khá tốt và không sai lệch nhiều so với giá trị thực tế.
Conf(table(predict(nth_probit, type="response") >=0.5,nth_probit$data$sal2 == '1'))
##
## Confusion Matrix and Statistics
##
## Reference
## Prediction FALSE TRUE
## FALSE 894 431
## TRUE 3199 12010
##
## Total n : 16'534
## Accuracy : 0.7805
## 95% CI : (0.7741, 0.7867)
## No Information Rate : 0.7524
## P-Value [Acc > NIR] : < 2.2e-16
##
## Kappa : 0.2377
## Mcnemar's Test P-Value : < 2.2e-16
##
## Sensitivity : 0.2184
## Specificity : 0.9654
## Pos Pred Value : 0.6747
## Neg Pred Value : 0.7897
## Prevalence : 0.2476
## Detection Rate : 0.0801
## Detection Prevalence : 0.0541
## Balanced Accuracy : 0.5919
## F-val Accuracy : 0.3300
## Matthews Cor.-Coef : 0.2921
##
## 'Positive' Class : FALSE
Từ kết quả trên ma trận cho ta biết có 894 số lần mô hình dự đoán đúng là mức lương thấp. Tuy nhiên có 431 lần mô hình dự đoán sai đối với mức lương trung và khoảng 3199 lần mô hình dự đoán sai mức lương thấp. Số lần mô hình dự đoán đúng ở mức lương trung là 12010 lần đạt cao nhất.
Ngoài ra ta có tỷ lệ dự đoán đúng của mô hình là 0.7805 và ở độ tin cậy là 95% thì tỷ lệ này nằm trong khoảng từ (0.7741, 0.7867)`
Trong phần này em lập bảng tần số cho biến company_size nhằm mục tiêu thống kê số lượng kỹ sư trong mỗi nhóm hình thức công việc và so sánh sự khác biệt về số lượng cũng như tỷ lệ giữa chúng. Từ đó giúp cho mọi người hiểu rõ hơn ý nghĩa của biến này trước khi đánh giá sự ảnh hưởng của hình thức công việc đến mức lương.
Trong đó:
Trước tiên chúng ta cần đổi tên biến từ company_size thành size điều này sẽ giúp câu lệnh được thực hiện một cách ngắn gọn và dễ nhận biết hơn trước.
d <- d %>% rename(size = company_size) # đổi tên biến
Tiếp đến, em sẽ lập bảng tần số cho biến size, kết quả sẽ được trình bày ở dưới đây.
table(d$size)
##
## L M S
## 1040 15306 188
Bảng tần số cho biết có 188 kỹ sư làm việc ở công ty có quy mô nhỏ (S). Thêm vào đó, số lượng kỹ sư làm việc ở công ty có quy mô nhỏ có số lượng thấp nhất.
Trong khi đó, số lượng kỹ sư làm việc ở công ty có quy mô trung bình (M) có số lượng lớn nhất trong bộ dữ liệu với 15306 kỹ sư.
Cuối cùng, có 1040 kỹ sư làm việc ở công ty có quy mô lớn (L)
# Tạo bảng tần số cho biến size
size <- as.data.frame(table(d$size))
colnames(size) <- c("Size", "Count")
ggplot(size, aes(x = Size, y = Count, fill = Size)) +
geom_bar(stat = "identity", width = 0.5) +
geom_text(aes(label = Count), vjust = -0.5, size = 3.5) + # Hiển thị số liệu bên trên cột
labs(title = "DO THI QUY MO CONG TY ", x = "Size", y = "Frequency") +
theme_minimal()
Dựa vào kết quả từ bảng và đồ thị ta nhận thấy hầu hết các kỹ sư (15306 người) sẽ làm việc ở công ty có quy mô M. Ngược lại, số kỹ sư làm việc ở công ty nhỏ là ít nhất (188 người). Số lượng kỹ sư làm việc ở công ty có quy mô mô lớn là 1040 kỹ sư
Trong trường hợp này, em vẫn lập bảng tần số giữa hai biến size (quy mô công ty) và sal (mức lương phân tổ thành ba nhóm). Bảng tần số giúp chúng ta hiểu rõ hơn về sự phân bố mức lương theo từng quy mô công ty.
table(d$size, d$sal) %>% addmargins()
##
## cao trung thap Sum
## L 257 454 329 1040
## M 2920 8748 3638 15306
## S 6 56 126 188
## Sum 3183 9258 4093 16534
Vẽ đồ thị cột phân theo từng nhóm quy mô công ty để đánh giá mức lương của từng nhóm kỹ sư.
table_data5 <- table(d$size, d$sal)
# Chuyển bảng chéo sang định dạng dài
long_data5 <- as.data.frame(table_data5)
colnames(long_data5) <- c("size", "salary", "count")
# Vẽ biểu đồ cột
ggplot(long_data5, aes(x = salary, y = count, fill = salary)) +
geom_bar(stat = "identity", position = position_dodge()) +
geom_text(aes(label = count), vjust = 1, size = 2.5, position = position_dodge(1)) +
facet_wrap(~ size, scales = "free_y", nrow = 2) +
labs(title = "Phân bố mức lương theo quy mô công ty",
x = "Mức lương",
y = "Số lượng kỹ sư",
fill = "Mức lương") +
theme_minimal()
Dựa vào đồ thị chúng ta có thể thấy rằng:
Số lượng kỹ sư làm không làm việc ở công ty có quy mô lớn (L) có mức lương cao là 257 kỹ sư. Trong khi đó, số lượng kĩ sư làm việc ở công ty có quy mô lớn có mức lương trung có số lượng lớn nhất là 454 kỹ sư. Số kỹ sư có mức lương thấp làm việc ở công ty lớn là 329 kỹ sư.
Cùng với đó, số lượng kỹ sư làm việc ở công ty có quy mô trung bình (M) có mức lương thấp có số lượng là 3638 kỹ sư. Số lượng kỹ sư làm việc ở công ty quy mô trung bình có mức lương trung bình có số lượng lớn nhất trong tất cả các nhóm là 8748. Tương tự, số lượng kỹ sư làm việc ở công ty quy mô trung bình có mức lương cao có số lượng lớn nhất trong tất cả các nhóm với 2920 kỹ sư.
Số kỹ sư làm việc ở công ty có quy mô nhỏ (S) có số lượng ít nhất so với nhóm kỹ sư làm việc ở công ty có quy mô lớn và quy mô trung bình. Do đó số lượng kỹ sư có mức lương làm việc ở công ty quy mô nhỏ chỉ có 126 kỹ sư. Số lượng này giảm dần ở mức lương trung (56) kỹ sư và mức lương cao chỉ có 6 kỹ sư đạt mức lương này.
Trong phần này em sẽ thực hiện tính Relative Risk nhằm phân tích sự ảnh hưởng của từng quy mô công ty lên các mức lương khác nhau.
RRsize <- matrix(c(329,454,257,
3638, 8748,2920,
126,56, 6),
nrow = 3,
byrow = TRUE,
dimnames = list(c("L", "M", "S"),
c("thap", "trung", "cao")))
RRsize
## thap trung cao
## L 329 454 257
## M 3638 8748 2920
## S 126 56 6
Dựa trên bảng tần số đã được giải thích rõ ở phần trước em sẽ tính toán Relative Risk trong đó biến độc lập là quy mô công ty và phụ thuộc là mức lương. Nhóm công việc được chọn làm tham chiếu là M
Trước tiên cần tính tỷ lệ mức lương thấp, trung, cao trong tổng số người ở mỗi quy mô công ty khác nhau. Sau đó lấy tỷ lệ vừa tính được chia cho tỷ lệ ở nhóm tham chiếu, trong ví dụ là nhóm làm việc ở công ty có quy mô trung bình
R_tablesize <- as.data.frame(RRsize)
# Tính tổng số người trong mỗi nhóm
R_tablesize$Total <- rowSums(R_tablesize)
# Tính tỷ lệ người ở mỗi mức lương trong mỗi nhóm
R_tablesize$Rate_thap <- R_tablesize$thap / R_tablesize$Total
R_tablesize$Rate_trung <- R_tablesize$trung / R_tablesize$Total
R_tablesize$Rate_cao <- R_tablesize$cao / R_tablesize$Total
# Tính tỷ lệ cho nhóm tham chiếu (FT)
rate_ref_thap5 <- R_tablesize["M", "Rate_thap"]
rate_ref_trung5 <- R_tablesize["M", "Rate_trung"]
rate_ref_cao5 <- R_tablesize["M", "Rate_cao"]
# Tính Relative Risk cho các nhóm so với nhóm tham chiếu (M)
R_tablesize$RR_thap <- R_tablesize$Rate_thap / rate_ref_thap5
R_tablesize$RR_trung <- R_tablesize$Rate_trung / rate_ref_trung5
R_tablesize$RR_cao <- R_tablesize$Rate_cao / rate_ref_cao5
print(R_tablesize)
## thap trung cao Total Rate_thap Rate_trung Rate_cao RR_thap RR_trung
## L 329 454 257 1040 0.3163462 0.4365385 0.24711538 1.330949 0.7637926
## M 3638 8748 2920 15306 0.2376846 0.5715406 0.19077486 1.000000 1.0000000
## S 126 56 6 188 0.6702128 0.2978723 0.03191489 2.819757 0.5211744
## RR_cao
## L 1.2953247
## M 1.0000000
## S 0.1672909
Nhận xét tỷ lệ Relative Risk ở mức lương thấp với quy mô công ty
Kết quả Relative Risk là 2.8 điều này cho thấy xác suất để những kỹ sư làm việc ở những công ty có quy mô nhỏ (S) nhận thức lương “thấp” cao hơn xác suất kỹ sư làm việc ở các công ty quy mô trung bình đạt được mức lương thấp là khoảng 2.8
Kết quả của nhóm kỹ sư làm việc ở công ty có quy mô lớn (L) là 1.3 cho thấy những kỹ sư trong nhóm làm việc ở công ty có quy mô lớn có khả năng nhận được mức lương “thấp” với tỷ lệ cao gấp 1.3 lần so với nhóm kỹ sư làm việc ở công ty có quy mô trung bình.
Nhận xét tỷ lệ Relative Risk ở mức lương trung với quy mô công ty
Khi phân tích mức lương trung với quy mô công ty, kết quả cho thấy RR ở nhóm kỹ sư làm việc ở công ty có quy mô nhỏ (S) là 0.52. Kết quả, này cho thấy những kỹ sư trong nhóm công ty nhỏ có khả năng nhận mức lương “trung” với tỷ lệ chỉ bằng 0.52 lần so với nhóm kỹ sư làm việc ở công ty có quy mô trung (M).
Ngoài ra nhóm kỹ sư làm việc ở công ty có quy mô lớn (L) là 0.76 cho thấy những kỹ sư trong nhóm công ty lớn có khả năng nhận mức lương trung với tỷ lệ chỉ bằng 0.76 lần so với nhóm kỹ sư làm việc ở công ty có quy mô trung (M).
Nhận xét tỷ lệ Relative Risk ở mức lương cao với quy mô công ty
Kỹ sư làm việc ở công ty có quy mô nhỏ có khả năng nhận mức lương “cao” bằng 0.167 lần so với nhóm kỹ sư làm việc ở công ty có quy mô trung.
Trong nhóm làm việc ở công ty có quy mô lớn thì khả năng nhận mức lương “cao” bằng 1.29 lần so với nhóm kỹ sư làm việc ở công ty mô trung.
Trong phần này em sẽ thực hiện khoảng ước lượng cho RR với khoảng tin cậy 95%. Tuy nhiên để cho kết quả dễ dàng phân tích hơn trước nên em sẽ thực hiện trên bảng có số lượng nhỏ hơn là 2x2.
Trước tiên em tính khoảng ước lượng Relative Risk đối với mức lương thấp và trung nhưng trong nhóm công ty quy mô nhỏ và quy mô trung.
RR5 <- matrix(c(329, 454,
3638, 8748),
nrow = 2,
byrow = TRUE,
dimnames = list(c( "L", "M"),
c("thap", "trung")))
epitab(RR5, method = 'riskratio', rev = 'c')
## $tab
## trung p0 thap p1 riskratio lower upper p.value
## L 454 0.5798212 329 0.4201788 1.0000000 NA NA NA
## M 8748 0.7062813 3638 0.2937187 0.6990327 0.6409833 0.7623392 3.625974e-13
##
## $measure
## [1] "wald"
##
## $conf.level
## [1] 0.95
##
## $pvalue
## [1] "fisher.exact"
Xác suất để có mức lương thấp đối với nhóm làm việc ở công ty có quy mô trung à khoảng 29.37% và xác suất để đạt mức lương thấp đối với nhóm kỹ sư làm việc ở công ty lớn là 42.0%, trong trường hợp này tỷ lệ Relative Risk (RR) là khoảng 0.699 điều này cho thấy khả năng để nhóm kỹ sư đang làm việc ở công ty có quy mô trung đạt được mức lương “thấp” khoảng 0.699 lần đối với nhóm làm việc ở công ty có quy mô lớn. Khoảng ước lượng Relative Risk có kết quả là 0.64 và 0.76 đối với độ tin cậy 95%. Giá trị p.value khoảng 3.625974e-13 nhỏ hơn mức ý nghĩa 5% vì vậy kết quả này có ý nghĩa thống kê cho thấy có sự khác biệt đối với hai quy mô này.
RR6 <- matrix(c( 3638, 8748,
126, 56 ),
nrow = 2,
byrow = TRUE,
dimnames = list(c( "M", "S"),
c("thap", "trung")))
epitab(RR6, method = 'riskratio', rev = 'c')
## $tab
## trung p0 thap p1 riskratio lower upper p.value
## M 8748 0.7062813 3638 0.2937187 1.000000 NA NA NA
## S 56 0.3076923 126 0.6923077 2.357043 2.131395 2.60658 4.637891e-28
##
## $measure
## [1] "wald"
##
## $conf.level
## [1] 0.95
##
## $pvalue
## [1] "fisher.exact"
xác suất để đạt mức lương thấp đối với nhóm kỹ sư làm việc ở công ty nhỏ là 69.2%, trong trường hợp này tỷ lệ Relative Risk (RR) là khoảng 2.35 điều này cho thấy khả năng để nhóm kỹ sư đang làm việc ở công ty có quy mô nhỏ đạt được mức lương “thấp” khoảng 2.35 lần đối với nhóm làm việc ở công ty có quy mô trung Khoảng ước lượng Relative Risk có kết quả là 2.13 và 2.60 đối với độ tin cậy 95%. Giá trị p.value khoảng 4.637891e-28 nhỏ hơn mức ý nghĩa 5% vì vậy kết quả này có ý nghĩa thống kê cho thấy có sự khác biệt đối với hai quy mô này.
Trong phần này em sẽ tính tỷ lệ Odd Ratio của hai biến lương và quy mô công ty. Để dễ dàng trong phân tích em sẽ tách số liệu từ bảng tần số ban đầu thành bảng tần số mới chỉ bao gồm nhóm lương thấp, cao và quy mô công ty nhỏ, quy mô công ty trung.
odd5 <- table(d$size,d$sal)
odd6 <- odd5[c("M", "S"), c("thap", "cao")]
epitab(odd6, method = 'oddsratio')
## $tab
##
## thap p0 cao p1 oddsratio lower upper
## M 3638 0.96652497 2920 0.997949419 1.00000000 NA NA
## S 126 0.03347503 6 0.002050581 0.05932811 0.02611879 0.1347622
##
## p.value
## M NA
## S 8.479312e-25
##
## $measure
## [1] "wald"
##
## $conf.level
## [1] 0.95
##
## $pvalue
## [1] "fisher.exact"
Kết quả Odd Ratio của S so với M khoảng 0.05932811 điều này có nghĩa là odd mức lương cao của nhóm kỹ sư làm việc ở công ty quy mô nhỏ chỉ bằng khoảng 5,9% so với odd kỹ sư làm việc ở công ty có quy mô trung.
Khoảng tin cậy của OR đối với độ tin cậy 95% là từ 0.0261 đến 0.1348 và giá trị p là 8.479312e-25, nhỏ hơn 0.05 cho thấy có sự khác biệt giữa odd ở hai nhóm S và M
Trái ngược lại với phần trên chúng ta có thể sử dụng tham số rev = ‘c’ thay đổi giữa mức lương cao và mức lương thấp để tính odd ratio của nhóm làm việc ở công ty có quy mô nhỏ với nhóm làm việc ở công ty có quy mô trung ở mức lương thấp.
epitab(odd6, method = 'oddsratio', rev = 'c' )
## $tab
##
## cao p0 thap p1 oddsratio lower upper p.value
## M 2920 0.997949419 3638 0.96652497 1.00000 NA NA NA
## S 6 0.002050581 126 0.03347503 16.85542 7.42048 38.28661 8.479312e-25
##
## $measure
## [1] "wald"
##
## $conf.level
## [1] 0.95
##
## $pvalue
## [1] "fisher.exact"
Kết quả cho thấy Odd ratio của S so với M là 16.85542 ta có thể giải thích rằng odd của nhóm kỹ sư làm việc ở công ty có quy mô nhỏ có mức lương “thấp” cao hơn gấp 16.85 lần so với odd của nhóm làm việc tại công ty quy mô trung.
Kết quả chỉ ra rằng giá trị p là 8.479312e-25, nhỏ hơn 0.05, do đó có sự khác biệt giữa odd của hai nhóm trên. Ngoài ra, khoảng ước lượng cho odd ratio đối với độ tin cậy 95% là 7.42 và 38.29.
Phân tích mức lương trung giữa nhóm kỹ sư làm việc tại quy mô nhỏ và trung
odd7 <- odd5[c("M", "S"), c("thap", "trung")]
epitab(odd7, method = 'oddsratio')
## $tab
##
## thap p0 trung p1 oddsratio lower upper p.value
## M 3638 0.96652497 8748 0.993639255 1.0000000 NA NA NA
## S 126 0.03347503 56 0.006360745 0.1848295 0.1345977 0.253808 4.637891e-28
##
## $measure
## [1] "wald"
##
## $conf.level
## [1] 0.95
##
## $pvalue
## [1] "fisher.exact"
Kết quả Odds ratio là 0.185, điều này có nghĩa là odds của kỹ sư có mức lương “trung” trong nhóm làm việc tại công ty nhỏ chỉ bằng khoảng 18.5% so với odds nhóm làm việc toàn thời gian có cùng mức lương này.
Ngoài ra khoảng tin cậy của OR là từ 0.1345 đến 0.2538 đối với khoảng tin cậy 95% và giá trị p là 4.637891e-28, nhỏ hơn 0.05, cho thấy sự khác biệt giữa odds của nhóm S và M.
epitab(odd7, method = 'oddsratio', rev = 'c')
## $tab
##
## trung p0 thap p1 oddsratio lower upper p.value
## M 8748 0.993639255 3638 0.96652497 1.00000 NA NA NA
## S 56 0.006360745 126 0.03347503 5.41039 3.939986 7.429549 4.637891e-28
##
## $measure
## [1] "wald"
##
## $conf.level
## [1] 0.95
##
## $pvalue
## [1] "fisher.exact"
Kết quả cho thấy odds kỹ sư có mức lương thấp đang làm việc ở công ty có quy mô nhỏ gấp khoảng 5.41 lần so với odds của kỹ sư đang làm việc ở công ty có quy mô trung.
Hơn thế kết quả cũng cho rằng khoảng tin cậy của OR là từ 3.94 đến 7.43 và giá trị p là 4.637891e-28, nhỏ hơn 0.05, cho thấy sự khác biệt giữa odds của hai nhóm là có ý nghĩa thống kê.
Kết quả bên dưới được sử dụng để kiểm tra tính độc lập giữa hai biến quy mô công ty và mức lương bằng phương pháp chi bình phương. Giả thuyết Ho được đưa ra là biến mức lương và quy mô công ty độc lập.
chisq.test(table(d$size,d$sal))
##
## Pearson's Chi-squared test
##
## data: table(d$size, d$sal)
## X-squared = 258.13, df = 4, p-value < 2.2e-16
Qua kết quả kiểm định cho ta thấy giá trị p−value là 2.2e−16 nhỏ hơn 0.05 , nên bác bỏ H0, nghĩa là mức lương và quy mô công ty là có liên quan với nhau.
Ước lượng tỷ lệ cho nhóm lương trung bình
Trước tiên ta cần lọc kết quả về mức lương theo nhóm làm việc ở công ty có quy trung bình. Kết quả thu được như sau.
M <- table(d[d$size == "M", ]$sal)
M
##
## cao trung thap
## 2920 8748 3638
prop.test(M["trung"], sum(M), p = 0.7)
##
## 1-sample proportions test with continuity correction
##
## data: M["trung"] out of sum(M), null probability 0.7
## X-squared = 1202.1, df = 1, p-value < 2.2e-16
## alternative hypothesis: true p is not equal to 0.7
## 95 percent confidence interval:
## 0.5636512 0.5793939
## sample estimates:
## p
## 0.5715406
Sau đó ước lượng tỷ lệ với giả thuyết không (H0): Tỷ lệ người có mức lương trung và làm việc tại công ty quy mô vừa là 0.7 (70%).
Kết quả từ kiểm định cho thấy giá trị p-value là 2.2e-16, rất gần với 0 và nhỏ hơn giá trị α=0.05. Điều này cho thấy có bằng chứng mạnh mẽ để bác bỏ giả thuyết rằng tỷ lệ người có mức lương trung trong nhóm toàn thời gian là 70%.
Kết quả cũng cho thấy tỷ lệ kỹ sư có mức lương trung trong công ty có quy mô trung là 57.15%. Khoảng ước lượng tỷ lệ ở độ tin cây 95% nằm trong khoảng từ 55.37% đến 57.94%.
Ước lượng tỷ lệ cho nhóm lương thấp
prop.test(M["thap"], sum(M), p = 0.2)
##
## 1-sample proportions test with continuity correction
##
## data: M["thap"] out of sum(M), null probability 0.2
## X-squared = 135.62, df = 1, p-value < 2.2e-16
## alternative hypothesis: true p is not equal to 0.2
## 95 percent confidence interval:
## 0.2309751 0.2445263
## sample estimates:
## p
## 0.2376846
Ta có giả thuyết không (H0): Tỷ lệ người có mức lương thấp trong nhóm làm việc ở công ty có quy mô trung là 0.2 (20%). Vì giá trị p rất nhỏ (p < 2.2e-16), chúng ta bác bỏ giả thuyết không (H0). Điều này có nghĩa là tỷ lệ kỹ sư có mức lương “thấp” trong nhóm làm việc ở công ty trung không bằng 20%. Kết quả cho rằng tỷ lệ này đạt khoảng 23.77% và ở mức độ tin cậy 95% kết quả ước lượng chỉ ra tỷ lệ kỹ sư có mức lương thấp trong nhóm toàn thời gian đạt khoảng từ 23.1% đến 24.45%.
prop.test(M["cao"], sum(M), p = 0.3)
##
## 1-sample proportions test with continuity correction
##
## data: M["cao"] out of sum(M), null probability 0.3
## X-squared = 869.02, df = 1, p-value < 2.2e-16
## alternative hypothesis: true p is not equal to 0.3
## 95 percent confidence interval:
## 0.1845959 0.1971098
## sample estimates:
## p
## 0.1907749
Ta có giả thuyết không (H0) là tỷ lệ kỹ sư có mức lương “cao” trong nhóm công ty có quy mô trung là 30%. Kết quả bên trên có giá trị p rất nhỏ (p < 2.2e-16), chúng ta bác bỏ giả thuyết không (H0). Điều này có nghĩa là tỷ lệ kỹ sư có mức lương “cao” trong nhóm nhóm làm công ty có quy mô trung sẽ không bằng 30%. Thực tế, tỷ lệ kỹ sư có mức lương “cao” chỉ khoảng 19%, thấp hơn so với tỷ lệ kỳ vọng là 30%. Ngoài ra, khoảng tin cậy 95% cho tỷ lệ kỹ sư có mức lương cao trong nhóm toàn thời gian là [0.1846, 0.1971]. Điều này có nghĩa là chúng ta có 95% tin tưởng rằng tỷ lệ này nằm trong khoảng từ 18.46% đến 19.71%.
Mô hình xác suất tuyến tính
Phần này sẽ thực hiện hồi quy xác suất tuyến tính giữa biến độc lập là quy mô công ty và biến phụ thuộc là lương. Tuy nhiên trong phần này ta cần tiếp tục thay đổi biến sal thành biến nhị phân với giá trị 1 tương ứng với mức lương trung và cao, ngược lại giá trị 0 tương ứng với mức lương thấp.
d$sal3 <- ifelse(d$sal %in% c("trung", "cao"), 1, 0)
d$size <- relevel(as.factor(d$size), ref = "M")
glm4 <- glm(sal3 ~ size, data = d)
summary(glm4)
##
## Call:
## glm(formula = sal3 ~ size, data = d)
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.762315 0.003466 219.935 < 2e-16 ***
## sizeL -0.078662 0.013741 -5.724 1.06e-08 ***
## sizeS -0.432528 0.031466 -13.746 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 0.1838835)
##
## Null deviance: 3079.8 on 16533 degrees of freedom
## Residual deviance: 3039.8 on 16531 degrees of freedom
## AIC: 18927
##
## Number of Fisher Scoring iterations: 2
\(\hat{\pi} = 0.762315 - 0.078662 \cdot sizeL - 0.432528 \cdot sizeS\)
Tất cả các hệ số đều có giá trị \(p\) nhỏ hơn mức ý nghĩa \(\alpha = 0.05\), ngụ ý rằng tất cả các hệ số đều có ảnh hưởng ý nghĩa thống kê đến mức lương.
Hệ số \(\hat{\beta}_0 = 0.765156\) có ý nghĩa rằng tỷ lệ để người kỹ sư đạt được mức lương trung bình là 76.51% trong trường hợp không xuất hiện các hình thức làm việc khác.
Đối với hệ số \(\hat{\beta}_1 = - 0.078662\), ta nhận thấy khi kỹ sư làm việc ở công ty có quy mô lớn khi không bị ảnh hưởng bởi các yếu tố khác thì chênh lệch tỷ lệ đạt mức lương trung trở lên giữa nhóm làm việc ở công ty có quy mô lớn so với làm việc ở công ty có quy mô vừa là -7.8%, hay có thể nói cách khác tỷ lệ bình quân để đạt mức lương trung trở lên đối với kỹ sư làm việc tại công ty vừa là 68.95%.
Tương tự, \(\hat{\beta}_2 = -0.432528\), ta nhận thấy khi kỹ sư làm việc ở công ty có quy mô nhỏ và không bị ảnh hưởng bởi các yếu tố khác thì chênh lệch tỷ lệ đạt mức lương trung trở lên giữa nhóm làm việc ở công ty nhỏ và nhóm là việc ở công ty trung là -43.25%, hay tỷ lệ đạt mức lương trung trở đối với kỹ sư làm việc tại công ty nhỏ là 33.56%.
Ta nhận thấy tác động ở cả ba nhóm quy mô công ty nhỏ, trung và lớn. Nhưng nếu không có sự tác động của ba nhóm trên thì xác suất để đạt mức lương trên là khá cao, khoảng 0.762315.
nsize <- data.frame(size = c("S", "M", "L") )
predict.glm(glm4, newdata = nsize, type="response")
## 1 2 3
## 0.3297872 0.7623154 0.6836538
Dựa vào kết quả ta có thể nhận thấy xác suất dự báo cho biểu hiện kỹ sư làm việc tại công ty có quy mô nhỏ (1) là khoảng 33%. Điều này có nghĩa một kỹ sư làm việc ở công ty có quy mô nhỏ có xác suất nhận mức lương trung đến cao là khoảng 33% .
Bên cạnh đó kết quả dự báo đối với nhóm kỹ sư làm việc ở công ty có quy mô trung (2) thì xác suất để đạt được mức lương trung trở lên đạt 76.23%, vì vậy nhóm này có xác suất rất cao để nhận được mức lương trung trở lên.
Trong nhóm kỹ sư làm việc ở công ty có quy mô lớn (3) thì xác suất để đạt mức lương trung khoảng 68.37%, cao hơn so với xác suất của nhóm làm ở quy mô nhỏ.
Mô hình Logit
logitsize <- glm(factor(sal2) ~ size, data = d, family = binomial(link = 'logit'))
summary(logitsize)
##
## Call:
## glm(formula = factor(sal2) ~ size, family = binomial(link = "logit"),
## data = d)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.16542 0.01899 61.373 < 2e-16 ***
## sizeL -0.39480 0.06933 -5.695 1.24e-08 ***
## sizeS -1.87456 0.15629 -11.994 < 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: 18506 on 16533 degrees of freedom
## Residual deviance: 18324 on 16531 degrees of freedom
## AIC: 18330
##
## Number of Fisher Scoring iterations: 4
Ta có kết quả từ hàm hồi quy thông qua mô hình logit như sau:
\[ \ln\left(\frac{\hat{\pi}}{1 - \hat{\pi}}\right) = 1.16542 - 0.39480 \cdot sizeL - 0.08898 \cdot sizeS \]
Hệ số chặn là 1.16542 điều này có nghĩa là khi không có sự tác động của các quy mô công ty khác thì \(\ln\left(\frac{\hat{\pi}}{1 - \hat{\pi}}\right)\) của kỹ sư làm việc trực tiếp đạt mức lương trung là 1.16542.
Đối với hệ số \(\hat{\beta}_1 = -0.39480\), điều này cho thấy chênh lệch \(\ln\left(\frac{\hat{\pi}}{1 - \hat{\pi}}\right)\) đạt mức lương trung trở lên giữa nhóm vừa làm viêc ở công ty có quy mô lớn so với nhóm làm việc ở công ty quy mô trung là -0.39480.
Đối với hệ số \(\hat{\beta}_2 = -0.08898\), điều này cho thấy chênh lệch \(\ln\left(\frac{\hat{\pi}}{1 - \hat{\pi}}\right)\) đạt mức lương trung trở lên giữa nhóm làm việc ở công ty có quy mô nhỏ so với làm việc ở công ty có quy mô trung là -0.08898
Sau đó từ mô hình logit ta dự báo xác suất của biến lương tương ứng với mỗi quy mô công ty. Ta cần một khung dữ liệu n với một cột size chứa các giá trị “S”, “M”, “L”.
size <- data.frame(size = c("S", "M", "L") )
predict.glm(logitsize, newdata = size, type="response")
## 1 2 3
## 0.3297872 0.7623154 0.6836538
Dựa vào kết quả ta có thể nhận thấy xác suất dự báo cho biểu hiện kỹ sư làm việc tại công ty có quy mô nhỏ (1) là khoảng 33%. Điều này có nghĩa một kỹ sư làm việc ở công ty có quy mô nhỏ có xác suất nhận mức lương trung đến cao là khoảng 33% .
Bên cạnh đó kết quả dự báo đối với nhóm kỹ sư làm việc ở công ty có quy mô trung (2) thì xác suất để đạt được mức lương trung trở lên đạt 76.23%, vì rất cao cho nên đối với nhóm này có xác suất rất cao để nhận được mức lương trung trở lên.
Trong nhóm kỹ sư làm việc ở công ty có quy mô lớn (3) thì xác suất để đạt mức lương trung khoảng 68.37%, cao hơn so với xác suất của nhóm làm ở quy mô nhỏ.
Mô hình Probit
Tiếp tục, thay vì sử dụng mô hình logit, ta có thể thay thế bằng mô hình probit dự đoán xác suất của một sự kiện xảy ra. Mô hình probit sử dụng hàm phân phối chuẩn. Trong đó hàm liên kết có dạng \(\text{probit}(\pi) = \Phi^{-1}(\pi)\) và hàm xác suất là \(\pi(x) = \Phi(\beta_0 + \beta_1 X + \beta_2 X_2 + \ldots + \beta_n X_n)\).
probitsize <- glm(factor(sal2) ~ size, data = d, family = binomial(link = 'probit'))
summary(probitsize)
##
## Call:
## glm(formula = factor(sal2) ~ size, family = binomial(link = "probit"),
## data = d)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.71377 0.01113 64.151 <2e-16 ***
## sizeL -0.23583 0.04202 -5.612 2e-08 ***
## sizeS -1.15427 0.09536 -12.105 <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: 18506 on 16533 degrees of freedom
## Residual deviance: 18324 on 16531 degrees of freedom
## AIC: 18330
##
## Number of Fisher Scoring iterations: 4
Ta có kết quả mô hình hồi quy như sau:
\[ \hat{\pi} = \Phi(0.71377 - 1.15427 \times sizeS - 0.23583 \times sizeL) \]
Hệ số chặn là 0.71377 và có ý nghĩa thống kê khi giá trị p_value nhỏ hơn mức ý nghĩa 5%. Điều này có nghĩa là khi không có sự tác động các biến độc lập khác thì xác suất để kỹ sư đạt mức lương trung là \(\Phi(0.71377)\).
Hệ số \(\hat{\beta}_1 = -1.15427\) có ý nghĩa là trong trường hợp các yếu tố khác không đổi thì chênh lệch tỷ lệ đạt mức lương trung giữa nhóm làm việc ở công ty có quy mô nhỏ so với nhóm làm việc ở công ty có quy mô trung là \(\Phi(-1.15427)\).
Hệ số \(\hat{\beta}_2 = -0.23583\) có ý nghĩa là trong trường hợp các yếu tố khác không đổi thì chênh lệch tỷ lệ đạt mức lương trung giữa nhóm làm ở công ty có quy mô lớn so với nhóm làm việc ở công ty có quy mô trung là \(\Phi(-0.23583)\).
Từ kết quả trên, em sẽ thực hiện dự báo cho tỷ lệ để đạt được mức lương trung ở các quy mô công ty. Kết quả được trình bày trong bảng bên dưới:
predict.glm(probitsize, newdata = size, type="response")
## 1 2 3
## 0.3297872 0.7623154 0.6836538
Dựa vào kết quả ta có thể nhận thấy xác suất dự báo cho biểu hiện kỹ sư làm việc tại công ty có quy mô nhỏ (1) là khoảng 33%. Điều này có nghĩa một kỹ sư làm việc ở công ty có quy mô nhỏ có xác suất nhận mức lương trung đến cao là khoảng 33% .
Bên cạnh đó kết quả dự báo đối với nhóm kỹ sư làm việc ở công ty có quy mô trung (2) thì xác suất để đạt được mức lương trung trở lên đạt 76.23%, vì tỉ lệ rất cao cho nên đối với nhóm này có xác suất rất cao để nhận được mức lương trung trở lên.
Trong nhóm kỹ sư làm việc ở công ty có quy mô lớn (3) thì xác suất để đạt mức lương trung khoảng 68.37%, cao hơn so với xác suất của nhóm làm ở quy mô nhỏ.
Chỉ số AIC
AIC(glm4)
## [1] 18926.91
AIC(logitsize)
## [1] 18330
AIC(probitsize)
## [1] 18330
Mô hình Logit và Probit đều tốt hơn mô hình xác suất tuyến tính trong việc giải thích mối quan hệ giữa quy mô công ty và xác suất đạt mức lương từ trung trở lên, Vì chỉ số AIC của hai mô hình này thấp hơn chỉ số AIC của mô hình xác suất tuyến tính. Logit và Probit có hiệu suất tương đương nhau vì chúng có cùng chỉ số AIC
Hệ số Brier Score
Phần này sẽ tính chỉ số Brier Score được sử dụng để đo lường độ chính xác của các dự đoán xác suất. Nó tính toán độ lệch bình phương giữa dự đoán xác suất và giá trị thực tế. Một giá trị Brier Score thấp cho thấy mô hình dự đoán tốt hơn. Vì thế ta sẽ tính và so sánh giá trì này ở hai mô hình hồi quy logit và probit.
BrierScore(logitsize)
## [1] 0.1838502
BrierScore(probitsize)
## [1] 0.1838502
Kết quả bên trên cho thấy Brier Score cho logit và probit đều là 0.1838502 Điều này có nghĩa là cả hai mô hình đều có hiệu suất tương đương khi dự đoán xác suất cho biến mục tiêu. Brier Score khá thấp (gần 0), cho thấy rằng cả hai mô hình đều dự đoán khá tốt và không sai lệch nhiều so với giá trị thực tế.
Ma trận nhầm lẫn
Conf(table(predict(probitsize, size="response") >=0.5,probitsize$data$sal2 == '1'))
##
## Confusion Matrix and Statistics
##
## Reference
## Prediction FALSE TRUE
## FALSE 455 773
## TRUE 3638 11668
##
## Total n : 16'534
## Accuracy : 0.7332
## 95% CI : (0.7264, 0.7399)
## No Information Rate : 0.7524
## P-Value [Acc > NIR] : 1.0000
##
## Kappa : 0.0641
## Mcnemar's Test P-Value : < 2.2e-16
##
## Sensitivity : 0.1112
## Specificity : 0.9379
## Pos Pred Value : 0.3705
## Neg Pred Value : 0.7623
## Prevalence : 0.2476
## Detection Rate : 0.0743
## Detection Prevalence : 0.0275
## Balanced Accuracy : 0.5245
## F-val Accuracy : 0.1710
## Matthews Cor.-Coef : 0.0807
##
## 'Positive' Class : FALSE
Từ kết quả trên ma trận cho ta biết có 455 số lần mô hình dự đoán đúng là mức lương thấp. Tuy nhiên có 773 lần mô hình dự đoán sai đối với mức lương trung và khoảng 3638 lần mô hình dự đoán sai mức lương thấp. Số lần mô hình dự đoán đúng ở mức lương trung là 11668 lần đạt cao nhất.
Ngoài ra ta có tỷ lệ dự đoán đúng của mô hình là 0.7332 và ở độ tin cậy là 95% thì tỷ lệ này nằm trong khoảng từ (0.7264, 0.7264)
Dưới đây em sẽ thực hiện hồi quy xác suất tuyến tính trong đó các biến độc lập là loại hình công việc (type), hình thức làm việc (remote), trình độ (remote), quy mô doanh nghiệp (company_size) sẽ tác động đến mức lương của những người kỹ sư. Tương tự với các phần trước đó em sẽ xem xét xác suất để xảy ra mức lương từ trung trở lên.
glm4 <- glm(sal2 ~ type + remote + el + size + job_group,
data = d,
family = binomial(link = "logit"))
summary(glm4)
##
## Call:
## glm(formula = sal2 ~ type + remote + el + size + job_group, family = binomial(link = "logit"),
## data = d)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.85669 0.04592 40.432 < 2e-16 ***
## typeCT -0.62447 0.45670 -1.367 0.17152
## typeFL -14.01992 122.87171 -0.114 0.90916
## typePT -0.55967 0.39618 -1.413 0.15775
## remote50 -1.86181 0.17640 -10.554 < 2e-16 ***
## remote100 -0.14113 0.04421 -3.192 0.00141 **
## elEN -2.09728 0.06901 -30.392 < 2e-16 ***
## elEX 0.92959 0.17524 5.305 1.13e-07 ***
## elMI -1.31945 0.04439 -29.725 < 2e-16 ***
## sizeL -0.17302 0.08823 -1.961 0.04988 *
## sizeS -1.37084 0.18752 -7.310 2.67e-13 ***
## job_groupDS 0.38443 0.05702 6.742 1.56e-11 ***
## job_groupML 1.16266 0.08265 14.067 < 2e-16 ***
## job_groupOT -0.75393 0.05067 -14.879 < 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: 18506 on 16533 degrees of freedom
## Residual deviance: 15020 on 16520 degrees of freedom
## AIC: 15048
##
## Number of Fisher Scoring iterations: 12
Từ kết quả mô hình hồi quy logistic với biến phụ thuộc là xác suất kỹ sư có mức lương từ trung bình trở lên, ta có phương trình ước lượng như sau:
\[ \hat{\pi} = 1.8567 - 0.6245*CT - 14.0199*FL - 0.5597*PT - 1.8618*remote50 - 0.1411*remote100 - 2.0973*EN + 0.9296*EX - 1.3195*MI - 0.1730*sizeL - 1.3708*sizeS + 0.3844*job\_groupDS + 1.1627*job\_groupML - 0.7539*job\_groupOT \]
Dựa trên kết quả mô hình hồi quy logistic đã bao gồm các yếu tố: hình thức làm việc, tỷ lệ làm việc từ xa, cấp độ kinh nghiệm, quy mô công ty và nhóm công việc, em tiến hành tạo bảng dữ liệu giả định để dự báo xác suất kỹ sư đạt mức lương từ trung bình trở lên.
# Tạo bộ dữ liệu giả định để dự báo
new_data <- data.frame(
type = c("FT", "FT", "FT", "FT"),
remote = factor(c("0", "0", "0", "0")),
el = factor(c("SE", "SE", "MI", "MI")),
size = factor(c("L", "S", "L", "S")),
job_group = factor(c("DE", "DE", "ML", "ML"))
)
new_data
## type remote el size job_group
## 1 FT 0 SE L DE
## 2 FT 0 SE S DE
## 3 FT 0 MI L ML
## 4 FT 0 MI S ML
# Tính xác suất dự báo
predictions <- predict(glm4, new_data, type = "response")
predictions
## 1 2 3 4
## 0.8433905 0.6191275 0.8215504 0.5815303
Trường hợp kỹ sư làm việc toàn thời gian, trực tiếp tại công ty, có chuyên môn cao cấp (SE) và làm việc tại công ty quy mô lớn có xác suất đạt mức lương từ trung trở lên là 84.34%. Đây là nhóm có khả năng cao nhất để đạt mức lương trung bình hoặc cao.
Khi giữ nguyên các điều kiện nhưng làm việc tại công ty quy mô nhỏ, xác suất giảm còn 61.91%, cho thấy quy mô doanh nghiệp có ảnh hưởng rõ rệt đến mức lương.
Đối với kỹ sư trung cấp (MI), làm việc toàn thời gian, trực tiếp tại công ty trong công ty quy mô lớn, xác suất để đạt mức lương trung trở lên là 82.16%. Điều này thể hiện tầm quan trọng của nhóm công việc và quy mô doanh nghiệp, ngay cả khi chuyên môn chưa ở mức cao cấp.
Cuối cùng, trường hợp kỹ sư trung cấp (MI) làm việc tại công ty quy mô nhỏ chỉ có xác suất 58.15%, đây là nhóm có khả năng đạt mức lương trung trở lên thấp nhất trong các trường hợp xét đến.
log4 <- glm(factor(sal2) ~ type + remote + el + size + job_group,
data = d,
family = binomial(link = 'logit'))
summary(log4)
##
## Call:
## glm(formula = factor(sal2) ~ type + remote + el + size + job_group,
## family = binomial(link = "logit"), data = d)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.85669 0.04592 40.432 < 2e-16 ***
## typeCT -0.62447 0.45670 -1.367 0.17152
## typeFL -14.01992 122.87171 -0.114 0.90916
## typePT -0.55967 0.39618 -1.413 0.15775
## remote50 -1.86181 0.17640 -10.554 < 2e-16 ***
## remote100 -0.14113 0.04421 -3.192 0.00141 **
## elEN -2.09728 0.06901 -30.392 < 2e-16 ***
## elEX 0.92959 0.17524 5.305 1.13e-07 ***
## elMI -1.31945 0.04439 -29.725 < 2e-16 ***
## sizeL -0.17302 0.08823 -1.961 0.04988 *
## sizeS -1.37084 0.18752 -7.310 2.67e-13 ***
## job_groupDS 0.38443 0.05702 6.742 1.56e-11 ***
## job_groupML 1.16266 0.08265 14.067 < 2e-16 ***
## job_groupOT -0.75393 0.05067 -14.879 < 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: 18506 on 16533 degrees of freedom
## Residual deviance: 15020 on 16520 degrees of freedom
## AIC: 15048
##
## Number of Fisher Scoring iterations: 12
Kết quả từ mô hình logit như sau:
\[ \log(\text{odds}) = 1.8567 - 0.6245\cdot CT - 14.0199\cdot FL - 0.5597\cdot PT - 1.8618\cdot remote50 - 0.1411\cdot remote100 - 2.0973\cdot EN + 0.9296\cdot EX - 1.3195\cdot MI - 0.1730\cdot sizeL - 1.3708\cdot sizeS + 0.3844\cdot DS + 1.1627\cdot ML - 0.7539\cdot OT \]
Hệ số \(\hat{\beta}_0 = 1.8567\): Khi các biến độc lập ở nhóm tham chiếu (làm toàn thời gian, không làm từ xa, kỹ sư cao cấp, công ty quy mô vừa, nhóm Data Engineering), log(odds) để có lương từ trung trở lên là 1.8567, tương ứng odds ≈ 6.40.
\(\hat{\beta}_{CT} = -0.6245\): Chênh lệch log(odds) giữa kỹ sư làm hợp đồng và toàn thời gian là -0.6245. Odds của nhóm hợp đồng chỉ bằng \(e^{-0.6245} \approx 0.535\) lần so với toàn thời gian.
\(\hat{\beta}_{FL} = -14.0199\): Nhóm làm tự do có odds gần như bằng 0 so với nhóm toàn thời gian, cho thấy khả năng đạt lương trung trở lên cực thấp.
\(\hat{\beta}_{PT} = -0.5597\): Kỹ sư bán thời gian có odds bằng \(e^{-0.5597} \approx 0.571\) lần so với toàn thời gian.
\(\hat{\beta}_{remote50} = -1.8618\): Nhóm làm việc kết hợp (50% từ xa) có odds chỉ bằng khoảng \(e^{-1.8618} \approx 0.155\) lần so với làm hoàn toàn tại công ty.
\(\hat{\beta}_{remote100} = -0.1411\): Làm hoàn toàn từ xa có odds đạt lương trung trở lên bằng \(e^{-0.1411} \approx 0.868\) lần so với làm trực tiếp.
\(\hat{\beta}_{EN} = -2.0973\): Kỹ sư mới vào nghề có odds bằng khoảng \(e^{-2.0973} \approx 0.122\) lần so với kỹ sư cao cấp.
\(\hat{\beta}_{EX} = 0.9296\): Kỹ sư điều hành có odds bằng \(e^{0.9296} \approx 2.53\) lần so với kỹ sư cao cấp.
\(\hat{\beta}_{MI} = -1.3195\): Kỹ sư trung cấp có odds bằng \(e^{-1.3195} \approx 0.267\) lần so với kỹ sư cao cấp.
\(\hat{\beta}_{sizeL} = -0.1730\): Kỹ sư làm ở công ty quy mô lớn có odds ≈ \(e^{-0.1730} \approx 0.841\) lần so với công ty vừa.
\(\hat{\beta}_{sizeS} = -1.3708\): Kỹ sư làm ở công ty quy mô nhỏ có odds ≈ \(e^{-1.3708} \approx 0.254\) lần so với công ty vừa.
\(\hat{\beta}_{DS} = 0.3844\): Nhóm Data Science có odds ≈ \(e^{0.3844} \approx 1.47\) lần so với Data Engineering.
\(\hat{\beta}_{ML} = 1.1627\): Nhóm Machine Learning/AI có odds ≈ \(e^{1.1627} \approx 3.20\) lần so với Data Engineering.
\(\hat{\beta}_{OT} = -0.7539\): Các công việc khác có odds ≈ \(e^{-0.7539} \approx 0.470\) lần so với Data Engineering.
Tương tự với bảng dữ liệu mới được tạo trong phần xác suất tuyến tính em sẽ dùng để thực hiện dự báo trong trường hợp mô hình logit.
predictions <- predict(log4, new_data, type="response")
predictions
## 1 2 3 4
## 0.8433905 0.6191275 0.8215504 0.5815303
Kết quả đầu tiên cho thấy những kỹ sư làm việc toàn thời gian, làm trực tiếp tại công ty, đạt chuyên môn cao cấp và công tác tại công ty quy mô lớn thuộc nhóm Data Engineering có xác suất mức lương từ trung trở lên là 84.34%.
Thứ hai, với những kỹ sư làm việc toàn thời gian, làm trực tiếp tại công ty, đạt chuyên môn cao cấp và công tác tại công ty quy mô nhỏ trong nhóm Data Engineering, xác suất này giảm xuống còn 61.91%.
Thứ ba, các kỹ sư làm việc toàn thời gian, làm trực tiếp tại công ty, đạt chuyên môn trung cấp và công tác tại công ty quy mô lớn trong nhóm Machine Learning/AI có xác suất mức lương từ trung trở lên là 82.16%.
Cuối cùng, những kỹ sư làm việc toàn thời gian, làm trực tiếp tại công ty, đạt chuyên môn trung cấp và công tác tại công ty quy mô nhỏ trong nhóm Machine Learning/AI có xác suất mức lương từ trung trở lên chỉ là 58.15%.
Ngoài ra em cũng thực hiện dự báo xác suất cho mức lương trung trở lên thông qua mô hình probit. Tương đồng với việc thực hiện mô hình này cho các phần trên, tuy nhiên tại đây em sẽ thực hiện hồi quy loại hình công việc, hình thức làm việc, trình độ, quy mô doanh nghiệp và các nhóm công việc tác động đến mức lương của kỹ sư.
probit4 <- glm(factor(sal2) ~ type + remote + el + size + job_group,
data = d, family = binomial(link = 'probit'))
summary(probit4)
##
## Call:
## glm(formula = factor(sal2) ~ type + remote + el + size + job_group,
## family = binomial(link = "probit"), data = d)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.10078 0.02563 42.944 < 2e-16 ***
## typeCT -0.30176 0.27064 -1.115 0.26485
## typeFL -5.77312 32.74833 -0.176 0.86007
## typePT -0.24066 0.23309 -1.033 0.30184
## remote50 -1.05155 0.10092 -10.419 < 2e-16 ***
## remote100 -0.06973 0.02530 -2.756 0.00585 **
## elEN -1.25241 0.04101 -30.542 < 2e-16 ***
## elEX 0.49209 0.08790 5.598 2.17e-08 ***
## elMI -0.77162 0.02589 -29.804 < 2e-16 ***
## sizeL -0.09543 0.05057 -1.887 0.05917 .
## sizeS -0.75460 0.10928 -6.906 5.00e-12 ***
## job_groupDS 0.21743 0.03205 6.785 1.16e-11 ***
## job_groupML 0.62969 0.04334 14.529 < 2e-16 ***
## job_groupOT -0.44813 0.02960 -15.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: 18506 on 16533 degrees of freedom
## Residual deviance: 15017 on 16520 degrees of freedom
## AIC: 15045
##
## Number of Fisher Scoring iterations: 12
Kết quả từ mô hình Probit được trình bày:
\[ \hat\pi = \Phi(1.1008 - 0.3018 \cdot CT - 5.7731 \cdot FL - 0.2407 \cdot PT - 1.0516 \cdot remote50 - 0.0697 \cdot remote100 - 1.2524 \cdot EN + 0.4921 \cdot EX - 0.7716 \cdot MI - 0.0954 \cdot sizeL - 0.7546 \cdot sizeS + 0.2174 \cdot job\_groupDS + 0.6297 \cdot job\_groupML - 0.4481 \cdot job\_groupOT) \]
Kỹ sư điển hình (nhóm tham chiếu)
Một kỹ sư làm full-time, tại công ty, thuộc nhóm Senior (SE), làm ở công ty quy mô vừa, và thuộc ngành Data Engineering, có:
1.1008Φ(1.1008) ≈ 86.4%Yếu tố làm giảm xác suất lương cao
knitr::kable(
data.frame(
`Yếu tố thay đổi` = c(
"type = CT (Contract)", "type = FL (Freelance)", "type = PT (Part-time)",
"remote = 50%", "remote = 100%",
"el = EN (Entry-level)", "el = MI (Mid-level)",
"size = S (Small company)",
"job_group = OT (Other)"
),
`Tác động đến z-score` = c(
"-0.30", "-5.77", "-0.24",
"-1.05", "-0.07",
"-1.25", "-0.77",
"-0.75",
"-0.45"
),
`Ảnh hưởng đến xác suất (ước lượng)` = c(
"86.4% → ~78.7%", "≈ 0%", "≈ 80.1%",
"86.4% → ~65.6%", "86.4% → ~83.6%",
"86.4% → ~59.1%", "86.4% → ~71.7%",
"86.4% → ~72.0%",
"86.4% → ~76.1%"
),
`Giải thích dễ hiểu` = c(
"Giảm ~7.7%, do làm hợp đồng nên ít ổn định hơn",
"Rất thấp, freelancer gần như không đạt mức lương trung trở lên",
"Giảm nhẹ so với full-time",
"Giảm mạnh ~21%, do thiếu sự hiện diện trực tiếp",
"Ảnh hưởng rất nhẹ, remote hoàn toàn ổn định hơn hybrid",
"Mới ra trường, lương thấp hơn rõ",
"Thấp hơn Senior nhưng vẫn ổn",
"Công ty nhỏ trả lương thấp hơn rõ rệt",
"Ngành không chuyên sâu về dữ liệu thì lương thấp hơn"
)
),
caption = "Các yếu tố làm giảm xác suất đạt mức lương cao"
)
| Yếu.tố.thay.đổi | Tác.động.đến.z.score | Ảnh.hưởng.đến.xác.suất..ước.lượng. | Giải.thích.dễ.hiểu |
|---|---|---|---|
| type = CT (Contract) | -0.30 | 86.4% → ~78.7% | Giảm ~7.7%, do làm hợp đồng nên ít ổn định hơn |
| type = FL (Freelance) | -5.77 | ≈ 0% | Rất thấp, freelancer gần như không đạt mức lương trung trở lên |
| type = PT (Part-time) | -0.24 | ≈ 80.1% | Giảm nhẹ so với full-time |
| remote = 50% | -1.05 | 86.4% → ~65.6% | Giảm mạnh ~21%, do thiếu sự hiện diện trực tiếp |
| remote = 100% | -0.07 | 86.4% → ~83.6% | Ảnh hưởng rất nhẹ, remote hoàn toàn ổn định hơn hybrid |
| el = EN (Entry-level) | -1.25 | 86.4% → ~59.1% | Mới ra trường, lương thấp hơn rõ |
| el = MI (Mid-level) | -0.77 | 86.4% → ~71.7% | Thấp hơn Senior nhưng vẫn ổn |
| size = S (Small company) | -0.75 | 86.4% → ~72.0% | Công ty nhỏ trả lương thấp hơn rõ rệt |
| job_group = OT (Other) | -0.45 | 86.4% → ~76.1% | Ngành không chuyên sâu về dữ liệu thì lương thấp hơn |
Các yếu tố làm tăng xác suất lương cao
knitr::kable(
data.frame(
`Yếu tố thay đổi` = c(
"el = EX (Executive)",
"job_group = DS (Data Science)",
"job_group = ML (Machine Learning/AI)"
),
`Tác động đến z-score` = c("+0.49", "+0.22", "+0.63"),
`Ảnh hưởng đến xác suất (ước lượng)` = c(
"86.4% → ~94.2%", "86.4% → ~89.3%", "86.4% → ~92.9%"
),
`Giải thích dễ hiểu` = c(
"Lãnh đạo thường lương cao hơn",
"Chuyên ngành có giá trị thị trường tốt",
"Ngành hot nhất, trả lương cao nhất"
)
),
caption = "Các yếu tố làm tăng xác suất đạt mức lương cao"
)
| Yếu.tố.thay.đổi | Tác.động.đến.z.score | Ảnh.hưởng.đến.xác.suất..ước.lượng. | Giải.thích.dễ.hiểu |
|---|---|---|---|
| el = EX (Executive) | +0.49 | 86.4% → ~94.2% | Lãnh đạo thường lương cao hơn |
| job_group = DS (Data Science) | +0.22 | 86.4% → ~89.3% | Chuyên ngành có giá trị thị trường tốt |
| job_group = ML (Machine Learning/AI) | +0.63 | 86.4% → ~92.9% | Ngành hot nhất, trả lương cao nhất |
Muốn đạt lương cao, kỹ sư nên:
Làm full-time
Làm hoàn toàn từ xa hoặc tại công ty (tránh hybrid)
Trở thành Senior hoặc Executive
Làm trong ngành Machine Learning hoặc Data Science
Làm việc ở công ty vừa hoặc lớn
Rủi ro nhận lương thấp hơn nếu bạn là:
Freelancer, part-time, hoặc contract
Làm hybrid (50% remote)
Mới ra trường
Làm ở công ty nhỏ
Làm công việc không chuyên sâu về dữ liệu
Tiếp theo, chúng ta dự báo xác suất đạt mức lương trung, cao của kỹ sư từ bảng số liệu được tạo trước đó. Kết quả xác suất được trình bày bên dưới:
predictions <- predict(probit4, new_data, type="response")
predictions
## 1 2 3 4
## 0.8426359 0.6353944 0.8060473 0.5809197
Đầu tiên, những kỹ sư làm việc toàn thời gian, làm trực tiếp tại công ty, đạt chuyên môn cao cấp và công tác tại công ty có quy mô lớn thì có xác suất mức lương từ trung trở lên khoảng 84.26%.
Thứ hai, đối với những kỹ sư làm việc toàn thời gian, làm trực tiếp tại công ty, trình độ chuyên môn cao cấp và công tác tại công ty có quy mô nhỏ thì xác suất này giảm xuống còn khoảng 63.54%.
Thứ ba, những kỹ sư làm việc toàn thời gian, làm trực tiếp tại công ty, đạt chuyên môn trung cấp và công tác tại công ty có quy mô lớn có xác suất mức lương từ trung trở lên khoảng 80.60%.
Cuối cùng, với những kỹ sư làm việc toàn thời gian, làm trực tiếp tại công ty, trình độ trung cấp và công tác tại công ty có quy mô nhỏ thì xác suất này chỉ còn khoảng 58.09%.
Chỉ số AIC
AIC(glm4)
## [1] 15048.42
AIC(log4)
## [1] 15048.42
Chỉ số AIC mô hình probit
AIC(probit4)
## [1] 15044.95
Kết quả trên cho chúng ta thấy rằng mô hình probit sẽ hiệu quả hơn các mô hình Logistic và xác suất tuyến tính. Bởi vì chỉ số AIC của probit đạt thấp nhất là 15044.95, trong khi đó của mô hình Logistic và xác suất tuyến tính đều là 15048.42
Hệ số Brier Score
BrierScore(log4)
## [1] 0.146026
BrierScore(probit4)
## [1] 0.1459798
Dựa vào kết quả của hệ số Brier Score em cho rằng mô hình probit vẫn sẽ là lựa chọn tốt nhất vì giá trị BrierScore thấp hơn đạt 0.1459798, kết quả này tương đồng với phương pháp AIC. Trong khi đó hệ số Brier Score của mô hình Logistic là 0.146026. Brier Score của hai mô hình khá thấp, cho thấy rằng cả hai mô hình đều dự đoán khá tốt và không sai lệch nhiều so với giá trị thực tế.
Ma trận nhầm lẫn
Conf(table(predict(glm4, type="response") >=0.5,glm4$data$sal2 == '1'))
##
## Confusion Matrix and Statistics
##
## Reference
## Prediction FALSE TRUE
## FALSE 1604 934
## TRUE 2489 11507
##
## Total n : 16'534
## Accuracy : 0.7930
## 95% CI : (0.7867, 0.7991)
## No Information Rate : 0.7524
## P-Value [Acc > NIR] : < 2.2e-16
##
## Kappa : 0.3631
## Mcnemar's Test P-Value : < 2.2e-16
##
## Sensitivity : 0.3919
## Specificity : 0.9249
## Pos Pred Value : 0.6320
## Neg Pred Value : 0.8222
## Prevalence : 0.2476
## Detection Rate : 0.1535
## Detection Prevalence : 0.0970
## Balanced Accuracy : 0.6584
## F-val Accuracy : 0.4838
## Matthews Cor.-Coef : 0.3793
##
## 'Positive' Class : FALSE
Kết quả ma trận nhầm lẫn cho thấy mô hình dự đoán đúng 1604 trường hợp mức lương thấp và 11,507 trường hợp mức lương trung trở lên. Độ chính xác tổng thể đạt 79.30%, với Specificity cao (92.49%) nhưng Sensitivity thấp (39.19%), cho thấy mô hình dự đoán tốt hơn ở nhóm lương trung trở lên so với nhóm lương thấp.
Conf(table(predict(log4, type="response") >=0.5,log4$data$sal2 == '1'))
##
## Confusion Matrix and Statistics
##
## Reference
## Prediction FALSE TRUE
## FALSE 1604 934
## TRUE 2489 11507
##
## Total n : 16'534
## Accuracy : 0.7930
## 95% CI : (0.7867, 0.7991)
## No Information Rate : 0.7524
## P-Value [Acc > NIR] : < 2.2e-16
##
## Kappa : 0.3631
## Mcnemar's Test P-Value : < 2.2e-16
##
## Sensitivity : 0.3919
## Specificity : 0.9249
## Pos Pred Value : 0.6320
## Neg Pred Value : 0.8222
## Prevalence : 0.2476
## Detection Rate : 0.1535
## Detection Prevalence : 0.0970
## Balanced Accuracy : 0.6584
## F-val Accuracy : 0.4838
## Matthews Cor.-Coef : 0.3793
##
## 'Positive' Class : FALSE
Conf(table(predict(probit4, type="response") >=0.5,probit4$data$sal2 == '1'))
##
## Confusion Matrix and Statistics
##
## Reference
## Prediction FALSE TRUE
## FALSE 1618 949
## TRUE 2475 11492
##
## Total n : 16'534
## Accuracy : 0.7929
## 95% CI : (0.7867, 0.7990)
## No Information Rate : 0.7524
## P-Value [Acc > NIR] : < 2.2e-16
##
## Kappa : 0.3646
## Mcnemar's Test P-Value : < 2.2e-16
##
## Sensitivity : 0.3953
## Specificity : 0.9237
## Pos Pred Value : 0.6303
## Neg Pred Value : 0.8228
## Prevalence : 0.2476
## Detection Rate : 0.1553
## Detection Prevalence : 0.0979
## Balanced Accuracy : 0.6595
## F-val Accuracy : 0.4859
## Matthews Cor.-Coef : 0.3802
##
## 'Positive' Class : FALSE
Từ kết quả ma trận nhầm lẫn của hai mô hình logit và probit có thể thấy chúng cho kết quả gần tương đồng nhau. Cả hai mô hình đều dự đoán đúng khoảng 1.6 nghìn trường hợp mức lương thấp và trên 11.4 nghìn trường hợp mức lương trung trở lên. Tỷ lệ chính xác tổng thể của cả hai mô hình đều đạt khoảng 79.3%, cho thấy khả năng dự đoán ở nhóm lương trung trở lên tốt hơn so với nhóm lương thấp.
Sau khi tiến hành phân tích các yếu tố ảnh hưởng đến mức lương của kỹ sư với bộ dữ liệu định tính đã thu thập, em rút ra được một số kết quả quan trọng. Trước hết, xét về loại hình công việc, kỹ sư làm việc toàn thời gian (FT) chiếm tỷ lệ áp đảo trong bộ dữ liệu và được chọn làm nhóm tham chiếu. Dựa trên thống kê tần số và phân tích Relative Risk, kỹ sư làm việc theo hợp đồng (CT), bán thời gian (PT) và tự do (FL) có xác suất nhận lương thấp cao hơn đáng kể so với kỹ sư toàn thời gian, lần lượt khoảng 2.33 lần, 2.79 lần và 4.08 lần. Đặc biệt, nhóm kỹ sư tự do hoàn toàn không có khả năng đạt mức lương cao. Với mức lương trung và cao, các tỷ lệ đều nhỏ hơn 1, cho thấy kỹ sư toàn thời gian vẫn có lợi thế lớn hơn các nhóm còn lại. Kết quả hồi quy với mô hình Logit và Probit cũng củng cố kết luận này khi hai mô hình đều cho chỉ số AIC thấp và Brier Score gần 0, chứng tỏ khả năng dự đoán tốt hơn so với mô hình xác suất tuyến tính.
Về hình thức làm việc, nhóm kỹ sư làm trực tiếp tại công ty được chọn làm tham chiếu. Thống kê mô tả và phân tích Relative Risk cho thấy kỹ sư làm việc kết hợp (50% từ xa) có khả năng nhận mức lương thấp cao gấp hơn 3 lần so với nhóm làm trực tiếp, trong khi nhóm làm việc hoàn toàn từ xa cũng có nguy cơ nhận lương thấp cao hơn 1.07 lần. Đối với mức lương cao, cả hai nhóm làm từ xa và kết hợp đều có khả năng đạt được mức này thấp hơn nhóm làm trực tiếp. Phân tích hồi quy xác nhận kết quả này, đồng thời cho thấy mô hình Logit và Probit đều có khả năng dự đoán tốt hơn với AIC thấp và Brier Score nhỏ.
Xét đến trình độ chuyên môn, kỹ sư cao cấp (SE) được chọn làm nhóm tham chiếu và chiếm tỷ lệ lớn nhất trong mẫu dữ liệu. Kết quả chỉ ra rằng kỹ sư mới vào nghề (EN) có nguy cơ nhận lương thấp gấp hơn 4.7 lần so với SE, trong khi kỹ sư trung cấp (MI) có nguy cơ nhận lương thấp gấp khoảng 2.8 lần. Ngược lại, kỹ sư điều hành (EX) là nhóm duy nhất có khả năng nhận lương cao vượt trội so với SE, với odds ratio lớn hơn 1.8. Kết quả này khẳng định trình độ chuyên môn là yếu tố có tác động mạnh mẽ nhất đến thu nhập, và được phản ánh rõ trong các mô hình hồi quy, đặc biệt là Logit và Probit với độ chính xác cao.
Yếu tố quy mô công ty cũng cho thấy ảnh hưởng đáng kể đến mức lương. Nhóm kỹ sư làm việc tại công ty có quy mô trung bình được chọn làm tham chiếu. Phân tích chỉ ra rằng kỹ sư làm việc tại công ty nhỏ có nguy cơ nhận lương thấp cao gấp khoảng 2.8 lần so với nhóm trung bình, trong khi khả năng nhận lương cao chỉ bằng khoảng 0.17 lần. Ở chiều ngược lại, kỹ sư làm việc tại công ty lớn có xác suất đạt mức lương cao nhỉnh hơn nhóm trung bình, với odds ratio khoảng 1.3.
Ngoài ra, khi thêm biến nhóm công việc vào mô hình, kết quả hồi quy Logit và Probit cho thấy các hệ số của nhóm ngành có ý nghĩa thống kê và làm tăng khả năng giải thích mô hình. Kỹ sư thuộc nhóm Machine Learning (ML) có xác suất đạt lương trung và cao cao hơn đáng kể so với nhóm tham chiếu, trong khi kỹ sư thuộc nhóm công việc khác (OT) lại có xu hướng nhận lương thấp hơn. Điều này chứng tỏ lĩnh vực chuyên môn là một trong những yếu tố quyết định mạnh mẽ mức lương của kỹ sư.
Tổng hợp lại, cả bốn yếu tố: loại hình công việc, hình thức làm việc, trình độ chuyên môn và quy mô công ty, cùng với nhóm công việc đều tác động rõ rệt đến mức lương. Trong đó, trình độ chuyên môn và nhóm ngành nghề nổi bật là các yếu tố có ảnh hưởng mạnh nhất. Các mô hình Logit và Probit tiếp tục thể hiện sự phù hợp vượt trội so với mô hình xác suất tuyến tính với độ chính xác dự đoán cao, chỉ số AIC thấp và Brier Score nhỏ, khẳng định khả năng giải thích tốt mối quan hệ giữa các yếu tố này và thu nhập của kỹ sư.
Dữ liệu sử dụng chủ yếu là định tính và được mã hóa, nên chưa phản ánh chi tiết các yếu tố ảnh hưởng đến lương kỹ sư. Một số biến quan trọng như số năm kinh nghiệm hay ngành nghề cụ thể chưa được đưa vào, làm giảm độ chính xác của kết quả.
Mô hình chỉ xem xét trong một thời điểm, chưa tính đến sự thay đổi của thị trường lao động theo thời gian. Điều này khiến kết quả khó phản ánh được xu hướng dài hạn.
Dữ liệu tập trung vào một nhóm kỹ sư nhất định, nên kết quả khó khái quát cho các ngành khác hoặc các khu vực lao động khác.
Các mô hình sử dụng như Logit, Probit và mô hình xác suất tuyến tính vẫn là mô hình tuyến tính cơ bản, chưa bắt được các mối quan hệ phức tạp giữa các yếu tố ảnh hưởng đến mức lương.
Hướng nghiên cứu tiếp theo có thể mở rộng theo một số khía cạnh. Trước hết, bổ sung thêm các biến định lượng như số năm kinh nghiệm, số giờ làm việc hoặc hiệu suất để đánh giá tác động cụ thể hơn đến mức lương. Bên cạnh đó, có thể thử nghiệm các mô hình nâng cao như Multinomial Logit/Probit để phân tích ba mức lương riêng biệt, hoặc áp dụng các thuật toán Machine Learning nhằm tăng khả năng dự báo. Ngoài ra, việc kết hợp dữ liệu theo thời gian (panel data) và đưa thêm yếu tố ngành nghề, khu vực hay điều kiện thị trường vào phân tích sẽ giúp mô hình phản ánh thực tế hơn. Cuối cùng, mở rộng so sánh giữa các ngành hoặc các quốc gia cũng là một hướng quan trọng để làm rõ các yếu tố ảnh hưởng đến mức lương kỹ sư.