Production Cost and Household Income: Some Insights from VHLSS 2018

R Daily Graph Series

Introduction

Bộ số liệu VHLSS (Household Living Standards Survey) được GSO - cơ quan thống kê quốc gia của Việt Nam tiến hành khảo sát cứ mỗi hai năm một lần. Gần đây nhất là năm 2018 gồm 44 files dữ liệu (đuôi là .dta của phần mềm Stata, tham khảo thêm thông tin về VHLSS 2018 tại đây). Mỗi một file dữ liệu là thông tin khảo sát ở cấp độ hộ gia đình hoặc cá nhân về các khía cạnh khác nhau của đời sống như Giáo Dục, Y Tế, Thu Nhập, Đầu Tư - Chi Tiêu. Post này sẽ hướng dẫn chi tiết cách thức xử lí bộ số liệu về Đầu Tư, Chi Tiêu và Thu Nhập của hộ gia đình được lưu trong file có tên HO3.dta gồm 366712 quan sát, 10 biến số. Tại thời điểm viết bài này, do không thể public bộ dữ liệu, bạn đọc quan tâm đến bộ số liệu này có thể email đến địa chỉ nguyenchidung0209@gmail.com.

R codes xử lí bộ số liệu này có thể được áp dụng để xử lí các bộ dữ liệu khác hoặc các bộ VHLSS của những năm trước đó.

Data Pre-processing

Trước hết load một số R package và đọc bộ số liệu HO3.dta:

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

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

# Define data paths: 

dir("F:/VHLSS2018_Household", full.names = TRUE) -> allPaths

# Import HO3.dta: 

read_dta(allPaths[3]) -> ho3

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ố. 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 for HO3.dta 
#==================================

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

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

Khai thác thông tin về mã tỉnh và tên tỉnh:

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

Một số mã tỉnh và tên tỉnh tướng ứng ở Table 2:

df_province %>% 
  head() %>% 
  select(province_eng, province_code) %>% 
  kbl(caption = "Table 2: Province Code", escape = TRUE) %>%
  kable_classic(full_width = TRUE, html_font = "Cambria")

Table 2: Province Code

province_eng	province_code
Ha Noi	01
Ha Giang	02
Cao Bang	04
Bac Kan	06
Tuyen Quang	08
Lao Cai	10

Joint hai bộ số liệu và xem một số quan sát (Table 3):

# Join the two data sets: 

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

ho3 %>% 
  filter(!duplicated(h_code)) %>% 
  group_by(province_eng) %>% 
  count(sort = TRUE) %>% 
  ungroup() %>% 
  mutate(total_obs = sum(n)) %>% 
  kbl(caption = "Table 4: Number of Observation by Province", escape = TRUE) %>%
  kable_classic(full_width = TRUE, html_font = "Cambria")

Table 4: Number of Observation by Province

province_eng	n	total_obs
Ho Chi Minh	1755	45839
Ha Noi	1575	45839
Thanh Hoa	1230	45839
Nghe An	1125	45839
Dong Nai	1035	45839
Nam Dinh	975	45839
Thai Binh	945	45839
Hai Phong	930	45839
Hai Duong	915	45839
Binh Duong	885	45839
Tien Giang	855	45839
Bac Giang	840	45839
Dong Thap	840	45839
Dak Lak	825	45839
Binh Dinh	810	45839
Kien Giang	810	45839
Quang Nam	795	45839
Long An	780	45839
Phu Tho	780	45839
Ben Tre	765	45839
Hung Yen	735	45839
Quang Ngai	735	45839
Quang Ninh	735	45839
Thai Nguyen	734	45839
Soc Trang	720	45839
An Giang	705	45839
Lam Dong	705	45839
Bac Ninh	690	45839
Ca Mau	690	45839
Can Tho	690	45839
Gia Lai	690	45839
Khanh Hoa	690	45839
Binh Thuan	675	45839
Tay Ninh	675	45839
Thua Thien Hue	675	45839
Vinh Long	675	45839
Vinh Phuc	675	45839
Ba Ria-Vung Tau	660	45839
Ninh Binh	645	45839
Tra Vinh	645	45839
Son La	630	45839
Da Nang	615	45839
Ha Nam	615	45839
Phu Yen	615	45839
Binh Phuoc	600	45839
Quang Binh	600	45839
Bac Lieu	570	45839
Hoa Binh	570	45839
Tuyen Quang	570	45839
Yen Bai	570	45839
Hau Giang	555	45839
Lang Son	540	45839
Ha Giang	525	45839
Bac Kan	510	45839
Dak Nong	510	45839
Dien Bien	510	45839
Kon Tum	510	45839
Lai Chau	510	45839
Lao Cai	510	45839
Ninh Thuan	510	45839
Quang Tri	510	45839
Cao Bang	495	45839
Ha Tinh	375	45839

Dễ thấy rằng có 45839 hộ được khảo sát (con số này bé hơn một chút so với 46995 trong báo cáo trên của GSO, có lẽ bộ data đang sử dụng là chưa đầy đủ). TP.HCM, Hà Nội và Thanh Hóa có số hộ được khảo sát là cao nhất. Đến đây chúng ta có thể tính toán các con số tổng, median, 25th và 75 percentiles (cho SXKD, chi khác, thu nhập) ở cấp độ hộ gia đình:

# Calculate some metrics at household level: 

ho3 %>% 
  group_by(h_code, province_eng) %>% 
  summarise(total_chi_SXKD = sum(chisxkd, na.rm = TRUE), 
            total_chi_khac = sum(chikhac, na.rm = TRUE), 
            total_thunhap = sum(thunhap, na.rm = TRUE)) %>% 
  ungroup() -> df_thu_chi


df_thu_chi %>% 
  group_by(province_eng) %>% 
  summarise(avg_income = mean(total_thunhap), 
            th25 = quantile(total_thunhap, 0.25), 
            th50 = quantile(total_thunhap, 0.50), 
            th75 = quantile(total_thunhap, 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

Gaps in Household Income

Thu nhập hộ gia đình (Figure 1) có mấy điểm chính sau:

• Với Đà Nẵng thì 25% số hộ có thu nhập hàng năm thấp hơn 250 triệu, 25% số hộ có thu nhập cao hơn 700 triệu. Đây cũng là tình có thu nhập trung vị (median) hộ gia đình hàng năm cao nhất cả nước.
• Điện Biên là tỉnh nghèo nhất nước: 50% số hộ ở tỉnh này có thu nhập bé hơn 80 triệu.

R codes cho 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, 800), breaks = seq(0, 800, 100), 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 = 8.5)) + 
  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()

Production Cost - Household Income Relationship

Để đánh giá ảnh hưởng của chi phí sản xuất (Production Cost) lên thu nhập hộ gia đình có tính đến sự khác biệt giữa các tỉnh, xét mô hình đơn giản sau:

\[Income_i = \beta_0 + \beta_1 Province_{i} + \beta_2 ProductionCost_{i} + u_i \ , \ i=1,\dots,n.\] Trong mô hình này, Province được coi là dummy variable mà hệ số hồi quy của chúng sẽ cho ta biết sự khác biệt về ảnh hưởng của chi cho sản xuất kinh doanh lên thu nhập hộ gia đình. Kết quả của mô hình này được trình bày ở Table 4 dưới đây:

# Regression: 
df_thu_chi %>% 
  mutate(Province = province_eng, 
         ProductionCost = total_chi_SXKD / 1000, 
         Income = total_thunhap / 1000) %>% # Rename for convinience. 
  lm(Income ~ Province + ProductionCost, data = .) -> reg

# Extract coefficients: 
broom::tidy(reg) %>% 
  mutate(Variable = str_replace_all(term, "Province", "")) -> df_coeff

# Show main regression results (Table 3): 
df_coeff %>% 
  mutate(Estimate = round(estimate, 1), Std.error = round(std.error, 3), P.value = round(p.value, 3)) %>% 
  select(Variable, Estimate, Std.error, P.value) -> df_coeff

df_coeff %>% 
  filter(Variable %in% c("(Intercept)", "ProductionCost")) %>% 
  bind_rows(df_coeff %>% 
              filter(!Variable %in% c("(Intercept)", "ProductionCost")) %>% 
              arrange(-Estimate)) %>% 
  kbl(caption = "Table 3: Income-Production Cost Relationship, R-squared = 86.92%", escape = TRUE) %>%
  kable_classic(full_width = TRUE, html_font = "Cambria")

Table 3: Income-Production Cost Relationship, R-squared = 86.92%

Variable	Estimate	Std.error	P.value
(Intercept)	213.8	9.678	0.000
ProductionCost	1.1	0.002	0.000
Da Nang	236.6	14.173	0.000
Ho Chi Minh	232.6	11.454	0.000
Ha Noi	197.4	11.640	0.000
Bac Ninh	183.8	13.765	0.000
Binh Duong	172.1	12.967	0.000
Dong Nai	140.0	12.544	0.000
Hai Phong	127.0	12.828	0.000
Quang Ninh	96.0	13.542	0.000
Vinh Phuc	71.8	13.833	0.000
Can Tho	65.8	13.756	0.000
Ba Ria-Vung Tau	54.2	13.913	0.000
Tay Ninh	40.7	13.833	0.003
Khanh Hoa	34.6	13.757	0.012
Hung Yen	26.5	13.542	0.050
Bac Giang	22.9	13.121	0.081
Tien Giang	20.6	13.068	0.116
Thua Thien Hue	20.4	13.833	0.141
Long An	19.9	13.349	0.137
Kien Giang	19.3	13.231	0.145
Hai Duong	17.7	12.873	0.169
Thai Nguyen	11.3	13.546	0.406
Quang Nam	5.1	13.289	0.699
Binh Thuan	0.3	13.833	0.983
Thai Binh	-2.1	12.784	0.870
Nam Dinh	-6.6	12.699	0.603
Ha Nam	-9.8	14.174	0.491
Ninh Binh	-11.1	13.996	0.428
Hau Giang	-12.0	14.577	0.410
Phu Tho	-15.7	13.349	0.240
Ca Mau	-15.9	13.757	0.247
Thanh Hoa	-17.3	12.135	0.155
Dong Thap	-19.0	13.120	0.147
Binh Dinh	-20.9	13.231	0.114
Ninh Thuan	-26.5	14.933	0.076
Quang Binh	-28.7	14.268	0.044
Bac Lieu	-28.9	14.469	0.046
Lao Cai	-30.9	14.932	0.038
Vinh Long	-32.6	13.833	0.018
Binh Phuoc	-33.4	14.267	0.019
Quang Tri	-33.7	14.932	0.024
Quang Ngai	-34.3	13.542	0.011
Lam Dong	-35.4	13.681	0.010
Ben Tre	-36.5	13.411	0.006
Soc Trang	-38.0	13.610	0.005
Nghe An	-44.1	12.340	0.000
Gia Lai	-45.0	13.756	0.001
Phu Yen	-46.2	14.174	0.001
Yen Bai	-50.6	14.469	0.000
Hoa Binh	-58.1	14.470	0.000
Ha Giang	-58.7	14.808	0.000
Tuyen Quang	-59.1	14.470	0.000
Ha Tinh	-63.1	16.418	0.000
Dak Lak	-64.8	13.175	0.000
Tra Vinh	-66.6	13.997	0.000
Lang Son	-66.8	14.690	0.000
Kon Tum	-75.0	14.933	0.000
Dien Bien	-87.3	14.933	0.000
Cao Bang	-87.9	15.064	0.000
Bac Kan	-95.3	14.933	0.000
Dak Nong	-98.3	14.932	0.000
Lai Chau	-102.3	14.933	0.000
Son La	-114.8	14.083	0.000

Một số kết quả chính từ Table 4:

• Chi 1 triệu cho SXKD thì trung bình thu nhập hộ gia đình tăng 1.1 triệu (ceteris paribus).
• Đà Nẵng, Tp Hồ Chí Minh, Hà Nội, Bắc Ninh, Bình Dương là những tỉnh tốt cho hoạt động đầu tư sản xuất kinh doanh. Đứng cuối bảng là Cao Bằng, Bắc Kạn, Dak Nông, Lai Châu và Sơn La.

Final Notes

Một số kết quả từ bộ dữ liệu này có thể tạo nên sự tranh cãi và bối rối. Bất cứ phản hồi nào đều hứu ích và được chào đón. Bạn đọc có thể liên hệ với địa chỉ email ở trên để lấy dữ liệu về kiểm tra - so sánh.