Prepare Panel Data for Econometric Analysis from VHLSS

R Econometrics Series

Introduction

Panel Data là form dữ liệu thường được sử dụng trong nhiều phân tích và nghiên cứu. Tuy vậy, việc xử lí - chuẩn bị (Data Pre-processing) dữ liệu cho các mô hình phân tích sử dụng Panel Data có thể gây bối rối và lúng túng. Post này hướng dẫn việc khai thác, chuẩn bị Panel Data từ các 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.

Giả sử chúng ta sử dụng bộ số liệu HO2.dta từ hai năm VHLSS 2016 và VHLSS 2014 với kí hiệu Y = m4dtn, X = m4atn cho mô hình nghiên cứu đề xuất sau:

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

Trước hết đọc HO2.dta:

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

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

# Ho2.dta, 2014-2016: 

read_dta("F:/VHLSS2016/HO2.dta") -> ho2_2016
read_dta("F:/VHLSS2014_Households/Ho2.dta") -> ho2_2014

Rồi viết hàm add_zero() để thêm vào các zeros như đã phân tích chi tiết ở post trước:

#========================
#  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, na.rm = TRUE), 
           delta = n_max - n_digits, 
           pre = strrep("0", times = delta), 
           full_code = str_c(pre, x_text)) %>% 
    pull(full_code) %>% 
    return()
}

Rồi sử dụng hàm đã viết để tạo ra mã định danh cho hộ gia đình:

# Create household codes: 

ho2_2016 %>% 
  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)) %>% 
  mutate(year = 2016) -> ho2_2016


ho2_2014 %>% 
  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)) %>% 
  mutate(year = 2014) -> ho2_2014

Prepare Panel Data

Ở cấp độ hộ gia đình thì GSO chọn mẫu khảo sát kiểu như sau: năm 2014 nếu chọn 1000 hộ gia đình khảo sát thì kì khảo sát kết tiếp, tức là vào năm 2016 sẽ chỉ giữ lại 50% số hộ đã khảo sát ở năm 2014 để tiếp tục khảo sát. Số thiếu còn lại là các hộ gia đình mới. Còn vào kì khảo sát 2018 thì số hộ gia đình đã được khảo sát ở các năm 2014, 2016 chỉ còn lại 50%. Số mới lại được bổ sung. Quy luật này được áp dụng như vậy cho tất cả các năm.

Tất nhiên con số 50% chỉ là tương đối và có thể sai lệch chút ít.

Để biết một hộ cụ thể ở năm 2016 có được khảo sát ở năm 2014 hay chưa thì phải căn cứ vào cột biến tinh_14 của bộ dữ liệu HO1.dta của năm 2016. Trước hết đọc bộ dữ liệu này và tạo ra cột biến h_code với hàm ý đây là mã hộ gia đình vào năm 2014 của hộ:

# HO1.dta, 2016: 
read_dta("F:/VHLSS2016/HO1.dta") -> ho1_2016


ho1_2016 %>% 
  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)) %>% 
  mutate(h_code14 = str_c(add_zero(tinh14), 
                          add_zero(huyen14), 
                          add_zero(xa14), 
                          add_zero(diaban14), 
                          add_zero(hoso14))) -> ho1_2016

Chúng ta xem một số quan sát ở Table 1:

ho1_2016 %>% 
  select(tinh, tinh14, h_code, h_code14) %>% 
  head() %>% 
  kbl(caption = "Table 1: Some Observations from HO1.dta, 2016", escape = TRUE) %>%
  kable_classic(full_width = FALSE, html_font = "Cambria")

Table 1: Some Observations from HO1.dta, 2016

tinh	tinh14	h_code	h_code14
1	1	010010000400814	010010000400814
1	1	010010000400815	010010000400815
1	0	010010000400819	NA
1	1	010010000700613	010010000700613
1	1	010010000700614	010010000700614
1	0	010010000700620	NA

Ở đây cột tinh14 có giá trị bằng 0 có nghĩa là quan sát này không được khảo sát ở năm 2014. Như vậy dòng 1 chỉ ra rằng hộ có h_code = 010010000400814 được khảo sát trong cả hai năm liên tiếp. Trong khi đó hộ có n_code = 010010000400819 thì không được khảo sát ở năm 2014. Như vậy muốn có Panel Data cho hai năm 2014 và 2016 thì ta chỉ chọn những quan sát mà xuất hiện ở cả hai năm, tức là phải thỏa mãn h_code == h_code14. Chúng ta xem xét một số quan sát vi phạm điều kiện này ở Table 2:

ho1_2016 %>% 
  filter(h_code != h_code14) %>% 
  select(tinh, tinh14, h_code, h_code14) %>% 
  head() %>% 
  kbl(caption = "Table 2: Some Observations from HO1.dta, 2016", escape = TRUE) %>%
  kable_classic(full_width = FALSE, html_font = "Cambria")

Table 2: Some Observations from HO1.dta, 2016

tinh	tinh14	h_code	h_code14
1	1	010010003401013	010010003401019
1	1	010010003401015	010010003401020
1	1	010020007901013	010020007901020
1	1	010030011201714	010030011201720
1	1	010060019006014	010060019006019
1	1	010060019901615	010060019901619

Theo logic này chúng ta chọn ra những hộ gia đình mà, theo lí thuyết, sẽ được khảo sát đồng thời cả ở hai năm 2014 và 2016:

# Common household codes: 
ho1_2016 %>% 
  filter(h_code == h_code14) %>% 
  pull(h_code) -> common_hcode

Đến đây có thể tạo ra Panel Data như sau:

# Select some columns: 
some_columns <- c("h_code", "m4atn", "m4dtn", "year")

ho2_2016 %>% 
  filter(h_code %in% common_hcode) %>% 
  select(some_columns) -> ho2_2016_mini

ho2_2014 %>% 
  filter(h_code %in% common_hcode) %>% 
  select(some_columns) -> ho2_2014_mini

#--------------------------
#  Unbalanced Panel Data
#--------------------------

# Join the two datasets: 

bind_rows(ho2_2016_mini, ho2_2014_mini) -> unbalanced_panel_data

Mặc dù chúng ta kì vọng rằng những hộ gia đình nào đã được khảo sát năm 2016 thì cũng sẽ có thông tin ở năm 2014 nhưng sự thực thì không phải vậy: năm 2016 có 3878 hộ trong khi năm 2014 có 3874 hộ. Trước hết chúng ta có thể chỉ ra 4 hộ không có mặt ở năm 2014 như sau:

# Observations not belong to 2014: 

unbalanced_panel_data %>% 
  group_by(h_code) %>% 
  count() %>% 
  filter(n == 1) %>% 
  pull(h_code) -> some_cases

some_cases

## [1] "012680954700416" "012781017401113" "383991600000913" "384041627900112"

Như thế có nghĩa là, Panel Data mà chúng ta thu được ở trên là loại không cân bằng (Unbalanced). Đến đây chúng ta có thể “chỉ thị” cho R làm việc với Panel Data bằng sử dụng hàm pdata.frame() của thư viện plm:

# Prepare panel data: 

library(plm)

pdata.frame(unbalanced_panel_data, 
            index = c("h_code","year"), 
            row.names = FALSE, 
            drop.index = FALSE) -> plm_data_unbalanced

# Some observations: 
head(plm_data_unbalanced) %>% 
  kbl(caption = "Table 2: Some Observations from Unbalanced Panel Data", escape = TRUE) %>%
  kable_classic(full_width = TRUE, html_font = "Cambria")

Table 2: Some Observations from Unbalanced Panel Data

h_code	m4atn	m4dtn	year
010010000400814	47280	20050	2014
010010000400814	42000	60000	2016
010010000400815	297600	295000	2014
010010000400815	415600	155000	2016
010010000700613	247500	5000	2014
010010000700613	232000	32000	2016

Panel Data Analysis

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(m4dtn ~ m4atn, 
                    data = plm_data_unbalanced, 
                    model = "pooling") 

# Fixed Effects Model: 
panel_fe <- plm(m4dtn ~ m4atn, 
                data = plm_data_unbalanced, 
                model = "within")  

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

Rồi so sánh kết quả của ba cách tiếp cận này:

# Compare results: 

library(stargazer)

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

## 
## Table 3: Compare Panel Data Regression Results
## ==========================================================================
##                                   Dependent variable:                     
##              -------------------------------------------------------------
##                                          m4dtn                            
##                       Pooled              Fixed Effects     Random Effects
##                        (1)                     (2)               (3)      
## --------------------------------------------------------------------------
## m4atn                0.042***                -0.025*           0.031***   
##                      (0.006)                 (0.013)           (0.006)    
##                                                                           
## Constant          10,430.220***                             10,997.740*** 
##                     (507.339)                                 (604.861)   
##                                                                           
## --------------------------------------------------------------------------
## Observations          7,752                   7,752             7,752     
## R2                    0.007                   0.001             0.003     
## Adjusted R2           0.007                  -0.999             0.003     
## F Statistic  56.086*** (df = 1; 7750) 3.647* (df = 1; 3873)   24.561***   
## ==========================================================================
## Note:                                          *p<0.1; **p<0.05; ***p<0.01

Nếu muốn chúng ta có thể chạy lại ba mô hình trên với Balanced Panel Data dưới đây:

plm_data_unbalanced %>% 
  filter(!h_code %in% some_cases) -> plm_data_balanced

Final Notes

Post này không trình bày mọi khía cạnh của khai thác - chuẩn bị dữ liệu từ các VHLSS. Để đơn giản, post này chỉ sử dụng dữ liệu của hai năm liên tiếp. Tuy nhiên những logic trình bày trong post này có thể áp dụng để tạo ra Panel Data cho khoảng thời gian dài hơn.

Việc chuẩn bị Panel Data được trình bày trong Post này cũng có thể áp dụng cho các bộ dữ liệu khác như bộ dữ liệu khảo sát doanh nghiệp vừa và nhỏ SME - Small and Medium Enterprise Survey. Với trường hợp của SME thì việc này là đơn giản hơn vì chúng ta chỉ cần sử dụng “Mã số thuế” của doanh nghiệp thay vì tạo cột biến h_code.