Ordinal Logistic Regression

R Econometrics Series

Data Cleaning

#==================================
#  R Code xử lí làm sạch dữ liệu 
#  cho Ordinal Logistic Regression
#==================================


# Clear môi trường R: 
rm(list = ls())

# Load một số packages cho xử lí - làm sạch dữ liệu: 

library(haven)
library(dplyr)
library(ggplot2)
library(stringi)
library(stringr)
library(writexl)
library(scales)
library(tidyr)
library(writexl)

# Chuẩn bị font chữ và color: 

colorFTU <- "#cc0000"

library(showtext) # -> Package for using extra fonts. 

font_subtitle <- "Roboto Condensed"  # -> Set Outfit font for our plot. 

font_add_google(name = font_subtitle, family = font_subtitle) # -> Load font for using. 

font_text <- "Oswald"

font_add_google(name = font_text, family = font_text)

font_main <- "Ubuntu"  

font_add_google(name = font_main, family = font_main) 

# Automatically render text: 
showtext_auto()


# Survey data về CMCN: 

readxl::read_xlsx("E:\\Project_LabourStructureThesis\\DATA-_Dieu tra Cong nghe 2019.xlsx") -> dataCMCN 

# Hiệu chỉnh mã số thuế: 

dataCMCN %>% mutate(Masothue2 = case_when(is.na(Masothue2) ~ "", TRUE ~ Masothue2)) -> dataCMCN

# Xử lí mã thuế: 

dataCMCN %>% 
  mutate(taxCode = str_c(Masothue, Masothue2)) -> dataCMCN


# Load phiếu 1A: 

read_dta("E:\\Project_LabourStructureThesis\\Data\\Data GSO_DN1620\\2020Dung/TQ_ASR_DN.dta") -> phieu1A

nrow(phieu1A) # Bộ dữ liệu này rất lớn. 

# Hàm chuẩn hóa mã số thuế: 

fullCodeTax <- function(dfSelected) {
  
  dfSelected$Masothue -> a
  
  dfSelected$Masothue2 -> b
  
  str_c(a, b) -> fullCode
  
  return(dfSelected %>% mutate(taxCode = fullCode))
  
}

# Hiệu chỉnh mã số thuế của phiếu 1A: 

phieu1A %>% fullCodeTax() -> phieu1A 

# Bổ sung kí hiệu "Phieu1A": 

str_c(names(phieu1A), "Phieu1A", sep = ".") -> newNames

names(phieu1A) <- newNames

phieu1A %>% rename(taxCode = taxCode.Phieu1A) -> phieu1A

# Combine dữ liệu: 

dataCMCN %>% left_join(phieu1A, by = "taxCode") -> dataCMCN


#---------------------------------------------------
#  Block 2: Recode và xử lí các biến số quan trọng
#---------------------------------------------------

#----  1. Loại hình DN theo tỉ lệ sở hữu --------


# Xem mô tả trong file Phieu1A, mục Questionaire: 

dataCMCN %>% 
  group_by(LoaihinhKT.Phieu1A) %>% 
  count(sort = TRUE)


dataCMCN %>% 
  mutate(LoaiDN = case_when(LoaihinhKT.Phieu1A %in% 1:3 ~ "NNTren50Von", 
                            LoaihinhKT.Phieu1A %in% c(4, 8, 10) ~ "NNDuoi50Von", 
                            LoaihinhKT.Phieu1A == 5 ~ "HopTacXa",
                            LoaihinhKT.Phieu1A %in% c(6, 7, 9) ~ "TuNhan", 
                            # is.na(LoaihinhKT.Phieu1A) ~ "Unknown", 
                            TRUE ~ "CoNN")) -> dataCMCN

# Tần suất của loại hình DN: 

dataCMCN %>% 
  group_by(LoaiDN) %>% 
  count(sort = TRUE)

#----- 2. Về triển  khai và thực thi chiến lược CMCN  --------

# Đọc kĩ Phieu2CMCN, Questionaire: 

dataCMCN %>% 
  group_by(ChienLuocCMCN) %>% 
  count()

# Recode lại ChienLuocCMCN: 

dataCMCN %>% 
  mutate(ChienLuocCMCN = case_when(ChienLuocCMCN == 1 ~ "ChuaHieu1", 
                                   ChienLuocCMCN == 2 ~ "DaHieu2", 
                                   ChienLuocCMCN == 3 ~ "YTuong3", 
                                   ChienLuocCMCN == 4 ~ "XayDung4", 
                                   TRUE ~ "DangTH5")) -> dataCMCN 

# Convert về factor: 


dataCMCN %>% 
  mutate(ChienLuocCMCN = factor(ChienLuocCMCN, levels = c("ChuaHieu1", 
                                                          "DaHieu2", 
                                                          "YTuong3", 
                                                          "XayDung4", 
                                                          "DangTH5"))) -> dataCMCN
# Kiểm tra lại: 


dataCMCN %>% 
  group_by(ChienLuocCMCN) %>% 
  count()


# Tạo biến mới ChienLuocAdj: 


dataCMCN <- dataCMCN %>%
  mutate(
    ChienLuocAdj = case_when(
      ChienLuocCMCN %in% c("ChuaHieu1", "DaHieu2", "YTuong3") ~ "ChuaTrienKhai",
      TRUE ~ "DangTrienKhai"
    )
  )

# Tính toán số lượng công nghệ: 

dataCMCN %>% 
  dplyr::select(contains("CongNghe")) %>% 
  dplyr::select(-CongNghe_KhongSuDung) %>% 
  rowSums() -> SoLuongCNSuDung


dataCMCN %>% 
  mutate(SoLuongCNSuDung = SoLuongCNSuDung) -> dataCMCN


# Số lượng công nghệ 4.0 đang sử dụng: 


dataCMCN %>% 
  dplyr::select(contains("HienTrang")) %>% 
  mutate_all(function(x) {case_when(x %in% c(1, 2, 6) ~ 0, TRUE ~ 1)}) %>% 
  rowSums() -> SoLuongCN4.0


dataCMCN %>% mutate(SoLuongCN4.0 = SoLuongCN4.0) -> dataCMCN


# Tính toán cơ cấu lao động theo giới tính: 

dataCMCN %>% 
  mutate(TiLeNu = SoLaodong_CuoiNam_Nu.Phieu1A / SoLaodong_CuoiNam.Phieu1A) %>% 
  mutate(TiLeNam = 1 - TiLeNu) -> dataCMCN


# Tỉ lệ lao động theo trình độ: 

dataCMCN %>% 
  select(DaiHoc_HienNay) %>% 
  summary() # Tỉ lệ

dataCMCN %>% 
  select(SoCap_HienNay) %>%  
  summary() # Tỉ lệ


dataCMCN %>% 
  select(ChuaDaoTao_HienNay) %>% 
  summary()


dataCMCN %>% select(contains("ienNay"))

dataCMCN %>% 
  mutate(DaiHocTiLe = DaiHoc_HienNay, 
         SoCapTiLe = SoCap_HienNay, 
         ChuaDaoTaoTiLe = ChuaDaoTao_HienNay) -> dataCMCN


# Kì vọng của doanh nghiệp về thay đổi cơ cấu lao động theo trình độ: 

dataCMCN %>% 
  select(DaiHoc_HienNay, DaiHoc_12Nam, SoCap_12Nam, ChuaDaoTao_12Nam)



dataCMCN %>% 
  mutate(DaiHocNamToi = case_when(DaiHoc_12Nam == 1 ~ "Tang", 
                                  DaiHoc_12Nam == 2 ~ "KhongDoi", 
                                  TRUE ~ "Giam")) %>% 
  mutate(SoCapNamToi = case_when(SoCap_12Nam == 1 ~ "Tang", 
                                 SoCap_12Nam == 2 ~ "KhongDoi", 
                                 TRUE ~ "Giam")) %>% 
  mutate(ChuaDTNamToi = case_when(ChuaDaoTao_12Nam == 1 ~ "Tang", 
                                  ChuaDaoTao_12Nam == 2 ~ "KhongDoi", 
                                  TRUE ~ "Giam")) -> dataCMCN

dataCMCN %>% 
  mutate(DaiHocNamToi = factor(DaiHocNamToi, levels = c("Tang", "KhongDoi", "Giam"))) %>% 
  mutate(SoCapNamToi = factor(SoCapNamToi, levels = c("Tang", "KhongDoi", "Giam"))) %>% 
  mutate(ChuaDTNamToi = factor(ChuaDTNamToi, levels = c("Tang", "KhongDoi", "Giam"))) -> dataCMCN



dataCMCN %>% 
  mutate(DaoTaoNV = case_when(CoDaoTao == 1 ~ "Co", 
                              CoDaoTao == 2 ~ "KeHoach", 
                              TRUE ~ "Khong")) -> dataCMCN


dataCMCN %>% 
  mutate(DaoTaoNV = factor(DaoTaoNV, levels = c("Co", "KeHoach", "Khong"))) -> dataCMCN

convertKyNang <- function(x) {
  
  case_when(x == 1 ~ "Tot", x == 2 ~ "DapUng", TRUE ~ "Yeu") -> x_converted
  
  factor(x_converted, levels = c("Tot", "DapUng", "Yeu")) %>% return()
  
}


dataCMCN %>% 
  mutate(KyNangCNTT = convertKyNang(KyNang_CNTT), 
         KyNangTuDongHoa = convertKyNang(KyNang_TuDongHoa), 
         KyNangPhanTichDuLieu = convertKyNang(KyNang_PhanTichDuLieu), 
         KyNangAnNinhDL = convertKyNang(KyNang_AnNinhDuLieu), 
         KyNangHTHoTro = convertKyNang(KyNang_HeThongHoTro), 
         KyNangPhanMemCT = convertKyNang(KyNang_PhanMemCongTac), 
         KyNangPhiKyThuat = convertKyNang(KyNang_PhiKyThuat)) -> dataCMCN



dataCMCN %>% 
  mutate(TrienKhai4.0 = ChienLuocAdj, 
         DoanhThuThuan = Doanhthu_Thuan_BH_CCDV.Phieu1A, 
         LoiNhuanSauThue = SR_Loinhuan_SauThue.Phieu1A, 
         ChiPhiLaoDong = Tongtien_Tra_NguoiLD.Phieu1A) -> dfOLR


dfOLR %>% 
  names() %>% 
  str_replace_all("\\.", "") %>% 
  str_sub(1, 32) -> newCols

names(dfOLR) <- newCols

dfOLR %>% 
  mutate(DoanhThuThuan = DoanhThuThuan / 1000, 
         LoiNhuanSauThue = LoiNhuanSauThue / 1000, 
         ChiPhiLaoDong = ChiPhiLaoDong / 1000) -> dfOLR


# Save Data: 

write_dta(dfOLR, "dfOLR.dta") 

saveRDS(dfOLR, "dfOLR.rds")

Ordinal Logistic Regression

#====================================
#     Ordinal Logistic Regression
#====================================

# Load data: 
rm(list = ls())
readRDS("dfOLR.rds") -> dfOLR

# Load một số package cần thiết: 

library(ordinal)
library(dplyr)
library(kableExtra)
library(broom)

# Convert về đơn vị Tỉ VND: 

dfOLR %>% 
  mutate(ChiPhiLaoDong = ChiPhiLaoDong / 1000, 
         DoanhThuThuan = DoanhThuThuan / 1000, 
         LoiNhuanSauThue = LoiNhuanSauThue / 1000) -> dfOLR

# Hàm tạo mô hình: 

model_OLR <- function(bienY) {
  
  biendoclap <- "~ LoaiDN + SoLuongCNSuDung + SoLuongCN40 + ChiPhiLaoDong + DoanhThuThuan + LoiNhuanSauThue + 
  DaoTaoNV + KyNangCNTT + KyNangTuDongHoa + KyNangPhanTichDuLieu + KyNangAnNinhDL + 
  KyNangHTHoTro + KyNangPhanMemCT + KyNangPhiKyThuat"
  
  formula(paste(bienY, biendoclap, sep = " ")) -> myModel
  
  return(myModel)
  
}

# Mô hình OLR: 

model_OLR("DaiHocNamToi") -> regDaiHoc

m_clm <- clm(regDaiHoc, data = dfOLR, link = "logit", Hess = TRUE)


# Hàm trình bày kết quả theo mẫu: 
# 1. https://www.ub.edu/rfa/R/RCDV_results2.pdf
# 2. https://pmc.ncbi.nlm.nih.gov/articles/PMC8831068
# 3. https://web.pdx.edu/~newsomj/cdaclass/ho_ordinal%20regression.pdf


trinhbayketqua <- function(objectOLR) {
  
  tidy(objectOLR) %>% 
  filter(coef.type == "location") %>% 
  mutate(odds_ratio = exp(estimate),
         conf.low   = exp(estimate - 1.96 * std.error),
         conf.high  = exp(estimate + 1.96 * std.error)) %>% 
  select(term, estimate, std.error, statistic, p.value, odds_ratio, conf.low, conf.high) %>% 
  mutate_if(is.numeric, function(x) {round(x, 3)}) %>% 
    return()
  
}

# Sử dụng hàm: 

m_clm %>% 
  trinhbayketqua() %>% 
  kbl(caption = "Table 1: OLR Results, DaiHocNamToi là biến phụ thuộc") %>%
  kable_classic(full_width = FALSE, html_font = "Cambria") 

Table 1: OLR Results, DaiHocNamToi là biến phụ thuộc

term	estimate	std.error	statistic	p.value	odds_ratio	conf.low	conf.high
LoaiDNHopTacXa	0.134	0.220	0.607	0.544	1.143	0.742	1.761
LoaiDNNNDuoi50Von	-0.457	0.179	-2.554	0.011	0.633	0.446	0.899
LoaiDNNNTren50Von	-0.542	0.217	-2.497	0.013	0.581	0.380	0.890
LoaiDNTuNhan	-0.402	0.179	-2.244	0.025	0.669	0.471	0.950
SoLuongCNSuDung	-0.081	0.030	-2.703	0.007	0.922	0.869	0.978
SoLuongCN40	-0.034	0.019	-1.809	0.071	0.967	0.932	1.003
ChiPhiLaoDong	-0.901	0.395	-2.279	0.023	0.406	0.187	0.881
DoanhThuThuan	0.035	0.016	2.217	0.027	1.035	1.004	1.068
LoiNhuanSauThue	0.015	0.083	0.182	0.856	1.015	0.863	1.195
DaoTaoNVKeHoach	0.098	0.072	1.362	0.173	1.103	0.958	1.270
DaoTaoNVKhong	0.973	0.074	13.167	0.000	2.645	2.288	3.057
KyNangCNTTDapUng	0.103	0.100	1.030	0.303	1.109	0.911	1.349
KyNangCNTTYeu	0.502	0.114	4.389	0.000	1.653	1.321	2.068
KyNangTuDongHoaDapUng	-0.364	0.163	-2.233	0.026	0.695	0.505	0.957
KyNangTuDongHoaYeu	-0.418	0.166	-2.515	0.012	0.659	0.476	0.912
KyNangPhanTichDuLieuDapUng	0.039	0.167	0.232	0.817	1.040	0.749	1.443
KyNangPhanTichDuLieuYeu	0.010	0.175	0.056	0.955	1.010	0.716	1.424
KyNangAnNinhDLDapUng	0.073	0.139	0.530	0.596	1.076	0.820	1.412
KyNangAnNinhDLYeu	-0.013	0.150	-0.084	0.933	0.988	0.736	1.326
KyNangHTHoTroDapUng	0.049	0.178	0.277	0.782	1.050	0.742	1.488
KyNangHTHoTroYeu	0.219	0.187	1.168	0.243	1.245	0.862	1.797
KyNangPhanMemCTDapUng	-0.096	0.157	-0.613	0.540	0.908	0.668	1.235
KyNangPhanMemCTYeu	-0.215	0.166	-1.294	0.196	0.807	0.583	1.117
KyNangPhiKyThuatDapUng	0.243	0.174	1.393	0.164	1.275	0.906	1.795
KyNangPhiKyThuatYeu	0.243	0.182	1.338	0.181	1.275	0.893	1.822

Diễn giải kết quả

Kết quả ước lượng từ mô hình Ordinal Logistic Regression cho thấy kỳ vọng của doanh nghiệp đối với lao động có trình độ đại học trong 1–2 năm tới chịu ảnh hưởng bởi nhiều yếu tố liên quan đến loại hình doanh nghiệp, mức độ sử dụng công nghệ, điều kiện tài chính, chiến lược đào tạo nhân lực và tình trạng kỹ năng hiện tại. Do biến phụ thuộc DaiHocNamToi được đo lường theo thứ bậc, các hệ số được diễn giải thông qua odds ratio, phản ánh khả năng doanh nghiệp kỳ vọng mức cao hơn (tăng so với không đổi/giảm), trong điều kiện các yếu tố khác không đổi.

Trước hết, xét theo loại hình doanh nghiệp, kết quả cho thấy sự khác biệt rõ rệt giữa doanh nghiệp có vốn nước ngoài và phần lớn các loại hình doanh nghiệp trong nước. Cụ thể, doanh nghiệp nhà nước có tỷ lệ vốn dưới 50% có odds kỳ vọng tăng lao động trình độ đại học chỉ bằng 0.63 lần so với doanh nghiệp có vốn nước ngoài, và sự khác biệt này có ý nghĩa thống kê ở mức 5%. Tương tự, doanh nghiệp nhà nước có tỷ lệ vốn trên 50% và doanh nghiệp tư nhân cũng có odds kỳ vọng tăng lao động đại học thấp hơn đáng kể, lần lượt bằng 0.58 lần và 0.67 lần so với doanh nghiệp có vốn nước ngoài. Trong khi đó, hợp tác xã không cho thấy sự khác biệt có ý nghĩa thống kê. Những kết quả này cho thấy doanh nghiệp có vốn nước ngoài có xu hướng lạc quan hơn trong việc mở rộng nhu cầu lao động trình độ đại học, phản ánh sự khác biệt về mô hình kinh doanh, cường độ công nghệ và tiêu chuẩn nhân lực.

Liên quan đến mức độ sử dụng công nghệ, mô hình chỉ ra sự phân hóa rõ ràng giữa công nghệ truyền thống và công nghệ 4.0. Biến SoLuongCNSuDung, đại diện cho số lượng công nghệ kiểu cũ mà doanh nghiệp đang sử dụng, có tác động âm và có ý nghĩa thống kê cao. Cụ thể, khi doanh nghiệp sử dụng thêm một công nghệ truyền thống, odds kỳ vọng tăng lao động trình độ đại học giảm khoảng 7.8%. Điều này cho thấy các doanh nghiệp phụ thuộc nhiều vào công nghệ cũ có xu hướng duy trì cơ cấu lao động hiện tại và ít kỳ vọng mở rộng lao động trình độ cao trong tương lai gần. Ngược lại, biến SoLuongCN40 – số lượng công nghệ thuộc nhóm 4.0 – không có ý nghĩa thống kê, cho thấy việc áp dụng công nghệ 4.0 về mặt số lượng chưa đủ để làm thay đổi rõ rệt kỳ vọng ngắn hạn về lao động đại học.

Xét đến các yếu tố tài chính, kết quả cho thấy những tác động khá đáng chú ý. ChiPhiLaoDong có hệ số âm và có ý nghĩa thống kê, với odds ratio bằng 0.41, hàm ý rằng doanh nghiệp có chi phí lao động cao hơn có xu hướng giảm kỳ vọng tăng lao động trình độ đại học. Ngược lại, DoanhThuThuan có tác động dương và có ý nghĩa thống kê, cho thấy khi doanh thu tăng, odds doanh nghiệp kỳ vọng tăng lao động đại học cũng tăng lên. Trong khi đó, LoiNhuanSauThue không có ý nghĩa thống kê, gợi ý rằng kỳ vọng tuyển dụng lao động trình độ cao trong ngắn hạn phụ thuộc nhiều hơn vào dòng doanh thu và quy mô hoạt động hơn là mức lợi nhuận ròng.

Về đào tạo nhân viên, kết quả cho thấy một phát hiện rất nổi bật. Doanh nghiệp không đào tạo và cũng không có kế hoạch đào tạo có odds kỳ vọng tăng lao động trình độ đại học cao gấp 2.65 lần so với doanh nghiệp đang đào tạo, và kết quả này có ý nghĩa thống kê rất cao. Điều này cho thấy các doanh nghiệp không đầu tư đào tạo nội bộ có xu hướng kỳ vọng tuyển mới lao động có trình độ đại học từ thị trường lao động bên ngoài, thay vì nâng cao kỹ năng cho lực lượng lao động hiện hữu. Ngược lại, doanh nghiệp có kế hoạch đào tạo nhưng chưa triển khai chưa cho thấy sự khác biệt có ý nghĩa so với doanh nghiệp đang đào tạo.

Cuối cùng, xét theo các nhóm kỹ năng, kết quả cho thấy tác động không đồng đều. Đối với kỹ năng CNTT, doanh nghiệp đánh giá kỹ năng ở mức yếu có odds kỳ vọng tăng lao động đại học cao hơn 1.65 lần so với doanh nghiệp đánh giá kỹ năng là tốt. Điều này phản ánh vai trò của lao động trình độ đại học trong việc bù đắp sự thiếu hụt kỹ năng CNTT hiện tại. Ngược lại, đối với kỹ năng tự động hóa, các doanh nghiệp có mức kỹ năng đáp ứng hoặc yếu lại có odds kỳ vọng tăng lao động đại học thấp hơn so với doanh nghiệp có kỹ năng tốt, cho thấy năng lực tự động hóa cao đi kèm với nhu cầu mở rộng lao động trình độ cao. Các nhóm kỹ năng khác như phân tích dữ liệu, an ninh dữ liệu, hệ thống hỗ trợ, phần mềm chuyên dụng và kỹ năng phi kỹ thuật không cho thấy tác động có ý nghĩa thống kê.

Tổng thể, kết quả mô hình cho thấy kỳ vọng của doanh nghiệp về lao động trình độ đại học trong tương lai gần chịu ảnh hưởng mạnh bởi loại hình sở hữu, mức độ phụ thuộc vào công nghệ truyền thống, điều kiện tài chính và chiến lược nhân sự, trong khi tác động của nhiều nhóm kỹ năng công nghệ khác vẫn còn hạn chế trong ngắn hạn.