Data Pre-processing: Case of VES (Vietnam Enterprise Survey)

Data Pre-processing Series

Our Problem

Tiền xử lí số liệu (Data Pre-processing) là một khâu mất nhiều thời gian. Post này sử dụng bộ dữ liệu VES - Vietnam Enterprise Survey (download tại đây). Giải nén và đọc bộ dữ liệu này:

#===================================================================
#  Data Processing Project with real-world data: VES 2015 data set
#===================================================================

# Clear work space: 
rm(list = ls())

# Import data: 

haven::read_dta("E:/dn2015.dta") -> dn2015

# Load dplyr gackage: 
library(dplyr)

# Number of columns/rows: 
dim(dn2015)

## [1] 455300    222

Bộ dữ liệu này là lớn với 455300 quan sát và 222 cột biến. Theo quy định của các cơ quan quản lí tại Việt Nam thì:

1. Mỗi một DN chỉ có một mã số thuế duy nhất và mã số thuế này là độc nhất.
2. Một mã số thuế hợp lệ dùng 10 chữ số.

Tuy nhiên một số MST (mã số thuế) trong bộ dữ liệu chỉ có 5 chữ số. Mặt khác, một số DN có MST trùng nhau. Chúng ta sẽ loại bỏ các quan sát này ra khỏi bộ dữ liệu ban đầu như sau:

#---------- Remove firms with duplicated tax code ---------------

dn2015 %>% 
  group_by(ma_thue) %>% 
  count() %>% 
  ungroup() %>% 
  filter(n > 1) %>% 
  pull(ma_thue) -> tax_codes_duplicated

dn2015 %>% filter(!ma_thue %in% tax_codes_duplicated) -> dn2015

#--------- Remove firms with tax code != 10 digits ----------------

library(stringr)

dn2015 %>% filter(str_count(ma_thue) == 10) -> dn2015

Theo Quyết định số 27/2018/QĐ-TTg thì mã ngành cấp 5 của doanh nghiệp sẽ phải có 5 chữ số. Tuy nhiên một số quan sát thì mã ngành chỉ có 4 chữ số. Ta có thể suy luận rằng nguyên nhân là do một số ngành có mã với số 0 đứng ở vị trí đầu tiên đã bị “biến mất” trong quá trình vào-lưu-gửi dữ liệu. Vì vậy chúng ta có thể khôi phục chứ số 0 này cho các quan sát với mã ngành chỉ có 4 chữ số như sau:

# Add 0 for 4-digits sector codes: 

dn2015 %>% 
  mutate(nganh_kd = as.character(nganh_kd)) %>% 
  mutate(nganh_kd  = case_when(str_count(nganh_kd) == 4 ~ str_c("0", nganh_kd), TRUE ~ nganh_kd)) -> dn2015

Chúng ta có thể xem qua dữ liệu sau khi xử lí:

library(knitr)

dn2015 %>% 
  head() %>% 
  select(1:7, contains("nganh")) %>% 
  kable(caption = "Table 1: Some Observations")

Table 1: Some Observations

tinh	capso	macs	ma_thue	huyen	xa	tthd	tennganhkd	nganh_kd
01	1	152	0100110528	001	00028	1	X©y Dùng	42200
01	1	465	0100100583	007	00283	1	DÖt Kim	14100
01	1	565	0100100752	020	00640	1	Sx Pin	27200
01	1	1825	0102071754	009	00355	1	Sx C¸c Sp Gç	16291
01	1	104752	0100886857	019	00625	1	X©y Dùng	42200
01	1	277	0100103785	006	00232	1	X©y Dùng	42200

Lấy dữ liệu về mã ngành của Quyết định số 27/2018/QĐ-TTg bằng đoạn R codes dưới đây:

# --------------------------------------------------------------------------------
# Extract sector code from https://dangkykinhdoanh.gov.vn/vn/Pages/NganhNghe.aspx
#---------------------------------------------------------------------------------

library(rvest)

link <- "https://dangkykinhdoanh.gov.vn/vn/Pages/NganhNghe.aspx"

link %>% 
  read_html() %>% 
  html_nodes(xpath = '//*[@id="ctl00_SPWebPartManager1_g_15a30b99_50d8_49e5_8f11_de265cb949e8"]/div/div[2]/table') %>% 
  html_table() %>% 
  .[[1]] -> df_sector_code

# Rename for columns: 

names(df_sector_code) <- c(str_c("code_level", 1:5), "sector_name")

# Remove/ filter data: 

df_sector_code %>% 
  slice(-1) %>% 
  filter(code_level1 != "21") %>% 
  mutate(sector_name_latin = stringi::stri_trans_general(sector_name, "Latin-ASCII")) -> df_sector_code

Có 21 ngành lớn (Cấp 1) theo quy định của chính phủ như ta thấy (chỉ hiển thị 6 ngành cấp 1):

df_sector_code %>% 
  filter(str_count(code_level1) != 0) %>% 
  select(code_level1, sector_name_latin) %>% 
  head() %>% 
  kable(caption = "Table 2: Some Level-1 Sectors")

Table 2: Some Level-1 Sectors

code_level1	sector_name_latin
A	NONG NGHIEP, LAM NGHIEP VA THUY SAN
B	KHAI KHOANG
C	CONG NGHIEP CHE BIEN, CHE TAO
D	SAN XUAT VA PHAN PHOI DIEN, KHI DOT, NUOC NONG, HOI NUOC VA DIEU HOA KHONG KHI
E	CUNG CAP NUOC; HOAT DONG QUAN LY VA XU LY RAC THAI, NUOC THAI
F	XAY DUNG

Để lấp đầy dữ liệu trống cho cột code_level1 chúng ta có thể sử dụng tidyr::fill() như sau:

df_sector_code %>% 
  mutate(code_level1 = case_when(str_count(code_level1) == 0 ~ NA_character_, TRUE ~ code_level1)) %>% 
  tidyr::fill(code_level1, .direction = "down") %>% 
  filter(str_detect(sector_name_latin, "[a-z]")) -> df_sector_code_filled

So sánh:

# Original data: 

df_sector_code %>% 
  select(1:5) %>% 
  head() %>% 
  kable(caption = "Table 3: Original Data")

Table 3: Original Data

code_level1	code_level2	code_level3	code_level4	code_level5
A
	01
		011
			0111	01110
			0112	01120
			0113	01130

# Filled data: 
df_sector_code_filled %>% 
  select(1:5) %>% 
  head() %>% 
  kable(caption = "Table 4: Filled Data")

Table 4: Filled Data

code_level1	code_level2	code_level3	code_level4	code_level5
A	01
A		011
A			0111	01110
A			0112	01120
A			0113	01130
A			0114	01140

Theo Quyết định số 27/2018/QĐ-TTg thì có tất cả 734 ngành cấp 5:

df_sector_code_filled %>% filter(str_count(code_level5) == 5) -> df_sector_code_level5

df_sector_code_level5 %>% pull(code_level5) -> code_level5

length(code_level5)

## [1] 734

Nhưng số lượng mã ngành cấp 5 của VES 2015 chỉ là 621:

dn2015 %>% filter(!is.na(nganh_kd)) -> dn2015

dn2015 %>% 
  pull(nganh_kd) %>% 
  unique() -> code_level5_ves

length(code_level5_ves)

## [1] 621

Dưới đây là 60 mã ngành ở VES 2015 không có mặt trong danh sách mã ngành của năm 2018:

condition_in2018 <- code_level5_ves %in% code_level5

code_level5_ves[!condition_in2018] -> not_in_2018

not_in_2018

##  [1] "42200" "42102" "41000" "22120" "42101" "42900" "10750" "13220" "35101"
## [10] "35102" "20110" "68100" "13290" "68200" "24200" "23100" "52219" "61200"
## [19] "61100" "79200" "10409" "58190" "85311" "85100" "58110" "85322" "72100"
## [28] "63290" "52211" "49200" "85200" "01190" "65129" "02210" "01290" "01300"
## [37] "01410" "02220" "01450" "03230" "03122" "13210" "02109" "47110" "63210"
## [46] "10204" "13240" "03210" "72200" "58130" "85312" "13230" "51100" "02300"
## [55] "85321" "35200" "01440" "01420" "58120" "03121"

Như vậy một số tên ngành của VES 2015 không có mặt trong Quyết định số 27/2018/QĐ-TTg mới nhất. Việc này có thể là do: một số mã ngành của năm 2015 sau này được phân loại chi tiết hơn thành nhiều hơn thành 2, 3 ngành nhỏ hơn hoặc một số mã ngành ở 2015 đã được recode lại theo một cách thức nào đó mà ta chưa rõ. Trước khi đi đến một giải pháp tinh tế hơn, thì tạm thời chúng ta phân loại những ngành với mã ngành không có mặt trong danh sách của Quyết định số 27/2018/QĐ-TTg thành một nhóm chung có tên Unclassified. Trước hết join các bộ dữ liệu lại:

# join the two data sets: 
dn2015 %>% full_join(
  df_sector_code_level5 %>% select(sector_name_latin, code_level5), 
  by = c("nganh_kd" = "code_level5")
  ) -> dn2015_joined

# Relabel: 

dn2015_joined %>% mutate(sector_name_latin = case_when(is.na(sector_name_latin) ~ "Unclassified", TRUE ~ sector_name_latin)) -> dn2015_joined

Hàm đánh giá tỉ lệ missing của dữ liệu:

missing_rate <- function(x) {sum(is.na(x)) / length(x)}

Hàm đánh giá tị lệ số âm:

negative_rate <- function(x) {sum(x < 0, na.rm = TRUE) / length(x)}

Tỉ lệ missing và số âm của cột biến kqkd1:

dn2015_joined$kqkd1 %>% missing_rate()

## [1] 0.03396914

dn2015_joined$kqkd1 %>% negative_rate()

## [1] 6.676325e-06

Giải pháp ngây thơ là loại bỏ các quan sát là số âm, missing ở cột biến kqkd1:

dn2015_joined %>% 
  filter(!is.na(kqkd1)) %>% 
  filter(kqkd1 >= 0) -> dn2015_joined1

Đến đây chúng ta có thể, ví dụ, tính tổng doanh thu theo ngành cấp 5:

dn2015_joined1 %>% 
  group_by(sector_name_latin) %>% 
  summarise(sales_by_sector = sum(kqkd1)) %>% 
  ungroup() %>% 
  arrange(-sales_by_sector) %>% 
  head() %>% 
  kable(caption = "Table 5: Total Revenue by Sector")

Table 5: Total Revenue by Sector

sector_name_latin	sales_by_sector
Unclassified	1476867424
San xuat thiet bi truyen thong	928331927
Ban buon xang dau va cac san pham lien quan	484522160
Ban buon sat, thep	329661497
May trang phuc (tru trang phuc tu da long thu)	255777993
San xuat giay, dep	219492935

Missing Data

Dữ liệu trống (Missing Data) là vấn đề luôn gặp ở gần như mọi dự án phân tích dữ liệu. Trước hết đánh giá tỉ lệ missing cho tất cả các cột biến của VES 2015:

sapply(dn2015_joined1, missing_rate) -> missing_columns

df_miss_rate <- tibble(var = names(dn2015_joined1), missing_rate_col = missing_columns)

df_miss_rate %>% arrange(-missing_rate_col) -> df_miss_rate

df_miss_rate %>% 
  head() %>% 
  kable(caption = "Table 5: Missing Rate by Column")

Table 5: Missing Rate by Column

var	missing_rate_col
kqkdn13	0.9873964
nkqkd13	0.9873872
nkd13	0.9853461
ldn131	0.9853115
nld13	0.9853069
ldn132	0.9853000

Từ table 5 có thể thấy tỉ lệ missing có thể là rất lớn: gần 100%. Dưới đây là danh sách (một số) cột biến có tỉ lệ missing bé hơn 50%:

df_miss_rate %>% 
  filter(missing_rate_col < 0.5) -> mising_under50

mising_under50 %>% 
  head() %>% 
  kable(caption = "Table 6: Columns with Missing Rate <= 50%")

Table 6: Columns with Missing Rate <= 50%

var	missing_rate_col
von_nn	0.4991638
khucn	0.4990601
tgnk_tt	0.4895711
kqkdn1	0.4794002
nkd1	0.4789095
ldn11	0.4787874

Chúng ta cũng có thể sử dụng gói naniar để đánh giá tỉ lệ missing của data một cách nhanh chóng:

library(naniar)

# Missing rate for 25 columns: 
gg_miss_var(dn2015_joined1 %>% select(1:25), show_pct = TRUE) 

Conclusion

Data Pre-processing/Data Wrangling là khâu mất nhiều thời gian của một dự án phân tích dữ liệu. Post này không thể cover hết mọi khía cạnh của quá trình này.