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

Data Pre-processing Series, FTU Short Course

Problem 1: Tax Code

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 được cung cấp bởi Dr. Mai Vu. Đọc và xem qua một số thông tin về bộ dữ liệu này:

#===================================================================
#  Data Processing Project with real-world data: VES 2015 data set
#  Data Source: maivp@ftu.edu.vn
#===================================================================

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

# Import data: 

haven::read_dta("F:\\VES_from_MaiVu_FTU\\Stata_2015\\dn2015.dta") -> dn2015

# Load dplyr gackage: 
library(dplyr)

# Number of columns/rows: 
dim(dn2015)

## [1] 455238    150

Bộ dữ liệu này là lớn với 455238 quan sát và 150 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) %>% 
  summarise(tansuat = n()) -> df_tansuat_mathue

df_tansuat_mathue %>% 
  filter(tansuat > 1) %>% 
  pull(ma_thue) -> tax_codes_duplicated

# Remove duplications: 

dn2015 %>% filter(!duplicated(ma_thue)) -> dn2015

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

library(stringr)

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

Problem 2: VSIC Code 5 Digits

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:

# Number of digits from nganh_kd: 

dn2015 %>% mutate(vsic_n_digits = str_count(nganh_kd)) -> dn2015

dn2015 %>% 
  group_by(vsic_n_digits) %>% 
  summarise(tansuat = n()) %>% 
  arrange(-tansuat)

## # A tibble: 3 × 2
##   vsic_n_digits tansuat
##           <int>   <int>
## 1             5  440397
## 2             4    9502
## 3            NA       2

Kết quả này chỉ ra có 9502 quan sát mà mã ngành (VSIC Code) có 4 kí tự. Chúng ta có thể khôi phục nguyên trạng VISIC Code như sau:

# Add 0 for 4-digits sector codes: 

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

# Remove missing: 

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

Chúng ta có thể viết một hàm mà xử lí đồng thời cả hai việc: (1) loại bỏ các observations mà tax code không hợp lệ, và (2) khôi phục - loại bỏ các quan sát mà VSIC codes không hợp lí như sau:

# Function for processing/cleaning data: 

cleaning_data_tax_vsic_code <- function(your_ves_data) {
  
  # your_ves_data <- dn2015
  
  # Convert to text: 
  
  your_ves_data %>% mutate(nganh_kd = as.character(nganh_kd)) -> ves_data 
  
  # Remove missing at nganh_kd: 
  
  ves_data %>% filter(!is.na(nganh_kd)) -> ves_data
  
  # Remove missing at ma_thue: 
  
  ves_data %>% filter(!is.na(ma_thue)) -> ves_data
  
  #---------------
  # VSIC 5 digits
  #---------------
  
  ves_data %>% 
    mutate(nganh_kd = case_when(str_count(nganh_kd) == 4 ~ str_c("0", nganh_kd), TRUE ~ nganh_kd)) %>% 
    filter(str_count(nganh_kd) == 5) -> ves_data
  
  #----------------------
  # Tax code processing
  #----------------------
  
  # Remove missing: 
  
  ves_data %>% filter(!is.na(ma_thue)) -> ves_data
  
  # Remove duplications: 
  
  ves_data %>% filter(!duplicated(ma_thue)) -> ves_data
  
  # Remove cases != 10 digits: 
  
  ves_data %>% filter(str_count(ma_thue) == 10) -> ves_data

  return(ves_data)
  
}

Với hàm đã có ở trên chúng ta có thể sử dụng để thực hiện việc xử lí tax-vsic code cho bất kì bộ dữ liệu nào. Chẳng hạn của năm 2014:

# Load VES data - 2014: 

haven::read_dta("F:\\VES_from_MaiVu_FTU\\Stata_2014\\dn2014.dta") -> dn2014

# Clean data for VES 2014 data: 

dn2014 %>% cleaning_data_tax_vsic_code() -> dn2014_processed

Problem 3: VSIC Code Level 3

Cũng theo quyết định 27/2018/QĐ-TTg thì mã ngành cấp 3 và cấp 4 có sự khác biệt. Do vậy chúng ta cần xử lí (hay đồng nhất) về các mã cấp 3 và cấp 4 này. Trước hết load dữ liệu chuyển đổi VSIC code:

# Load data (download from https://dangkykinhdoanh.gov.vn/vn/Pages/NganhNghe.aspx): 

readxl::read_xls("C://Users//Admin//Documents//Bảng chuyển đổi VSIC 2018 - VSIC 2007.xls") -> vsic_converted

vsic_converted %>% slice(-c(1:4)) -> vsic_converted

names(vsic_converted) <- c("vsic_new_c3", "vsic_new_c4", "sector_name_new", 
                           "vsic_old_c3", "vsic_old_c4", "sector_name_old")

vsic_converted %>% 
  filter(vsic_new_c3 != vsic_old_c3) %>% 
  select(vsic_new_c3, vsic_old_c3, sector_name_new, sector_name_old) -> vsic_diff_c3

# Show: 

vsic_diff_c3 %>% 
  select(vsic_new_c3, vsic_old_c3)

## # A tibble: 3 × 2
##   vsic_new_c3 vsic_old_c3
##   <chr>       <chr>      
## 1 023         022        
## 2 139         132        
## 3 799         792

Điều này có nghĩa là chúng ta cần phải, ví dụ, convert mã ngành (cấp 3) 022 về 023 (tên ngành: Khai thác, thu nhặt lâm sản khác trừ gỗ). Dưới đây là R codes xử lí:

# Extract old VSIC codes: 

dn2014_processed %>% 
  mutate(code_l2 = str_sub(nganh_kd, start = 1, end = 2), 
         code_l3 = str_sub(nganh_kd, start = 1, end = 3), 
         code_l4 = str_sub(nganh_kd, start = 1, end = 4), 
         code_5_end = str_sub(nganh_kd, start = 5, end = 5)) -> dn2014_processed

#-------------------------------------
#      Adjust VSIC code level 3
#-------------------------------------

# Solution 1 (use case_when() funtion): 

vsic_diff_c3$vsic_old_c3 -> vsic_old_c3

vsic_diff_c3$vsic_new_c3 -> vsic_new_c3

dn2014_processed %>% 
  mutate(code_l3_adj = case_when(code_l3 == vsic_old_c3[1] ~ vsic_new_c3[1], 
                                 code_l3 == vsic_old_c3[2] ~ vsic_new_c3[2], 
                                 code_l3 == vsic_old_c3[3] ~ vsic_new_c3[3], 
                                 TRUE ~ code_l3)) -> dn2014_processed_c1

# Check output: 

dn2014_processed_c1 %>% 
  filter(code_l3 %in% vsic_old_c3) %>% 
  select(nganh_kd, contains("code")) %>% 
  sample_n(6)

## # A tibble: 6 × 6
##   nganh_kd code_l2 code_l3 code_l4 code_5_end code_l3_adj
##   <chr>    <chr>   <chr>   <chr>   <chr>      <chr>      
## 1 13220    13      132     1322    0          139        
## 2 13290    13      132     1329    0          139        
## 3 13220    13      132     1322    0          139        
## 4 13220    13      132     1322    0          139        
## 5 13220    13      132     1322    0          139        
## 6 13220    13      132     1322    0          139

# Solution 2: 

full_join(dn2014_processed, 
          vsic_diff_c3 %>% select(vsic_new_c3, code_l3 = vsic_old_c3),
          by = c("code_l3")) -> dn2014_processed_c2

# Compare: 

dn2014_processed_c2 %>% 
  select(nganh_kd, contains("code"), vsic_new_c3) %>% 
  head()

## # A tibble: 6 × 6
##   nganh_kd code_l2 code_l3 code_l4 code_5_end vsic_new_c3
##   <chr>    <chr>   <chr>   <chr>   <chr>      <chr>      
## 1 42101    42      421     4210    1          <NA>       
## 2 23941    23      239     2394    1          <NA>       
## 3 41000    41      410     4100    0          <NA>       
## 4 33140    33      331     3314    0          <NA>       
## 5 42200    42      422     4220    0          <NA>       
## 6 21001    21      210     2100    1          <NA>

dn2014_processed_c2 %>% 
  select(nganh_kd, contains("code"), vsic_new_c3) %>% 
  filter(!is.na(vsic_new_c3)) %>% 
  head()

## # A tibble: 6 × 6
##   nganh_kd code_l2 code_l3 code_l4 code_5_end vsic_new_c3
##   <chr>    <chr>   <chr>   <chr>   <chr>      <chr>      
## 1 13220    13      132     1322    0          139        
## 2 13220    13      132     1322    0          139        
## 3 13220    13      132     1322    0          139        
## 4 13220    13      132     1322    0          139        
## 5 13220    13      132     1322    0          139        
## 6 13220    13      132     1322    0          139

# Process for missing: 

dn2014_processed_c2 %>% 
  mutate(code_l3_adj = case_when(is.na(vsic_new_c3) ~ code_l3, TRUE ~ vsic_new_c3)) -> dn2014_processed_c2

# Final output: 

dn2014_processed_c2 %>% 
  select(nganh_kd, contains("code"), vsic_new_c3) %>% 
  head()

## # A tibble: 6 × 7
##   nganh_kd code_l2 code_l3 code_l4 code_5_end code_l3_adj vsic_new_c3
##   <chr>    <chr>   <chr>   <chr>   <chr>      <chr>       <chr>      
## 1 42101    42      421     4210    1          421         <NA>       
## 2 23941    23      239     2394    1          239         <NA>       
## 3 41000    41      410     4100    0          410         <NA>       
## 4 33140    33      331     3314    0          331         <NA>       
## 5 42200    42      422     4220    0          422         <NA>       
## 6 21001    21      210     2100    1          210         <NA>

# Check again: 

dn2014_processed_c2 %>% 
  filter(code_l3 %in% vsic_old_c3) %>% 
  select(nganh_kd, contains("code"), vsic_new_c3) %>% 
  head()

## # A tibble: 6 × 7
##   nganh_kd code_l2 code_l3 code_l4 code_5_end code_l3_adj vsic_new_c3
##   <chr>    <chr>   <chr>   <chr>   <chr>      <chr>       <chr>      
## 1 13220    13      132     1322    0          139         139        
## 2 13220    13      132     1322    0          139         139        
## 3 13220    13      132     1322    0          139         139        
## 4 13220    13      132     1322    0          139         139        
## 5 13220    13      132     1322    0          139         139        
## 6 13220    13      132     1322    0          139         139

Problem 4: VSIC Code Level 4

Hiệu chỉnh VSIC cấp 4 nảy sinh một vấn đề: đó là sự không rõ ràng và tường minh. Một ví dụ là mã ngành cấp 4 (theo hệ thống cũ) 0130 (Nhân và chăm sóc cây giống nông nghiệp) sẽ được convert sang mã 0131 (Nhân và chăm sóc cây giống hàng năm) hoặc 0132 (Nhân và chăm sóc cây giống lâu năm). Tuy nhiên luật chuyển đổi không tuyên bố rõ tình huống nào thì 0130 sẽ chuyển đổi sang mã 0131/0132 như ta có thể thấy:

vsic_converted %>% 
  filter(vsic_new_c4 != vsic_old_c4) -> vsic_diff_c4

library(knitr)

vsic_diff_c4 %>% 
  select(-1, -4) %>% 
  slice(1:2) %>%  
  kable()

vsic_new_c4	sector_name_new	vsic_old_c4	sector_name_old
0131	Nhân và chăm sóc cây giống hàng năm	0130	Nhân và chăm sóc cây giống nông nghiệp
0132	Nhân và chăm sóc cây giống lâu năm	0130	Nhân và chăm sóc cây giống nông nghiệp

Trong khi chưa rõ ràng về luật chuyển đổi mã cấp 4 cho những cases trùng nhau này, chúng ta tạm chấp nhận giải pháp sau để hiệu chỉnh mã ngành cấp 4:

#-------------------------------------
#      Adjust VSIC code level 4
#-------------------------------------

vsic_diff_c4 %>% 
  filter(!duplicated(vsic_old_c4)) %>% 
  select(vsic_new_c4, code_l4 = vsic_old_c4) -> convert_table_c4

full_join(dn2014_processed_c2, convert_table_c4, by = "code_l4") -> dn2014_processed_c2

# Our data: 

dn2014_processed_c2 %>% 
  select(nganh_kd, contains("code"), vsic_new_c4) %>% 
  head()

## # A tibble: 6 × 7
##   nganh_kd code_l2 code_l3 code_l4 code_5_end code_l3_adj vsic_new_c4
##   <chr>    <chr>   <chr>   <chr>   <chr>      <chr>       <chr>      
## 1 42101    42      421     4210    1          421         4211       
## 2 23941    23      239     2394    1          239         <NA>       
## 3 41000    41      410     4100    0          410         4101       
## 4 33140    33      331     3314    0          331         <NA>       
## 5 42200    42      422     4220    0          422         4221       
## 6 21001    21      210     2100    1          210         <NA>

# Adjust missing: 

dn2014_processed_c2 %>% 
  mutate(code_l4_adj = case_when(is.na(vsic_new_c4) ~ code_l4, TRUE ~ vsic_new_c4)) -> dn2014_processed_c2

# Compare: 

dn2014_processed_c2 %>% 
  select(nganh_kd, contains("code"), vsic_new_c4) %>% 
  head()

## # A tibble: 6 × 8
##   nganh_kd code_l2 code_l3 code_l4 code_5_end code_l3_adj code_l4_adj vsic_new…¹
##   <chr>    <chr>   <chr>   <chr>   <chr>      <chr>       <chr>       <chr>     
## 1 42101    42      421     4210    1          421         4211        4211      
## 2 23941    23      239     2394    1          239         2394        <NA>      
## 3 41000    41      410     4100    0          410         4101        4101      
## 4 33140    33      331     3314    0          331         3314        <NA>      
## 5 42200    42      422     4220    0          422         4221        4221      
## 6 21001    21      210     2100    1          210         2100        <NA>      
## # … with abbreviated variable name ¹​vsic_new_c4

Đến đây chúng ta đã hiệu chỉnh xong mã ngành cấp 4. Chúng ta có thể khôi phục lại mã ngành đầy đủ theo quyết định 27/2018/QĐ-TTg cho dữ liệu của năm 2014 rồi lưu lại để sử dụng hoặc gửi cho đồng nghiệp như sau:

# Create new VSIC code: 

dn2014_processed_c2 %>% 
  mutate(vsic_code = str_c(code_l4_adj, code_5_end)) -> dn2014_processed_c2

# Check again: 

dn2014_processed_c2 %>% 
  select(nganh_kd, contains("code"), contains("vsic")) %>% 
  head() %>% 
  kable()

nganh_kd	code_l2	code_l3	code_l4	code_5_end	code_l3_adj	code_l4_adj	vsic_code	vsic_new_c3	vsic_new_c4
42101	42	421	4210	1	421	4211	42111	NA	4211
41000	41	410	4100	0	410	4101	41010	NA	4101
42200	42	422	4220	0	422	4221	42210	NA	4221
41000	41	410	4100	0	410	4101	41010	NA	4101
42102	42	421	4210	2	421	4211	42112	NA	4211
41000	41	410	4100	0	410	4101	41010	NA	4101

# Save data: 

haven::write_dta(dn2014_processed_c2, "dn2014_cleaned.dta")

Problem 5: VSIC Code Level 1

Các paper thường chỉ sử dụng đến mã ngành cấp 1 hoặc cấp 2. Trước hết lấy dữ liệu về mã ngành cấp 1 (và cả cấp 2) như sau:

# --------------------------------------------------------------------------------
# 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

# Save for using laler: 

writexl::write_xlsx(df_sector_code, "df_sector_code.xlsx")

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 = "Some Level-1 Sectors")

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 1 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

Tương tự:

df_sector_code_filled %>% 
  mutate(code_level2 = case_when(str_count(code_level2) == 0 ~ NA_character_, TRUE ~ code_level2)) %>% 
  tidyr::fill(code_level2, .direction = "down") -> df_sector_code_filled

Với data đã xử lí chúng ta có thể, ví dụ, báo cáo một số thống kê cho 21 ngành cấp 1:

df_sector_code %>% 
  select(code_level1, sector_name_latin) %>% 
  filter(str_count(code_level1) != 0) -> df_sector_l1

Rồi join các data sets:

full_join(df_sector_l1, 
          df_sector_code_filled %>% select(code_level1, code_level2)) -> df_sector_l1

df_sector_l1 %>% 
  rename(code_l2 = code_level2) %>% 
  filter(!duplicated(code_l2)) -> df_sector_l1


inner_join(dn2014_processed_c2, 
           df_sector_l1, 
           by = c("code_l2")) -> dn2014_processed_c2_updated 

Với data cuối cùng này chúng ta có thể, ví dụ, tạo ra báo cáo như sau:

dn2014_processed_c2_updated %>% 
  filter(!is.na(kqkd1)) %>% 
  group_by(sector_name_latin) %>% 
  summarise(Min = min(kqkd1), 
            Max = max(kqkd1), 
            SD = sd(kqkd1), 
            Median = median(kqkd1), 
            N_obs = n()) %>% 
  arrange(-N_obs) -> df_report

df_report %>% 
  mutate(sector_name_short = str_sub(sector_name_latin, 1, 20)) %>% 
  select(-sector_name_latin) %>% 
  select(sector_name_short, everything()) %>% 
  kable()

sector_name_short	Min	Max	SD	Median	N_obs
XAY DUNG	-514	7085110.0	104805.13	1692.0	41516
CONG NGHIEP CHE BIEN	0	6930426.0	257591.51	2517.0	3008
GIAO DUC VA DAO TAO	0	1647168.0	67346.63	364.0	1379
VAN TAI KHO BAI	-220	3179655.0	131597.57	1209.5	1008
SAN XUAT VA PHAN PHO	0	52630921.0	2991102.80	3159.0	691
HOAT DONG HANH CHINH	0	435217.0	23495.42	348.5	584
HOAT DONG CHUYEN MON	0	179989.0	14849.89	210.0	247
BAN BUON VA BAN LE;	0	245064.8	30109.35	1433.5	80
HOAT DONG TAI CHINH,	0	6587334.0	1239477.31	23633.0	52

Task for You

Bằng các bước xử lí như trên chúng ta có thể hiệu chỉnh mã ngành cấp 3, cấp 4 và tạo ra mã ngành VSIC mới cho bất kì bộ dữ liệu nào. Tuy nhiên sẽ thuận lợi hơn khi viết thành hàm để còn tái sử dụng nhằm đảm bảo nguyên lí DRY. Hãy viết hai hàm lần lượt có tên là process_vsic3() và process_vsic4() để hiệu chỉnh mã cấp 3 và cấp 4 nếu input của hàm là một data frame.

Dưới đây là phương án đề xuất để viết hai hàm có mô tả ở trên:

# Function corrects VSIC Level 3: 

process_vsic3 <- function(your_df) {
  
  your_df %>% 
    mutate(code_l2 = str_sub(nganh_kd, start = 1, end = 2), 
           code_l3 = str_sub(nganh_kd, start = 1, end = 3), 
           code_l4 = str_sub(nganh_kd, start = 1, end = 4), 
           code_5_end = str_sub(nganh_kd, start = 5, end = 5)) -> your_df
  
  full_join(your_df, 
            vsic_diff_c3 %>% select(vsic_new_c3, code_l3 = vsic_old_c3),
            by = c("code_l3")) -> your_df
  
  your_df %>% 
    mutate(code_l3_adj = case_when(is.na(vsic_new_c3) ~ code_l3, TRUE ~ vsic_new_c3)) -> final_df
  
  return(final_df)
}

# Function corrects VSIC Level 4:


process_vsic4 <- function(your_df) {
  
  your_df %>% 
    mutate(code_l2 = str_sub(nganh_kd, start = 1, end = 2), 
           code_l3 = str_sub(nganh_kd, start = 1, end = 3), 
           code_l4 = str_sub(nganh_kd, start = 1, end = 4), 
           code_5_end = str_sub(nganh_kd, start = 5, end = 5)) -> your_df
  
  full_join(your_df, convert_table_c4, by = "code_l4") -> your_df
  
  your_df %>% 
    mutate(code_l4_adj = case_when(is.na(vsic_new_c4) ~ code_l4, TRUE ~ vsic_new_c4)) -> final_df
  
  your_df %>% 
    mutate(code_l4_adj = case_when(is.na(vsic_new_c4) ~ code_l4, TRUE ~ vsic_new_c4)) %>% 
    mutate(nganh_kd_adj = str_c(code_l4_adj, code_5_end)) %>% 
    mutate(code_l2_adj = str_sub(nganh_kd_adj, start = 1, end = 2)) -> final_df
  
  return(final_df)
  
}

Với hai hàm này chúng ta có thể xử lí dữ liệu VES cho bất cứ năm nào. Ví dụ:

# Load VES 2015: 

haven::read_dta("F:\\VES_from_MaiVu_FTU\\Stata_2015\\dn2015.dta") -> dn2015

# Process VES 2015 data: 

dn2015 %>% 
  cleaning_data_tax_vsic_code() %>% 
  process_vsic3() %>% 
  process_vsic4() -> dn2015_after_processing

Tương tự chúng ta viết hàm join tên ngành cấp 1:

# Funtions join sector name: 

add_sector_name_at_level_1 <- function(your_df) {
  
  inner_join(your_df, 
             df_sector_l1, 
             by = c("code_l2_adj" = "code_l2")) -> final_df
  
  return(final_df)
  
}

Rồi sử dụng hàm này:

dn2015_after_processing %>% add_sector_name_at_level_1() -> dn2015_after_processing

Tính toán một số thống kê cho kqkd1 theo ngành cấp 1 cho năm 2015:

dn2015_after_processing %>% 
  filter(!is.na(kqkd1)) %>% 
  group_by(sector_name_latin) %>% 
  summarise(Min = min(kqkd1), 
            Max = max(kqkd1), 
            SD = sd(kqkd1), 
            Median = median(kqkd1), 
            N_obs = n()) %>% 
  arrange(-N_obs) %>% 
  mutate(sector_name_short = str_sub(sector_name_latin, 1, 20)) %>% 
  select(-sector_name_latin) %>% 
  select(sector_name_short, everything()) %>% 
  kable()

sector_name_short	Min	Max	SD	Median	N_obs
BAN BUON VA BAN LE;	0	105559749	395966.60	2102.00	169612
CONG NGHIEP CHE BIEN	0	375812792	2153560.26	2816.00	66510
XAY DUNG	-26	11966123	133276.12	1717.00	58651
HOAT DONG CHUYEN MON	-95	3404914	39534.42	435.00	36787
VAN TAI KHO BAI	0	10965323	130059.25	2175.00	26425
DICH VU LUU TRU VA A	0	4694252	62100.38	586.00	16056
HOAT DONG HANH CHINH	0	4815353	68874.28	545.00	15457
THONG TIN VA TRUYEN	0	10098702	145922.91	231.15	9448
HOAT DONG KINH DOANH	0	7292986	169680.21	290.50	8624
NONG NGHIEP, LAM NGH	0	13207539	190797.08	944.00	6300
GIAO DUC VA DAO TAO	0	1136408	35180.19	106.00	5283
HOAT DONG TAI CHINH,	-292	22942233	785044.73	1699.00	3175
HOAT DONG DICH VU KH	0	366965	13330.10	165.00	3155
KHAI KHOANG	0	57474075	1808576.59	3378.50	2396
NGHE THUAT, VUI CHOI	0	5615034	275200.03	153.00	2374
CUNG CAP NUOC; HOAT	0	2826496	118502.49	1002.00	1702
SAN XUAT VA PHAN PHO	0	77951360	3103259.73	2445.00	1495
Y TE VA HOAT DONG TR	0	969653	59694.20	753.90	1420

Prepare Panel Data

Giả sử chúng ta áp dụng Panel Data cho các năm 2014, 2015, 2016 để ước lượng mô hình sau:

\[\begin{equation} Y_{it}=\beta_{1}+\beta_{2}X_{it}+e_{it} \label{eq:panelgeneq15} \end{equation}\]

Trong đó Y là kqkd1 và X là ts11. Trước hết xử lí dữ liệu theo chuỗi như chúng ta đã biết cho năm 2014 và 2016:

# For 2014: 

haven::read_dta("F:\\VES_from_MaiVu_FTU\\Stata_2014\\dn2014.dta") %>% 
  cleaning_data_tax_vsic_code() %>% 
  process_vsic3() %>% 
  process_vsic4() %>% 
  add_sector_name_at_level_1() -> dn2014_after_processing

# For 2016: 

haven::read_dta("F:\\VES_from_MaiVu_FTU\\New data (2016-2020)\\Stata_2016\\dn2016.dta") %>% 
  cleaning_data_tax_vsic_code() %>% 
  process_vsic3() %>% 
  process_vsic4() %>% 
  add_sector_name_at_level_1() -> dn2016_after_processing

Tạo Panel Data:

dn2014_after_processing %>% 
  select(nganh_kd, ma_thue, sector_name_latin, kqkd1, ts11) %>% 
  mutate(year = 2014) %>% 
  na.omit() -> df2014

dn2015_after_processing %>% 
  select(nganh_kd, ma_thue, sector_name_latin, kqkd1, ts11) %>% 
  mutate(year = 2015) %>% 
  na.omit() -> df2015

dn2016_after_processing %>% 
  select(nganh_kd, ma_thue, sector_name_latin, kqkd1, ts11) %>% 
  mutate(year = 2016) %>% 
  na.omit() -> df2016

df2014 %>% 
  bind_rows(df2015) %>% 
  bind_rows(df2016) -> my_panel

# Prepare panel data: 

library(plm)

pdata.frame(my_panel, 
            index = c("ma_thue","year"), 
            row.names = FALSE, 
            drop.index = FALSE) -> panel_data_unbalanced

Tính toán một số thống kê thường được trình bày trong các paper:

library(tidyr)

my_panel %>% 
  select(kqkd1, ts11, year) %>% 
  pivot_longer(cols = c("kqkd1", "ts11"), names_to = "variable", values_to = "value") %>% 
  group_by(variable) %>% 
  summarise(Mean = mean(value), 
            Median = median(value), 
            SD = sd(value), 
            Min = min(value), 
            Max = max(value), 
            N_obs = n()) %>% 
  kable()

variable	Mean	Median	SD	Min	Max	N_obs
kqkd1	33309.63	1540	934370.5	-575.0	457221004	1281370
ts11	38280.46	4058	974257.8	-48923.9	350857811	1281370

Với dữ liệu đã có, chúng ta có thể chạy các mô hình phân tích cho Panel Data. Chẳng hạn ba cách tiếp cận được chọn là Pooled, Fixed Effects và Random Effects sử dụng ước lượng Amemiya:

#========================
# Some Panel Data Models
#========================

# Pooled Model: 
panel_pooled <- plm(kqkd1 ~ ts11, 
                    data = panel_data_unbalanced, 
                    model = "pooling") 

# Fixed Effects Model: 
panel_fe <- plm(kqkd1 ~ ts11, 
                data = panel_data_unbalanced, 
                model = "within")  

# Random Effects Models using Amemiya estimators (1971): 
panel_re <- plm(kqkd1 ~ ts11, 
                data = panel_data_unbalanced, 
                model = "random", 
                random.method = "amemiya")

Rồi so sánh các ước lượng từ ba mô hình này:

# Compare results: 

library(stargazer)

stargazer(panel_pooled, 
          panel_fe, 
          panel_re, 
          title = "Compare Panel Data Regression Results - Unbalanced",
          column.labels = c("Pooled","Fixed Effects", "Random Effects"), 
          type = "text", 
          align = TRUE)

## 
## Compare Panel Data Regression Results - Unbalanced
## ===========================================================================================
##                                           Dependent variable:                              
##              ------------------------------------------------------------------------------
##                                                  kqkd1                                     
##                           Pooled                      Fixed Effects          Random Effects
##                            (1)                             (2)                    (3)      
## -------------------------------------------------------------------------------------------
## ts11                     0.391***                        0.208***               0.240***   
##                          (0.001)                         (0.001)                (0.001)    
##                                                                                            
## Constant              18,342.110***                                          20,272.890*** 
##                         (754.302)                                             (1,086.506)  
##                                                                                            
## -------------------------------------------------------------------------------------------
## Observations            1,281,370                       1,281,370              1,281,370   
## R2                        0.166                           0.109                  0.120     
## Adjusted R2               0.166                           -0.568                 0.120     
## F Statistic  255,430.000*** (df = 1; 1281368) 89,487.630*** (df = 1; 727852) 174,088.900***
## ===========================================================================================
## Note:                                                           *p<0.1; **p<0.05; ***p<0.01

Balanced Panel Data

Để chuẩn bị Balanced Panel Data chúng ta làm như sau:

my_panel %>% 
  group_by(ma_thue) %>% 
  summarise(tansuat = n()) %>% 
  filter(tansuat == 3) %>% 
  pull(ma_thue) -> firms_3_times

pdata.frame(my_panel %>% filter(ma_thue %in% firms_3_times), 
            index = c("ma_thue","year"), 
            row.names = FALSE, 
            drop.index = FALSE) -> panel_data_balanced

Rồi thực hiện lại các bước phân tích như trên:

# Pooled Model: 
panel_pooled_b <- plm(kqkd1 ~ ts11, 
                      data = panel_data_balanced, 
                      model = "pooling") 

# Fixed Effects Model: 
panel_fe_b <- plm(kqkd1 ~ ts11, 
                  data = panel_data_balanced, 
                  model = "within")  

# Random Effects Models using Amemiya estimators (1971): 
panel_re_b <- plm(kqkd1 ~ ts11, 
                  data = panel_data_balanced, 
                  model = "random", 
                  random.method = "amemiya")

stargazer(panel_pooled_b, 
          panel_fe_b, 
          panel_re_b, 
          title = "Compare Panel Data Regression Results - Balanced",
          column.labels = c("Pooled","Fixed Effects", "Random Effects"), 
          type = "text", 
          align = TRUE)

## 
## Compare Panel Data Regression Results - Balanced
## ==========================================================================================
##                                           Dependent variable:                             
##              -----------------------------------------------------------------------------
##                                                  kqkd1                                    
##                          Pooled                      Fixed Effects          Random Effects
##                            (1)                            (2)                    (3)      
## ------------------------------------------------------------------------------------------
## ts11                    0.549***                        0.253***               0.308***   
##                          (0.001)                        (0.001)                (0.001)    
##                                                                                           
## Constant              17,471.040***                                         28,241.140*** 
##                        (1,009.208)                                           (1,747.594)  
##                                                                                           
## ------------------------------------------------------------------------------------------
## Observations             917,379                        917,379                917,379    
## R2                        0.213                          0.114                  0.132     
## Adjusted R2               0.213                          -0.329                 0.132     
## F Statistic  248,668.000*** (df = 1; 917377) 78,610.910*** (df = 1; 611585) 139,304.500***
## ==========================================================================================
## Note:                                                          *p<0.1; **p<0.05; ***p<0.01

Final Notes

Sử dụng VES (cùng với VHLSS) cho research paper là một công việc rất thách thức, đặt biệt là xử lí đồng thời nhiều bộ dữ liệu cùng lúc trong một khoảng thời gian dài. Về chỉ riêng điểm này thôi thì R tốt hơn Stata nhiều.

Appendix

#======================================
#      R codes for explanation
#======================================

# Clear R environment: 

rm(list = ls())

# Create fake data sets: 

ves_data <- data.frame(ma_thue = c("564", "213", "981", "", "213", "83", "745"), 
                       nganh_kd = c("64910", "13251", "2145", "17941", "13251", "27513", "1494")) # VES data


vsic_update <- data.frame(vsic3_new = c("139", "017"), 
                          vsic3_old = c("132", "014")) # Data for correcting. 

# Show our data: 

ves_data

vsic_update

#-------------------------
#  Duplicated Tax Code
#-------------------------

library(dplyr)

ves_data %>% 
  group_by(ma_thue) %>% 
  summarise(tansuat = n()) 

ves_data %>% 
  filter(!duplicated(ma_thue)) -> ves_not_dup

ves_not_dup

#----------------------
#  Taxt Code 3-digits
#----------------------

library(stringr)

ves_not_dup %>% 
  mutate(slkt_tax = str_count(ma_thue)) %>% 
  filter(slkt_tax == 3) # Solution 1. 

ves_not_dup %>% 
  filter(str_count(ma_thue) == 3) # Solution 2. 

#----------------------------
#  Add zero for sector code
#----------------------------

ves_not_dup %>% 
  filter(str_count(ma_thue) == 3) %>% 
  mutate(slkt_code = str_count(nganh_kd)) %>% 
  mutate(nganh_kd_adj = case_when(slkt_code == 4 ~ str_c("0", nganh_kd), 
                                  TRUE ~ nganh_kd))

ves_not_dup %>% 
  filter(str_count(ma_thue) == 3) %>% 
  mutate(slkt_code = str_count(nganh_kd)) %>% 
  mutate(nganh_kd = case_when(slkt_code == 4 ~ str_c("0", nganh_kd), 
                              TRUE ~ nganh_kd)) -> ves_cleaned

#-----------------------
#  Adjust sector code
#-----------------------

# About str_sub() function: 

ves_cleaned %>% pull(nganh_kd) -> nganh_vec

str_sub(nganh_vec, start = 1, end = 3)

str_sub(nganh_vec, start = 5, end = 5)

ves_cleaned %>% 
  mutate(code_l3 = str_sub(nganh_kd, start = 1, end = 3), 
         code_5th = str_sub(nganh_kd, start = 5, end = 5))

ves_cleaned %>% 
  mutate(code_l3 = str_sub(nganh_kd, start = 1, end = 3), 
         code_5th = str_sub(nganh_kd, start = 5, end = 5)) -> ves_cleaned

# Compare: 

vsic_update

vsic_update %>% pull(vsic3_old) -> old_code3

vsic_update$vsic3_new -> new_code3

# ----- Soluion 1 for correcting sector code -----

ves_cleaned %>% 
  mutate(code_l3_adj = case_when(code_l3 == old_code3[1] ~ new_code3[1], 
                                 code_l3 == old_code3[2] ~ new_code3[2], 
                                 TRUE ~ code_l3))

# ----  Soluion 2 for correcting sector code -----

full_join(ves_cleaned, vsic_update, by = c("code_l3" = "vsic3_old"))

inner_join(ves_cleaned, vsic_update, by = c("code_l3" = "vsic3_old"))

full_join(ves_cleaned, vsic_update, by = c("code_l3" = "vsic3_old")) -> ves_cleaned_new

ves_cleaned_new %>% 
  mutate(code_l3_adj = case_when(is.na(vsic3_new) ~ code_l3, 
                                 TRUE ~ vsic3_new))