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:
::read_dta("F:\\VES_from_MaiVu_FTU\\Stata_2015\\dn2015.dta") -> dn2015
haven
# 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ì:
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:
%>% filter(!duplicated(ma_thue)) -> dn2015
dn2015
#--------- Remove firms with tax code != 10 digits ----------------
library(stringr)
%>% filter(str_count(ma_thue) == 10) -> dn2015 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:
# Number of digits from nganh_kd:
%>% mutate(vsic_n_digits = str_count(nganh_kd)) -> dn2015
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:
%>% filter(!is.na(nganh_kd)) -> dn2015 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:
<- function(your_ves_data) {
cleaning_data_tax_vsic_code
# your_ves_data <- dn2015
# Convert to text:
%>% mutate(nganh_kd = as.character(nganh_kd)) -> ves_data
your_ves_data
# Remove missing at nganh_kd:
%>% filter(!is.na(nganh_kd)) -> ves_data
ves_data
# Remove missing at ma_thue:
%>% filter(!is.na(ma_thue)) -> ves_data
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:
%>% filter(!is.na(ma_thue)) -> ves_data
ves_data
# Remove duplications:
%>% filter(!duplicated(ma_thue)) -> ves_data
ves_data
# Remove cases != 10 digits:
%>% filter(str_count(ma_thue) == 10) -> ves_data
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:
::read_dta("F:\\VES_from_MaiVu_FTU\\Stata_2014\\dn2014.dta") -> dn2014
haven
# Clean data for VES 2014 data:
%>% cleaning_data_tax_vsic_code() -> dn2014_processed dn2014
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):
::read_xls("C://Users//Admin//Documents//Bảng chuyển đổi VSIC 2018 - VSIC 2007.xls") -> vsic_converted
readxl
%>% slice(-c(1:4)) -> vsic_converted
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_old_c3 -> vsic_old_c3
vsic_diff_c3
$vsic_new_c3 -> vsic_new_c3
vsic_diff_c3
%>%
dn2014_processed mutate(code_l3_adj = case_when(code_l3 == vsic_old_c3[1] ~ vsic_new_c3[1],
== vsic_old_c3[2] ~ vsic_new_c3[2],
code_l3 == vsic_old_c3[3] ~ vsic_new_c3[3],
code_l3 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,
%>% select(vsic_new_c3, code_l3 = vsic_old_c3),
vsic_diff_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
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:
::write_dta(dn2014_processed_c2, "dn2014_cleaned.dta") haven
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)
<- "https://dangkykinhdoanh.gov.vn/vn/Pages/NganhNghe.aspx"
link
%>%
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:
::write_xlsx(df_sector_code, "df_sector_code.xlsx") writexl
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")
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)) %>%
::fill(code_level1, .direction = "down") %>%
tidyrfilter(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)) %>%
::fill(code_level2, .direction = "down") -> df_sector_code_filled tidyr
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,
%>% select(code_level1, code_level2)) -> df_sector_l1
df_sector_code_filled
%>%
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 |
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:
<- function(your_df) {
process_vsic3
%>%
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,
%>% select(vsic_new_c3, code_l3 = vsic_old_c3),
vsic_diff_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:
<- function(your_df) {
process_vsic4
%>%
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:
::read_dta("F:\\VES_from_MaiVu_FTU\\Stata_2015\\dn2015.dta") -> dn2015
haven
# 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:
<- function(your_df) {
add_sector_name_at_level_1
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:
%>% add_sector_name_at_level_1() -> dn2015_after_processing 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 |
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:
::read_dta("F:\\VES_from_MaiVu_FTU\\Stata_2014\\dn2014.dta") %>%
havencleaning_data_tax_vsic_code() %>%
process_vsic3() %>%
process_vsic4() %>%
add_sector_name_at_level_1() -> dn2014_after_processing
# For 2016:
::read_dta("F:\\VES_from_MaiVu_FTU\\New data (2016-2020)\\Stata_2016\\dn2016.dta") %>%
havencleaning_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:
<- plm(kqkd1 ~ ts11,
panel_pooled data = panel_data_unbalanced,
model = "pooling")
# Fixed Effects Model:
<- plm(kqkd1 ~ ts11,
panel_fe data = panel_data_unbalanced,
model = "within")
# Random Effects Models using Amemiya estimators (1971):
<- plm(kqkd1 ~ ts11,
panel_re 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
Để 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:
<- plm(kqkd1 ~ ts11,
panel_pooled_b data = panel_data_balanced,
model = "pooling")
# Fixed Effects Model:
<- plm(kqkd1 ~ ts11,
panel_fe_b data = panel_data_balanced,
model = "within")
# Random Effects Models using Amemiya estimators (1971):
<- plm(kqkd1 ~ ts11,
panel_re_b 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
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.
#======================================
# R codes for explanation
#======================================
# Clear R environment:
rm(list = ls())
# Create fake data sets:
<- data.frame(ma_thue = c("564", "213", "981", "", "213", "83", "745"),
ves_data nganh_kd = c("64910", "13251", "2145", "17941", "13251", "27513", "1494")) # VES data
<- data.frame(vsic3_new = c("139", "017"),
vsic_update 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:
%>% pull(nganh_kd) -> nganh_vec
ves_cleaned
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
%>% pull(vsic3_old) -> old_code3
vsic_update
$vsic3_new -> new_code3
vsic_update
# ----- Soluion 1 for correcting sector code -----
%>%
ves_cleaned mutate(code_l3_adj = case_when(code_l3 == old_code3[1] ~ new_code3[1],
== old_code3[2] ~ new_code3[2],
code_l3 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))