Gaps in Household Income From VHLSS 2018

R Daily Graph Series

Introduction

Post trước cho thấy Đà Nẵng là tỉnh có thu nhập trung vị hộ gia đình cao nhất nước (450 triệu). Kết quả đó được tính toán dựa vào biến thunhap từ bộ dữ liệu HO3.dta:

Theo dữ liệu ở bảng này thì hộ ở tinh = 96, huyen = 973, xa = 32239 , diaban = 19, hoso = 1 xuất hiện ở 8 dòng và do vậy giả định rằng mỗi một dòng tương ứng với một thành viên của gia đình thì tổng thu nhập của hộ này sẽ là 169250 + 600 + 17435. Tính toán tương tự theo logic này cho các hộ còn lại.

Post này trình bày một cách tính toán khác căn cứ vào biến thubq - tức là thu nhập bình quân đầu người của hộ. Như vậy theo cách tính này thì chúng ta chỉ cần biết số thành viên trong gia đình của hộ, nhân với thubq rồi nhân tiếp với 12 chúng ta sẽ có tổng thu nhập của hộ theo năm.

Thông tin về số thành viên của hộ gia đình có thể khai thác từ bộ dữ liệu HO1.dta. Cách tính này chúng ta có kết quả như Figure 1 ở trên.

Data Pre-processing

Trước hết đọc dữ liệu vào R:

# Clear R environment: 
rm(list = ls())

# Load some R packages: 
library(dplyr)
library(stringr)
library(stringi)
library(kableExtra) # For presenting table. 
library(haven)

# Import HO3.dta: 

read_dta("D:/VHLSS2018_Household/HO3.dta") -> ho3

# HO1.dta: 

read_dta("D:/VHLSS2018_Household/HO1.dta") -> ho1 

Mỗi một hộ gia đình có ID được hình thành từ mã của 5 cột biến: tỉnh, huyện, xã, địa bàn, hộ số (cái này ban đầu em và nhiều người nữa vẫn hiểu là “hồ sơ” vì chữ không dấu, thực ra nó là “hộ số”). Vấn đề đầu tiên là chúng ta cần chuẩn hóa các mã này một cách thống nhất cho những mục đích sử dụng xa hơn sau này. Chẳng hạn, tỉnh có mã là 1 thì cần phải chuẩn hóa về 01 (tương ứng với Hà Nội). Mã tỉnh code chuẩn sử dụng 2 chữ số và do vậy với các tỉnh là số tự nhiên bé hơn 9 thì chúng ta phải thêm 1 số 0 đằng trước. Tương tự là Huyện sử dụng 3 chữ số đễ mã hóa. Do vậy với huyện mà chỉ sử dụng 1 chữ số thì chúng ta chuẩn hóa bằng cách thêm 2 chữ số 0 đằng trước, với huyện có mã là 2 chữ số thì thêm 1 chữ số 0 đằng trước. Cách thức chuẩn hóa này áp dụng tương tự cho các biến còn lại là xã, địa bàn và hộ số.

Chúng ta viết hàm có tên add_zero() để chuẩn hóa mã hành chính:

#========================
#  Data Pre-processing 
#========================

# Function creates full code by adding zeros: 

add_zero <- function(x) {
  
  tibble(x_text = as.character(x)) %>% 
    mutate(n_digits = str_count(x_text),
           n_max = max(n_digits), 
           delta = n_max - n_digits, 
           pre = strrep("0", times = delta), 
           full_code = str_c(pre, x_text)) %>% 
    pull(full_code) %>% 
    return()
}

Để thuận lợi cho đối chiếu và tra cứu, viết hàm có tên extract_description() để lấy các mô tả về các biến số:

# Function extracts variable description: 

extract_description <- function(df_selected) {
  
  sapply(df_selected, function(x) {attributes(x) %>% .$label}) %>% 
    data.frame() %>% 
    mutate(description = stri_trans_general(`.`, "Latin-ASCII")) -> df_des
  
  df_des %>% 
    mutate(var_name = row.names(df_des)) %>% 
    select(var_name, description) -> df_des
  
  row.names(df_des) <- NULL
  
  return(df_des)
  
}

Ý nghĩa của các biến số thuộc HO3.dta được trình bày ở Table 1 dưới đây:

# Use the function: 

extract_description(df_selected = ho3) -> des_ho3

des_ho3 %>% 
  kbl(caption = "Table 1: Variable Description for HO3.dta", escape = TRUE) %>%
  kable_classic(full_width = TRUE, html_font = "Cambria")

Table 1: Variable Description for HO3.dta

var_name	description
tinh	Tinh
huyen	Huyen
xa	Xa
diaban	Dia ban
hoso	Ho so
thunhap	Thu nhap
thubq	Thu nhap binh quan/nguoi/thang
tongthu	Tong thu
chisxkd	Chi phi SXKD
chikhac	Chi tieu va chi khac

Tạo một biến mới có tên h_code là mã của các hộ gia đình và mã này được hình thành từ mã chuẩn hóa của tỉnh, huyện, xã, địa bàn và hồ sơ:

#------------------------
#  Create Household ID
#------------------------

# Create h_code column for ho3: 

ho3 %>% 
  mutate(tinh_n = add_zero(tinh), 
         huyen_n = add_zero(huyen), 
         xa_n = add_zero(xa), 
         diaban_n = add_zero(diaban), 
         hoso_n = add_zero(hoso)) %>% 
  mutate(h_code = str_c(tinh_n, huyen_n, xa_n, diaban_n, hoso_n)) -> ho3

# For ho1:

ho1 %>% 
  mutate(tinh_n = add_zero(tinh), 
         huyen_n = add_zero(huyen), 
         xa_n = add_zero(xa), 
         diaban_n = add_zero(diaban), 
         hoso_n = add_zero(hoso)) %>% 
  mutate(h_code = str_c(tinh_n, huyen_n, xa_n, diaban_n, hoso_n)) -> ho1

# Remove duplications: 

ho1 %>% filter(!duplicated(h_code)) -> ho1

# Join the two data sets: 

ho3 %>% 
  select(h_code, tinh_n, thubq) %>% 
  filter(!is.na(thubq)) %>% 
  full_join(ho1 %>% select(h_code, tsnguoi), by = "h_code") -> df_monthly_income

# Calculate yearly-household income: 

df_monthly_income %>% 
  mutate(total_year_income = thubq*tsnguoi*12) -> df_monthly_income

#--------------------------
#  Extract province names
#--------------------------

# Extract province info:

ho3$tinh %>% 
  attributes() %>% 
  .$labels %>% data.frame() -> df_province

# Rename for DF: 
names(df_province) <- "province_code"

# Create some columns and relabel for provinces: 

df_province %>% 
  mutate(province_vie = row.names(df_province)) %>% 
  mutate(province_eng = stri_trans_general(province_vie, "Latin-ASCII")) %>% 
  mutate(province_eng = str_replace_all(province_eng, "Tinh |Thanh pho ", "")) %>% 
  mutate(province_eng = str_replace_all(province_eng, " - ", "-")) %>% 
  mutate(province_code = add_zero(province_code)) -> df_province

row.names(df_province) <- NULL

# Join province name: 

df_monthly_income %>% 
  full_join(df_province %>% select(province_eng, province_code), by = c("tinh_n" = "province_code")) -> df_household_income

# Calculate some percentiles by household level: 

df_household_income %>% 
  filter(!is.na(province_eng)) %>% 
  group_by(province_eng) %>% 
  summarise(th25 = quantile(total_year_income, 0.25), 
            th50 = quantile(total_year_income, 0.50), 
            th75 = quantile(total_year_income, 0.75)) %>% 
  mutate_if(is.numeric, function(x) {round(x / 1000, 1)}) %>% 
  ungroup() %>% 
  arrange(th50) %>% 
  mutate(province_eng = factor(province_eng, province_eng)) -> df_thunhap

Nếu chọn mean để tính toán và so sánh mức độ chênh lệnh về thu nhập của hộ gia đình có một điểm yếu là mean bị ảnh hưởng rất lớn bởi Outliers. Do vậy chọn thu nhập trung vị (Median) để đánh giá và so sánh chúng ta có thể khắc phục nhược điểm này.

R codes for Data Visualization

Sau bước xử lí dữ liệu thô và tính toán các thống kê, dưới đây là R Codes tạo ra Figure 1:

#====================
# Data Visualization
#====================

# Load some R packages for Data Visualization: 

library(ggeconodist) # install.packages("ggeconodist", repos = "https://cinc.rud.is")
library(ggplot2)
library(showtext)

# Select Ubuntu Condensed font: 
showtext.auto()
font_add_google(name = "Ubuntu Condensed", family = "ubu")
my_font <- "ubu"


df_thunhap %>% 
  ggplot(aes(x = province_eng)) + 
  geom_econodist(aes(ymin = th25, median = th50, ymax = th75), 
                 median_col = "firebrick", 
                 stat = "identity", 
                 median_point_size = 1.3, 
                 show.legend = TRUE) + 
  coord_flip() +
  theme_econodist() + 
  scale_y_continuous(expand = c(0, 0), limits = range(0, 350), breaks = seq(0, 350, 50), position = "right") + 
  labs(title = "Figure 1: Gaps in Household Income (millions VND) by Province, 2018", 
       caption = "Data Source: VHLSS 2018, GSO") +  
  theme(plot.margin = unit(c(0.3, 1, 0.3, 0.5), "cm")) + 
  theme(axis.title.y = element_blank()) + 
  theme(axis.text.y = element_text(family = my_font, size = 9)) + 
  theme(axis.text.x = element_text(family = my_font)) + 
  theme(plot.caption = element_text(family = my_font, size = 8, face = "italic")) + 
  theme(plot.title = element_text(family = my_font, size = 15)) -> f1 

grid.newpage()

f1 %>% 
  left_align(c("title", "caption")) %>% 
  add_econodist_legend(
    econodist_legend_grob(
      tenth_lab = "25th Percentile", 
      ninetieth_lab = "75th Percentile", 
      med_lab = "Median", 
      med_col = "firebrick", 
      family = my_font, 
      label_size = 11,
    ), 
    below = "title"
  ) %>% 
  grid.draw()