# ==============================================================================
# KỊCH BẢN R CHẠY TỪNG BƯỚC (INTERACTIVE STEP-BY-STEP WORKFLOW)
# Đề tài: Tác động của chất lượng thể chế và văn hoá đến sự phát triển TMĐT
# Dữ liệu: Panel_R_Complete_All_Variables_2012_2024.xlsx (Sheet: Panel_R_complete)
# Phương pháp chính: Correlated Random Effects (CRE / Mundlak), Year Fixed Effects,
# Sai số chuẩn cụm theo quốc gia (HC1) & CR2 Satterthwaite.
#
# HƯỚNG DẪN SỬ DỤNG TRÊN RSTUDIO:
# - Kịch bản này đã được chia thành 11 BƯỚC độc lập, hiển thị toàn bộ bảng biểu
# và đồ thị trực tiếp trên màn hình RStudio (Console & Plots), KHÔNG xuất ra
# thư mục R_results bên ngoài.
# - Cách chạy: Bôi đen từng khối lệnh (từ dòng BẮT ĐẦU BƯỚC đến HẾT BƯỚC)
# rồi nhấn tổ hợp phím Ctrl + Enter (trên Windows) hoặc Cmd + Enter (trên Mac)
# để chạy và quan sát ngay kết quả ở cửa sổ Console / Plots bên dưới.
# - Mẹo trong RStudio: Bạn có thể nhấn nút "Outline" (góc trên bên phải cửa sổ
# code) để nhảy nhanh giữa các Bước 0 -> Bước 11.
# ==============================================================================
# ==============================================================================
# BƯỚC 0: CÀI ĐẶT GÓI, NẠP THƯ VIỆN & THIẾT LẬP MÔI TRƯỜNG ====
# -> Cách chạy: Bôi đen toàn bộ Bước 0 và nhấn Ctrl + Enter.
# ==============================================================================
options(stringsAsFactors = FALSE, scipen = 5, width = 120, digits = 4)
set.seed(20261110)
# Tuỳ chỉnh loại trừ quốc gia thủ công (nếu cần) và cờ tự động loại dữ liệu lỗi:
EXCLUDE_ISO3 <- character() # Ví dụ: c("MEX") nếu muốn tự loại Mexico
AUTO_DROP_INVALID <- TRUE # TRUE: Tự động loại quốc gia vi phạm [FIX 1]
required_packages <- c(
"readxl", "dplyr", "janitor", "plm", "lmtest", "sandwich", "car",
"clubSandwich", "ggplot2", "corrplot", "readr"
)
new_packages <- setdiff(required_packages, rownames(utils::installed.packages()))
if (length(new_packages) > 0L) {
message("Đang cài đặt các gói còn thiếu: ", paste(new_packages, collapse = ", "))
install.packages(new_packages, repos = "https://cloud.r-project.org")
}
invisible(lapply(required_packages, function(pkg) {
if (!requireNamespace(pkg, quietly = TRUE)) {
stop("Không thể nạp gói bắt buộc: ", pkg)
}
}))
## Registered S3 method overwritten by 'clubSandwich':
## method from
## bread.mlm sandwich
# Hàm hỗ trợ chạy an toàn và ghi thông báo trực tiếp lên màn hình Console
log_note <- function(...) {
s <- paste0(..., collapse = "")
cat(sprintf("\n[THÔNG BÁO] %s\n", s))
invisible(s)
}
safe_call <- function(expr, label) {
tryCatch(expr, error = function(e) {
cat(sprintf("\n[LỖI KHI CHẠY] %s: %s\n", label, conditionMessage(e)))
NULL
})
}
cat("\n>>> HOÀN TẤT BƯỚC 0: Tất cả thư viện đã sẵn sàng!\n")
##
## >>> HOÀN TẤT BƯỚC 0: Tất cả thư viện đã sẵn sàng!
# ==============================================================================
# BƯỚC 1: ĐỌC DỮ LIỆU EXCEL & KIỂM TRA CHẤT LƯỢNG DỮ LIỆU BẢNG (AUDIT) ====
# -> Lưu ý: Đặt file Excel cùng thư mục làm việc (Session -> Set Working Directory)
# hoặc nếu R không tìm thấy trong thư mục hiện tại, cửa sổ chọn file sẽ hiện ra.
# -> Cách chạy: Bôi đen toàn bộ Bước 1 và nhấn Ctrl + Enter.
# ==============================================================================
base_dir <- getwd()
input_files <- list.files(
base_dir,
pattern = "^Panel_R_Complete_All_Variables_2012_2024.*[.]xlsx$",
full.names = TRUE
)
if (!length(input_files)) {
cat("\nKhông tìm thấy file Excel trong thư mục hiện tại (", base_dir, ").\n", sep = "")
if (interactive()) {
cat("Vui lòng chọn file Excel 'Panel_R_Complete_All_Variables_2012_2024.xlsx' trong hộp thoại...\n")
input_path <- file.choose()
} else {
stop("Không tìm thấy file Excel trong thư mục: ", base_dir)
}
} else {
if (length(input_files) > 1L) {
log_note("Tìm thấy nhiều file Excel khớp tên. Đang sử dụng: ", basename(input_files[1]))
}
input_path <- input_files[1]
}
log_note("Đang đọc dữ liệu từ file: ", basename(input_path))
##
## [THÔNG BÁO] Đang đọc dữ liệu từ file: Panel_R_Complete_All_Variables_2012_2024.xlsx.xlsx
df <- readxl::read_excel(input_path, sheet = "Panel_R_complete") |>
janitor::clean_names() |>
dplyr::mutate(
iso3 = toupper(trimws(as.character(iso3))),
country = as.character(country),
year = as.integer(year)
) |>
dplyr::arrange(iso3, year) |>
as.data.frame()
# Khai báo danh sách các nhóm biến trong mô hình
culture_vars <- c("pdi", "idv", "mas", "uai", "lto", "ivr")
wgi_vars <- c("va", "pv", "ge", "rq", "rl", "cc")
controls_full <- c(
"ln_gdppc", "ln_population", "internet_users_pct",
"private_credit_pct_gdp", "trade_open_pct_gdp", "secondary_enrollment_pct"
)
controls_stage1 <- c("ln_gdppc", "ln_population")
controls_stage2 <- c(controls_stage1, "internet_users_pct")
controls_stage3 <- c(controls_stage2, "private_credit_pct_gdp", "trade_open_pct_gdp")
controls_stage4 <- controls_full
model_vars <- unique(c("ln_ecom", "iq_avg", wgi_vars, culture_vars, controls_full))
other_needed <- c(
"iso3", "country", "year", "ecom_musd", "population_total",
"observed_only_sample", "conservative_sample", "high_risk_cell_count"
)
# Kiểm tra cấu trúc cột và khóa chính (iso3 - year)
missing_cols <- setdiff(c(model_vars, other_needed), names(df))
if (length(missing_cols)) stop("Thiếu các cột trong Excel: ", paste(missing_cols, collapse = ", "))
if (anyNA(df[c("iso3", "year")])) stop("Cột khóa iso3 hoặc year có giá trị NA.")
if (anyDuplicated(df[c("iso3", "year")])) stop("Có dòng trùng lặp cặp quốc gia - năm (iso3, year).")
if (any(df$ecom_musd <= 0, na.rm = TRUE)) stop("ecom_musd có giá trị <= 0, không thể lấy logarit.")
if (any(df$population_total <= 0, na.rm = TRUE)) stop("population_total có giá trị <= 0.")
if (any(!vapply(df[model_vars], is.numeric, logical(1)))) {
stop("Có biến mô hình không phải dạng số (numeric). Hãy kiểm tra định dạng ô trong Excel.")
}
# [FIX 1] Kiểm tra tính hợp lý của các biến kiểm soát dạng tỷ lệ phần trăm (%)
pct_vars <- c("internet_users_pct", "private_credit_pct_gdp",
"trade_open_pct_gdp", "secondary_enrollment_pct")
invalid_tbl <- do.call(rbind, lapply(pct_vars, function(v) {
non_positive <- df[[v]] <= 0
n_unique <- stats::ave(df[[v]], df$iso3, FUN = function(z) length(unique(z)))
n_years <- stats::ave(df[[v]], df$iso3, FUN = length)
flat <- n_unique == 1 & n_years > 3
hit <- non_positive | flat
if (!any(hit)) return(NULL)
data.frame(
iso3 = df$iso3[hit], year = df$year[hit], variable = v,
value = df[[v]][hit],
problem = ifelse(non_positive[hit], "Tỷ lệ phần trăm <= 0", "Không đổi suốt các năm")
)
}))
if (!is.null(invalid_tbl)) {
bad_iso <- sort(unique(invalid_tbl$iso3))
cat("\n=== [CẢNH BÁO DỮ LIỆU - FIX 1] PHÁT HIỆN GIÁ TRỊ KIỂM SOÁT BẤT THƯỜNG ===\n")
print(invalid_tbl, row.names = FALSE)
if (AUTO_DROP_INVALID) {
EXCLUDE_ISO3 <- union(EXCLUDE_ISO3, bad_iso)
log_note("Đã tự động loại quốc gia bị lỗi khỏi phân tích: ", paste(bad_iso, collapse = ", "))
}
} else {
cat("\n[KIỂM TRA FIX 1]: Tất cả các biến kiểm soát tỷ lệ % đều hợp lệ (PASS).\n")
}
##
## [KIỂM TRA FIX 1]: Tất cả các biến kiểm soát tỷ lệ % đều hợp lệ (PASS).
# [FIX 7] Kiểm tra đơn vị của ecom_musd (triệu USD hiện hành)
if ("gdppc_const2015_usd" %in% names(df)) {
ecom_share <- df$ecom_musd * 1e6 / (df$gdppc_const2015_usd * df$population_total)
ecom_step <- stats::ave(log(df$ecom_musd), df$iso3, FUN = function(z) c(NA, diff(z)))
unit_bad <- ecom_share < 0.001 | (is.finite(ecom_step) & abs(ecom_step) > 2.3)
if (any(unit_bad)) {
cat("\n=== [LỖI ĐƠN VỊ - FIX 7] CÁC DÒNG CÓ ECOM_MUSD LỆCH ĐƠN VỊ ===\n")
print(df[unit_bad, c("iso3", "year", "ecom_musd")], row.names = FALSE)
stop("ecom_musd có giá trị nhỏ/lớn bất thường (lệch khoảng 1.000 lần) tại: ",
paste(sort(unique(df$iso3[unit_bad])), collapse = ", "))
} else {
cat("[KIỂM TRA FIX 7]: Đơn vị biến phụ thuộc ecom_musd đồng nhất (PASS).\n")
}
}
## [KIỂM TRA FIX 7]: Đơn vị biến phụ thuộc ecom_musd đồng nhất (PASS).
if (length(EXCLUDE_ISO3)) {
log_note("Các quốc gia bị loại khỏi toàn bộ phân tích: ", paste(EXCLUDE_ISO3, collapse = ", "))
df <- df[!(df$iso3 %in% EXCLUDE_ISO3), , drop = FALSE]
}
# Tổng hợp thông tin cấu trúc bảng (Panel Audit) và in ra màn hình
pdata <- plm::pdata.frame(df, index = c("iso3", "year"), drop.index = FALSE)
panel_n <- nrow(df)
panel_countries <- dplyr::n_distinct(df$iso3)
panel_years <- dplyr::n_distinct(df$year)
balanced <- plm::is.pbalanced(pdata)
consecutive <- all(plm::is.pconsecutive(pdata))
consistency <- data.frame(
Chi_tieu_kiem_tra = c(
"Tổng số quan sát (N)",
"Số quốc gia (N_countries)",
"Số năm (T)",
"Bảng cân bằng (Balanced Panel: 1=Có, 0=Không)",
"Chuỗi năm liên tục (Consecutive Years: 1=Có, 0=Không)",
"Số dòng trùng khóa (iso3 - year)",
"Sai số tuyệt đối lớn nhất của công thức IQ_avg",
"Sai số tuyệt đối lớn nhất của công thức ln_ecom"
),
Gia_tri = c(
panel_n, panel_countries, panel_years, as.numeric(balanced),
as.numeric(consecutive), 0,
max(abs(df$iq_avg - rowMeans(df[wgi_vars]))),
max(abs(df$ln_ecom - log(df$ecom_musd)))
)
)
cat("\n================ BẢNG 1A: KẾT QUẢ KIỂM TRA CẤU TRÚC PANEL ================\n")
##
## ================ BẢNG 1A: KẾT QUẢ KIỂM TRA CẤU TRÚC PANEL ================
print(consistency, row.names = FALSE)
## Chi_tieu_kiem_tra Gia_tri
## Tổng số quan sát (N) 4.810e+02
## Số quốc gia (N_countries) 3.700e+01
## Số năm (T) 1.300e+01
## Bảng cân bằng (Balanced Panel: 1=Có, 0=Không) 1.000e+00
## Chuỗi năm liên tục (Consecutive Years: 1=Có, 0=Không) 1.000e+00
## Số dòng trùng khóa (iso3 - year) 0.000e+00
## Sai số tuyệt đối lớn nhất của công thức IQ_avg 3.553e-15
## Sai số tuyệt đối lớn nhất của công thức ln_ecom 5.151e-14
missing_tbl <- data.frame(
variable = model_vars,
missing_n = vapply(df[model_vars], function(x) sum(is.na(x)), integer(1)),
missing_pct = 100 * vapply(df[model_vars], function(x) mean(is.na(x)), numeric(1))
)
cat("\n================ BẢNG 1B: KIỂM TRA GIÁ TRỊ KHUYẾT THIẾU (MISSING) ================\n")
##
## ================ BẢNG 1B: KIỂM TRA GIÁ TRỊ KHUYẾT THIẾU (MISSING) ================
print(missing_tbl, row.names = FALSE)
## variable missing_n missing_pct
## ln_ecom 0 0
## iq_avg 0 0
## va 0 0
## pv 0 0
## ge 0 0
## rq 0 0
## rl 0 0
## cc 0 0
## pdi 0 0
## idv 0 0
## mas 0 0
## uai 0 0
## lto 0 0
## ivr 0 0
## ln_gdppc 0 0
## ln_population 0 0
## internet_users_pct 0 0
## private_credit_pct_gdp 0 0
## trade_open_pct_gdp 0 0
## secondary_enrollment_pct 0 0
country_count <- df |>
dplyr::group_by(iso3, country) |>
dplyr::summarise(
n_years = dplyr::n_distinct(year),
first_year = min(year),
last_year = max(year),
.groups = "drop"
)
cat("\n================ BẢNG 1C: DANH SÁCH QUỐC GIA TRONG MẪU (37 NƯỚC) ================\n")
##
## ================ BẢNG 1C: DANH SÁCH QUỐC GIA TRONG MẪU (37 NƯỚC) ================
print(as.data.frame(country_count), row.names = FALSE)
## iso3 country n_years first_year last_year
## AUS Australia 13 2012 2024
## AUT Austria 13 2012 2024
## BEL Belgium 13 2012 2024
## BGR Bulgaria 13 2012 2024
## CAN Canada 13 2012 2024
## CHN China 13 2012 2024
## COL Colombia 13 2012 2024
## CZE Czechia 13 2012 2024
## DEU Germany 13 2012 2024
## DNK Denmark 13 2012 2024
## ESP Spain 13 2012 2024
## EST Estonia 13 2012 2024
## FIN Finland 13 2012 2024
## FRA France 13 2012 2024
## GBR United Kingdom 13 2012 2024
## GRC Greece 13 2012 2024
## HKG China, Hong Kong SAR 13 2012 2024
## HRV Croatia 13 2012 2024
## HUN Hungary 13 2012 2024
## JPN Japan 13 2012 2024
## KOR Republic of Korea 13 2012 2024
## LTU Lithuania 13 2012 2024
## LUX Luxembourg 13 2012 2024
## MEX Mexico 13 2012 2024
## MLT Malta 13 2012 2024
## MYS Malaysia 13 2012 2024
## NLD Netherlands (Kingdom of the) 13 2012 2024
## NOR Norway 13 2012 2024
## PHL Philippines 13 2012 2024
## POL Poland 13 2012 2024
## PRT Portugal 13 2012 2024
## SGP Singapore 13 2012 2024
## SRB Serbia 13 2012 2024
## SVK Slovakia 13 2012 2024
## SVN Slovenia 13 2012 2024
## THA Thailand 13 2012 2024
## USA United States of America 13 2012 2024
# ==============================================================================
# BƯỚC 2: THỐNG KÊ MÔ TẢ, MA TRẬN TƯƠNG QUAN, VIF & ĐỒ THỊ ====
# -> Cách chạy: Bạn có thể chạy từng mục nhỏ (2A, 2B, 2C, 2D) hoặc bôi đen
# cả Bước 2 rồi nhấn Ctrl + Enter. Đồ thị sẽ hiện ngay ở tab "Plots".
# ==============================================================================
# --- 2A. Thống kê mô tả toàn mẫu và theo từng năm ---
descriptive_vars <- unique(c("ecom_musd", "ln_ecom", "iq_avg", wgi_vars,
culture_vars, controls_full))
descriptive_stats <- do.call(rbind, lapply(descriptive_vars, function(v) {
x <- df[[v]]
data.frame(
variable = v,
N = sum(!is.na(x)),
Mean = round(mean(x, na.rm = TRUE), 3),
SD = round(sd(x, na.rm = TRUE), 3),
Min = round(min(x, na.rm = TRUE), 3),
P25 = round(unname(quantile(x, .25, na.rm = TRUE)), 3),
Median = round(median(x, na.rm = TRUE), 3),
P75 = round(unname(quantile(x, .75, na.rm = TRUE)), 3),
Max = round(max(x, na.rm = TRUE), 3)
)
}))
cat("\n================ BẢNG 2A: THỐNG KÊ MÔ TẢ CÁC BIẾN TRONG MÔ HÌNH ================\n")
##
## ================ BẢNG 2A: THỐNG KÊ MÔ TẢ CÁC BIẾN TRONG MÔ HÌNH ================
print(descriptive_stats, row.names = FALSE)
## variable N Mean SD Min P25 Median P75 Max
## ecom_musd 481 584557.343 1568384.033 869.000 20988.000 103807.000 285470.510 12407741.000
## ln_ecom 481 11.439 1.991 6.767 9.952 11.550 12.562 16.334
## iq_avg 481 0.884 0.675 -0.626 0.416 0.987 1.449 1.900
## va 481 0.817 0.766 -1.741 0.295 1.031 1.418 1.870
## pv 481 0.515 0.639 -1.444 0.195 0.713 0.969 1.418
## ge 481 1.113 0.654 -0.197 0.610 1.150 1.694 2.241
## rq 481 1.050 0.641 -0.267 0.514 1.083 1.632 2.194
## rl 481 0.908 0.800 -1.115 0.350 1.083 1.543 2.097
## cc 481 0.901 0.893 -0.961 0.158 0.825 1.675 2.326
## pdi 481 57.568 21.700 11.000 40.000 60.000 70.000 104.000
## idv 481 50.568 23.365 13.000 27.000 55.000 69.000 91.000
## mas 481 50.541 22.248 8.000 39.000 50.000 64.000 110.000
## uai 481 66.703 24.535 8.000 50.000 70.000 86.000 112.000
## lto 481 54.306 21.838 13.098 36.020 52.141 70.025 100.000
## ivr 481 47.038 19.858 15.625 29.464 45.536 62.723 97.321
## ln_gdppc 481 10.075 0.820 7.837 9.484 10.157 10.710 11.616
## ln_population 481 16.614 1.672 12.948 15.535 16.193 17.762 21.069
## internet_users_pct 481 81.448 14.236 26.460 74.700 85.001 92.000 99.769
## private_credit_pct_gdp 481 89.026 44.406 19.057 50.384 86.661 116.009 264.447
## trade_open_pct_gdp 481 122.059 84.790 23.140 66.429 100.359 150.966 442.620
## secondary_enrollment_pct 481 107.903 16.419 73.246 99.329 104.321 112.536 164.080
by_year <- df |>
dplyr::group_by(year) |>
dplyr::summarise(
N_countries = dplyr::n_distinct(iso3),
median_ecom_musd = round(median(ecom_musd, na.rm = TRUE), 2),
mean_ln_ecom = round(mean(ln_ecom, na.rm = TRUE), 3),
mean_iq_avg = round(mean(iq_avg, na.rm = TRUE), 3),
.groups = "drop"
)
cat("\n================ BẢNG 2B: XU HƯỚNG TRUNG VỊ TMĐT & THỂ CHẾ THEO NĂM ================\n")
##
## ================ BẢNG 2B: XU HƯỚNG TRUNG VỊ TMĐT & THỂ CHẾ THEO NĂM ================
print(as.data.frame(by_year), row.names = FALSE)
## year N_countries median_ecom_musd mean_ln_ecom mean_iq_avg
## 2012 37 77287 11.00 0.896
## 2013 37 88988 11.09 0.882
## 2014 37 90203 11.20 0.900
## 2015 37 89612 11.14 0.904
## 2016 37 76039 11.19 0.895
## 2017 37 92999 11.34 0.895
## 2018 37 108732 11.48 0.902
## 2019 37 127184 11.50 0.904
## 2020 37 120861 11.51 0.877
## 2021 37 110236 11.68 0.884
## 2022 37 155009 11.75 0.868
## 2023 37 167098 11.87 0.853
## 2024 37 155923 11.95 0.831
# --- 2B. Ma trận tương quan Pearson & Biểu đồ Heatmap trên cửa sổ Plots ---
cor_vars <- unique(c("iq_avg", wgi_vars, culture_vars, controls_full))
cor_matrix <- stats::cor(df[cor_vars], use = "pairwise.complete.obs")
cat("\n================ BẢNG 2C: MA TRẬN TƯƠNG QUAN PEARSON ================\n")
##
## ================ BẢNG 2C: MA TRẬN TƯƠNG QUAN PEARSON ================
print(round(cor_matrix, 3))
## iq_avg va pv ge rq rl cc pdi idv mas uai lto ivr
## iq_avg 1.000 0.841 0.825 0.950 0.940 0.990 0.965 -0.702 0.633 -0.248 -0.230 0.183 0.150
## va 0.841 1.000 0.664 0.692 0.695 0.833 0.728 -0.680 0.692 -0.236 0.138 0.009 0.254
## pv 0.825 0.664 1.000 0.723 0.692 0.796 0.716 -0.485 0.427 -0.168 -0.036 0.364 -0.175
## ge 0.950 0.692 0.723 1.000 0.926 0.943 0.956 -0.630 0.543 -0.213 -0.369 0.194 0.213
## rq 0.940 0.695 0.692 0.926 1.000 0.924 0.949 -0.635 0.573 -0.205 -0.432 0.136 0.182
## rl 0.990 0.833 0.796 0.943 0.924 1.000 0.956 -0.729 0.623 -0.275 -0.212 0.186 0.133
## cc 0.965 0.728 0.716 0.956 0.949 0.956 1.000 -0.683 0.604 -0.251 -0.366 0.155 0.184
## pdi -0.702 -0.680 -0.485 -0.630 -0.635 -0.729 -0.683 1.000 -0.649 0.265 0.091 0.017 -0.273
## idv 0.633 0.692 0.427 0.543 0.573 0.623 0.604 -0.649 1.000 0.020 -0.132 -0.045 0.298
## mas -0.248 -0.236 -0.168 -0.213 -0.205 -0.275 -0.251 0.265 0.020 1.000 0.040 0.120 0.018
## uai -0.230 0.138 -0.036 -0.369 -0.432 -0.212 -0.366 0.091 -0.132 0.040 1.000 0.025 -0.100
## lto 0.183 0.009 0.364 0.194 0.136 0.186 0.155 0.017 -0.045 0.120 0.025 1.000 -0.591
## ivr 0.150 0.254 -0.175 0.213 0.182 0.133 0.184 -0.273 0.298 0.018 -0.100 -0.591 1.000
## ln_gdppc 0.910 0.718 0.702 0.911 0.869 0.901 0.912 -0.658 0.597 -0.152 -0.235 0.162 0.278
## ln_population -0.368 -0.410 -0.559 -0.196 -0.286 -0.347 -0.259 0.180 -0.062 0.363 -0.163 -0.101 0.167
## internet_users_pct 0.715 0.561 0.623 0.693 0.658 0.708 0.703 -0.489 0.476 -0.144 -0.179 0.219 0.115
## private_credit_pct_gdp 0.333 -0.005 0.153 0.446 0.427 0.361 0.446 -0.179 -0.040 -0.139 -0.477 0.101 -0.017
## trade_open_pct_gdp 0.236 -0.020 0.335 0.235 0.332 0.206 0.253 0.091 -0.136 -0.060 -0.230 0.282 -0.239
## secondary_enrollment_pct 0.547 0.550 0.366 0.476 0.493 0.543 0.556 -0.484 0.494 -0.309 0.017 -0.080 0.242
## ln_gdppc ln_population internet_users_pct private_credit_pct_gdp trade_open_pct_gdp
## iq_avg 0.910 -0.368 0.715 0.333 0.236
## va 0.718 -0.410 0.561 -0.005 -0.020
## pv 0.702 -0.559 0.623 0.153 0.335
## ge 0.911 -0.196 0.693 0.446 0.235
## rq 0.869 -0.286 0.658 0.427 0.332
## rl 0.901 -0.347 0.708 0.361 0.206
## cc 0.912 -0.259 0.703 0.446 0.253
## pdi -0.658 0.180 -0.489 -0.179 0.091
## idv 0.597 -0.062 0.476 -0.040 -0.136
## mas -0.152 0.363 -0.144 -0.139 -0.060
## uai -0.235 -0.163 -0.179 -0.477 -0.230
## lto 0.162 -0.101 0.219 0.101 0.282
## ivr 0.278 0.167 0.115 -0.017 -0.239
## ln_gdppc 1.000 -0.237 0.765 0.410 0.259
## ln_population -0.237 1.000 -0.264 0.138 -0.612
## internet_users_pct 0.765 -0.264 1.000 0.294 0.210
## private_credit_pct_gdp 0.410 0.138 0.294 1.000 0.202
## trade_open_pct_gdp 0.259 -0.612 0.210 0.202 1.000
## secondary_enrollment_pct 0.508 -0.186 0.358 0.172 -0.040
## secondary_enrollment_pct
## iq_avg 0.547
## va 0.550
## pv 0.366
## ge 0.476
## rq 0.493
## rl 0.543
## cc 0.556
## pdi -0.484
## idv 0.494
## mas -0.309
## uai 0.017
## lto -0.080
## ivr 0.242
## ln_gdppc 0.508
## ln_population -0.186
## internet_users_pct 0.358
## private_credit_pct_gdp 0.172
## trade_open_pct_gdp -0.040
## secondary_enrollment_pct 1.000
safe_call({
corrplot::corrplot(
cor_matrix, type = "upper", method = "color", tl.col = "black",
tl.cex = 0.75, tl.srt = 45, addCoef.col = "black", number.cex = 0.50,
title = "Ma tran tuong quan giua cac bien giai thich", mar = c(0, 0, 2, 0)
)
}, "Vẽ biểu đồ nhiệt ma trận tương quan (Correlation Heatmap)")
# --- 2C. Biểu đồ xu hướng doanh số TMĐT trung vị qua các năm (2012-2024) ---
safe_call({
p_trend <- ggplot2::ggplot(by_year, ggplot2::aes(x = year, y = median_ecom_musd)) +
ggplot2::geom_line(linewidth = 1, color = "#1f77b4") +
ggplot2::geom_point(size = 2.5, color = "#d62728") +
ggplot2::scale_x_continuous(breaks = sort(unique(by_year$year))) +
ggplot2::labs(
title = "Doanh số thương mại điện tử trung vị qua các năm (2012 - 2024)",
subtitle = paste0("Mẫu gồm ", panel_countries, " quốc gia (", panel_n, " quan sát)"),
x = "Năm", y = "Doanh số TMĐT trung vị (Triệu USD)"
) +
ggplot2::theme_minimal(base_size = 12)
print(p_trend)
}, "Vẽ biểu đồ xu hướng TMĐT")
# --- 2D. Kiểm tra đa cộng tuyến (VIF) ---
vif_formula <- stats::reformulate(
c("iq_avg", culture_vars, controls_full), response = "ln_ecom"
)
vif_lm <- stats::lm(vif_formula, data = df)
vif_raw <- safe_call(car::vif(vif_lm), "Tính hệ số phóng đại phương sai VIF")
if (!is.null(vif_raw)) {
if (is.matrix(vif_raw)) {
vif_tbl <- data.frame(variable = rownames(vif_raw), round(vif_raw, 3),
row.names = NULL, check.names = FALSE)
} else {
vif_tbl <- data.frame(variable = names(vif_raw), VIF = round(as.numeric(vif_raw), 3))
}
cat("\n================ BẢNG 2D: HỆ SỐ PHÓNG ĐẠI PHƯƠNG SAI (VIF) ================\n")
cat("(Lưu ý: Không đưa đồng thời iq_avg và 6 biến WGI vào cùng 1 mô hình VIF vì iq_avg là trung bình cộng)\n")
print(vif_tbl, row.names = FALSE)
}
##
## ================ BẢNG 2D: HỆ SỐ PHÓNG ĐẠI PHƯƠNG SAI (VIF) ================
## (Lưu ý: Không đưa đồng thời iq_avg và 6 biến WGI vào cùng 1 mô hình VIF vì iq_avg là trung bình cộng)
## variable VIF
## iq_avg 10.300
## pdi 2.749
## idv 3.110
## mas 1.615
## uai 1.897
## lto 2.160
## ivr 2.493
## ln_gdppc 11.893
## ln_population 3.924
## internet_users_pct 2.622
## private_credit_pct_gdp 2.212
## trade_open_pct_gdp 3.496
## secondary_enrollment_pct 1.821
# ==============================================================================
# BƯỚC 3: KHỞI TẠO CÁC HÀM ƯỚC LƯỢNG & HÀM HIỂN THỊ BẢNG KẾT QUẢ ====
# -> Bước này nạp các hàm tính toán mô hình CRE/Mundlak, sai số chuẩn cụm HC1,
# sai số chuẩn CR2 Satterthwaite và hàm in bảng tổng hợp lên màn hình R.
# -> Cách chạy: Bôi đen toàn bộ Bước 3 và nhấn Ctrl + Enter.
# ==============================================================================
make_formula <- function(y = "ln_ecom", institution = "iq_avg", culture = culture_vars,
controls = controls_full, add_mundlak = FALSE) {
explanatory <- unique(c(institution, culture, controls))
if (add_mundlak) {
explanatory <- c(explanatory, paste0(c(institution, controls), "_m"))
}
stats::reformulate(c(explanatory, "factor(year)"), response = y)
}
prepare_cre_data <- function(dat, institution = "iq_avg", culture = culture_vars,
controls = controls_full, y = "ln_ecom") {
needed <- unique(c("iso3", "year", y, institution, culture, controls))
if (length(setdiff(needed, names(dat)))) {
stop("Thiếu biến trong dữ liệu: ", paste(setdiff(needed, names(dat)), collapse = ", "))
}
dat <- dat[stats::complete.cases(dat[needed]), , drop = FALSE]
if (nrow(dat) < 30L || dplyr::n_distinct(dat$iso3) < 6L) {
stop("Số quan sát hoặc số quốc gia quá ít để ước lượng CRE.")
}
dat |>
dplyr::arrange(iso3, year) |>
dplyr::group_by(iso3) |>
dplyr::mutate(
dplyr::across(dplyr::all_of(c(institution, controls)),
~mean(.x), .names = "{.col}_m")
) |>
dplyr::ungroup() |>
as.data.frame()
}
fit_cre <- function(dat, institution = "iq_avg", culture = culture_vars,
controls = controls_full, y = "ln_ecom") {
dat <- prepare_cre_data(dat, institution, culture, controls, y)
pdat <- plm::pdata.frame(dat, index = c("iso3", "year"), drop.index = FALSE)
f <- make_formula(y, institution, culture, controls, add_mundlak = TRUE)
# [FIX 2] Tự động chuyển sang Wallace-Hussain nếu mẫu con không cân bằng làm Swamy-Arora bị suy biến
tryCatch(
plm::plm(f, data = pdat, model = "random", effect = "individual",
random.method = "swar"),
error = function(e) {
message(" [Lưu ý FIX 2] Swamy-Arora không khả dụng trên mẫu không cân bằng; tự động chuyển sang Wallace-Hussain.")
plm::plm(f, data = pdat, model = "random", effect = "individual",
random.method = "walhus")
}
)
}
fit_noncre <- function(dat, method = c("pooling", "within", "random"),
with_culture = TRUE) {
method <- match.arg(method)
cult <- if (with_culture) culture_vars else character()
f <- make_formula(culture = cult)
needed <- unique(c("iso3", "year", "ln_ecom", "iq_avg", controls_full, cult))
dat <- dat[stats::complete.cases(dat[needed]), , drop = FALSE]
pdat <- plm::pdata.frame(dat, index = c("iso3", "year"), drop.index = FALSE)
if (method == "random") {
plm::plm(f, data = pdat, model = method, effect = "individual", random.method = "swar")
} else {
plm::plm(f, data = pdat, model = method, effect = "individual")
}
}
vcov_country <- function(m) {
plm::vcovHC(m, method = "arellano", type = "HC1", cluster = "group")
}
model_n <- function(m) if (is.null(m)) NA_integer_ else length(stats::residuals(m))
model_countries <- function(m) {
if (is.null(m)) return(NA_integer_)
length(unique(as.character(plm::index(m)[[1]])))
}
empty_tidy <- function() {
data.frame(term = character(), estimate = numeric(), std_error = numeric(),
statistic = numeric(), df = numeric(), p_value = numeric())
}
get_coefs <- function(m, inference = c("cluster", "cr2", "dk", "double")) {
inference <- match.arg(inference)
if (is.null(m)) return(empty_tidy())
if (inference == "cr2") {
x <- safe_call(as.data.frame(clubSandwich::coef_test(
m, vcov = "CR2", cluster = "individual", test = "Satterthwaite"
)), "CR2 estimation")
if (is.null(x)) return(empty_tidy())
return(data.frame(
term = as.character(x$Coef), estimate = as.numeric(x$beta),
std_error = as.numeric(x$SE), statistic = as.numeric(x$tstat),
df = as.numeric(x$df_Satt), p_value = as.numeric(x$p_Satt)
))
}
V <- safe_call(switch(inference,
cluster = vcov_country(m),
dk = plm::vcovSCC(m, type = "HC1", maxlag = 2),
double = plm::vcovDC(m, type = "HC1")
), paste0("Covariance: ", inference))
if (is.null(V)) return(empty_tidy())
test <- safe_call(lmtest::coeftest(m, vcov. = V),
paste0("Coefficient test: ", inference))
if (is.null(test)) return(empty_tidy())
data.frame(
term = rownames(test), estimate = unname(test[, 1]),
std_error = unname(test[, 2]), statistic = unname(test[, 3]),
df = NA_real_, p_value = unname(test[, 4]), row.names = NULL
)
}
collect_coefs <- function(model, model_label, inference = "cluster", keep = NULL) {
x <- get_coefs(model, inference)
if (!nrow(x)) return(data.frame())
if (!is.null(keep)) x <- x[x$term %in% keep, , drop = FALSE]
if (!nrow(x)) return(data.frame())
cbind(
data.frame(model = model_label, inference = inference,
N = model_n(model), countries = model_countries(model)),
x
)
}
signif_mark <- function(p) {
if (!is.finite(p)) return("")
if (p < .01) return("***")
if (p < .05) return("**")
if (p < .10) return("*")
""
}
# Hàm in bảng hồi quy tổng hợp trực tiếp ra màn hình Console (thay cho xuất file HTML)
print_model_table <- function(models, heading, keep = NULL, inference = "cluster") {
models <- models[!vapply(models, is.null, logical(1))]
if (!length(models)) return(invisible(NULL))
tables <- lapply(models, get_coefs, inference = inference)
all_terms <- unique(unlist(lapply(tables, function(x) x$term)))
# Ẩn các biến giả năm factor(year) cho bảng gọn đẹp dễ nhìn trên Console
all_terms <- all_terms[!grepl("^factor\\(year\\)", all_terms)]
if (!is.null(keep)) all_terms <- intersect(keep, all_terms)
if (!length(all_terms)) return(invisible(NULL))
out_df <- data.frame(Variable = c(all_terms, "Observations (N)", "Countries"),
stringsAsFactors = FALSE)
for (j in seq_along(models)) {
tb <- tables[[j]]
col_vals <- character(length(all_terms) + 2L)
for (i in seq_along(all_terms)) {
z <- tb[tb$term == all_terms[i], , drop = FALSE]
if (nrow(z) == 1L) {
col_vals[i] <- sprintf("%.4f%s (%.4f)", z$estimate,
signif_mark(z$p_value), z$std_error)
} else {
col_vals[i] <- ""
}
}
col_vals[length(all_terms) + 1L] <- as.character(model_n(models[[j]]))
col_vals[length(all_terms) + 2L] <- as.character(model_countries(models[[j]]))
out_df[[names(models)[j]]] <- col_vals
}
cat("\n==============================================================================\n")
cat(heading, "\n")
cat("Ghi chú: Hệ số ước lượng kèm sai số chuẩn (", toupper(inference), ") trong ngoặc đơn.\n", sep = "")
cat("Mức ý nghĩa: * p < 0.10; ** p < 0.05; *** p < 0.01. Đã kiểm soát hiệu ứng năm.\n")
cat("==============================================================================\n")
print(out_df, row.names = FALSE, right = FALSE)
invisible(out_df)
}
cat("\n>>> HOÀN TẤT BƯỚC 3: Đã nạp toàn bộ hàm ước lượng và hàm hiển thị bảng!\n")
##
## >>> HOÀN TẤT BƯỚC 3: Đã nạp toàn bộ hàm ước lượng và hàm hiển thị bảng!
# ==============================================================================
# BƯỚC 4: ƯỚC LƯỢNG CRE VỚI 4 TẦNG BIẾN KIỂM SOÁT (M1 -> M4) ====
# -> M1: Chỉ kiểm soát GDP bình quân & Dân số
# -> M2: Thêm tỷ lệ người dùng Internet
# -> M3: Thêm Tín dụng tư nhân & Độ mở thương mại
# -> M4: Mô hình đầy đủ (thêm Giáo dục trung học)
# -> Cách chạy: Bôi đen toàn bộ Bước 4 và nhấn Ctrl + Enter.
# ==============================================================================
control_stages <- list(
M1_GDP_Pop = controls_stage1,
M2_Internet = controls_stage2,
M3_Finance_Trade = controls_stage3,
M4_Education_Full = controls_stage4
)
nested_models <- lapply(names(control_stages), function(lbl) {
safe_call(fit_cre(df, controls = control_stages[[lbl]]),
paste0("Nested CRE: ", lbl))
})
names(nested_models) <- names(control_stages)
# In bảng so sánh 4 tầng mô hình M1 -> M4 trực tiếp lên Console
print_model_table(
nested_models,
heading = "BẢNG 4A: KẾT QUẢ ƯỚC LƯỢNG CRE QUA 4 TẦNG BIẾN KIỂM SOÁT (M1 - M4)",
keep = c("(Intercept)", "iq_avg", culture_vars, controls_full)
)
##
## ==============================================================================
## BẢNG 4A: KẾT QUẢ ƯỚC LƯỢNG CRE QUA 4 TẦNG BIẾN KIỂM SOÁT (M1 - M4)
## Ghi chú: Hệ số ước lượng kèm sai số chuẩn (CLUSTER) trong ngoặc đơn.
## Mức ý nghĩa: * p < 0.10; ** p < 0.05; *** p < 0.01. Đã kiểm soát hiệu ứng năm.
## ==============================================================================
## Variable M1_GDP_Pop M2_Internet M3_Finance_Trade M4_Education_Full
## (Intercept) -22.7402*** (4.8011) -23.9933*** (3.9925) -26.4692*** (3.7917) -27.4373*** (4.1206)
## iq_avg 0.0812 (0.2823) 0.0521 (0.2699) 0.0060 (0.2805) 0.0307 (0.2801)
## pdi 0.0040 (0.0083) 0.0010 (0.0055) -0.0028 (0.0056) -0.0036 (0.0061)
## idv -0.0005 (0.0043) -0.0042 (0.0047) -0.0094** (0.0043) -0.0119** (0.0052)
## mas 0.0025 (0.0037) 0.0033 (0.0028) 0.0011 (0.0031) 0.0024 (0.0032)
## uai 0.0044 (0.0046) 0.0054 (0.0042) 0.0049 (0.0039) 0.0036 (0.0036)
## lto 0.0122** (0.0048) 0.0055 (0.0049) 0.0010 (0.0044) 0.0006 (0.0045)
## ivr -0.0032 (0.0062) -0.0058 (0.0058) -0.0091* (0.0048) -0.0103* (0.0055)
## ln_gdppc 1.6398*** (0.4092) 1.2161*** (0.3567) 1.2034*** (0.4232) 1.1367*** (0.4195)
## ln_population 0.6650 (0.8476) 0.4836 (0.8947) 0.5104 (0.8058) 0.4700 (0.8030)
## internet_users_pct 0.0110** (0.0047) 0.0117*** (0.0044) 0.0123*** (0.0047)
## private_credit_pct_gdp -0.0024 (0.0022) -0.0022 (0.0022)
## trade_open_pct_gdp 0.0019 (0.0029) 0.0017 (0.0029)
## secondary_enrollment_pct 0.0049* (0.0029)
## Observations (N) 481 481 481 481
## Countries 37 37 37 37
# In bảng tổng hợp riêng các biến chính (iq_avg và 6 chiều văn hoá)
nested_coefs <- do.call(rbind, lapply(names(nested_models), function(lbl) {
collect_coefs(nested_models[[lbl]], lbl, "cluster", keep = c("iq_avg", culture_vars))
}))
cat("\n================ BẢNG 4B: CHI TIẾT HỆ SỐ THỂ CHẾ & VĂN HOÁ QUA M1 - M4 ================\n")
##
## ================ BẢNG 4B: CHI TIẾT HỆ SỐ THỂ CHẾ & VĂN HOÁ QUA M1 - M4 ================
print(nested_coefs, digits = 4, row.names = FALSE)
## model inference N countries term estimate std_error statistic df p_value
## M1_GDP_Pop cluster 481 37 iq_avg 0.0811878 0.282302 0.28759 NA 0.77379
## M1_GDP_Pop cluster 481 37 pdi 0.0040213 0.008292 0.48499 NA 0.62792
## M1_GDP_Pop cluster 481 37 idv -0.0004845 0.004338 -0.11168 NA 0.91113
## M1_GDP_Pop cluster 481 37 mas 0.0025383 0.003692 0.68753 NA 0.49210
## M1_GDP_Pop cluster 481 37 uai 0.0044019 0.004614 0.95407 NA 0.34055
## M1_GDP_Pop cluster 481 37 lto 0.0121959 0.004840 2.51999 NA 0.01208
## M1_GDP_Pop cluster 481 37 ivr -0.0031921 0.006238 -0.51170 NA 0.60911
## M2_Internet cluster 481 37 iq_avg 0.0521153 0.269917 0.19308 NA 0.84698
## M2_Internet cluster 481 37 pdi 0.0010430 0.005452 0.19131 NA 0.84837
## M2_Internet cluster 481 37 idv -0.0042078 0.004743 -0.88717 NA 0.37545
## M2_Internet cluster 481 37 mas 0.0033279 0.002819 1.18070 NA 0.23834
## M2_Internet cluster 481 37 uai 0.0054023 0.004178 1.29315 NA 0.19662
## M2_Internet cluster 481 37 lto 0.0055064 0.004874 1.12975 NA 0.25918
## M2_Internet cluster 481 37 ivr -0.0057780 0.005804 -0.99556 NA 0.32000
## M3_Finance_Trade cluster 481 37 iq_avg 0.0060220 0.280510 0.02147 NA 0.98288
## M3_Finance_Trade cluster 481 37 pdi -0.0027960 0.005635 -0.49617 NA 0.62001
## M3_Finance_Trade cluster 481 37 idv -0.0094082 0.004275 -2.20067 NA 0.02827
## M3_Finance_Trade cluster 481 37 mas 0.0010901 0.003054 0.35689 NA 0.72134
## M3_Finance_Trade cluster 481 37 uai 0.0048554 0.003866 1.25586 NA 0.20982
## M3_Finance_Trade cluster 481 37 lto 0.0009800 0.004364 0.22458 NA 0.82241
## M3_Finance_Trade cluster 481 37 ivr -0.0090878 0.004786 -1.89869 NA 0.05825
## M4_Education_Full cluster 481 37 iq_avg 0.0306525 0.280092 0.10944 NA 0.91290
## M4_Education_Full cluster 481 37 pdi -0.0035597 0.006098 -0.58375 NA 0.55968
## M4_Education_Full cluster 481 37 idv -0.0118636 0.005211 -2.27682 NA 0.02327
## M4_Education_Full cluster 481 37 mas 0.0023906 0.003204 0.74606 NA 0.45602
## M4_Education_Full cluster 481 37 uai 0.0035679 0.003627 0.98383 NA 0.32573
## M4_Education_Full cluster 481 37 lto 0.0006345 0.004484 0.14150 NA 0.88754
## M4_Education_Full cluster 481 37 ivr -0.0102893 0.005498 -1.87154 NA 0.06192
# ==============================================================================
# BƯỚC 5: SO SÁNH POOLED OLS, FE, RE, CRE/MUNDLAK & KIỂM ĐỊNH H1 ====
# -> So sánh 4 phương pháp ước lượng dữ liệu bảng và kiểm định giả thuyết H1
# (cả tác động trong nước within-country và tác động giữa các nước between-country).
# -> Cách chạy: Bôi đen toàn bộ Bước 5 và nhấn Ctrl + Enter.
# ==============================================================================
pool_full <- safe_call(fit_noncre(df, method = "pooling", with_culture = TRUE), "Pooled OLS")
fe_model <- safe_call(fit_noncre(df, method = "within", with_culture = FALSE), "Country FE")
re_model <- safe_call(fit_noncre(df, method = "random", with_culture = TRUE), "Random Effects")
cre_full <- nested_models[["M4_Education_Full"]]
main_models <- list(
Pooled_OLS = pool_full,
Country_FE = fe_model,
Random_Effects = re_model,
CRE_Mundlak = cre_full
)
# 5A. Bảng đối chiếu 4 mô hình chính (Pooled OLS, FE, RE, CRE)
print_model_table(
main_models,
heading = "BẢNG 5A: ĐỐI CHIẾU POOLED OLS, FIXED EFFECTS (FE), RANDOM EFFECTS (RE) & CRE/MUNDLAK",
keep = c("(Intercept)", "iq_avg", culture_vars, controls_full, paste0(c("iq_avg", controls_full), "_m"))
)
##
## ==============================================================================
## BẢNG 5A: ĐỐI CHIẾU POOLED OLS, FIXED EFFECTS (FE), RANDOM EFFECTS (RE) & CRE/MUNDLAK
## Ghi chú: Hệ số ước lượng kèm sai số chuẩn (CLUSTER) trong ngoặc đơn.
## Mức ý nghĩa: * p < 0.10; ** p < 0.05; *** p < 0.01. Đã kiểm soát hiệu ứng năm.
## ==============================================================================
## Variable Pooled_OLS Country_FE Random_Effects CRE_Mundlak
## (Intercept) -24.2747*** (3.8917) -22.2111*** (2.9968) -27.4373*** (4.1206)
## iq_avg -0.1094 (0.5481) 0.0307 (0.2758) 0.0785 (0.2557) 0.0307 (0.2801)
## pdi 0.0006 (0.0074) 0.0005 (0.0080) -0.0036 (0.0061)
## idv -0.0058 (0.0039) -0.0021 (0.0057) -0.0119** (0.0052)
## mas 0.0024 (0.0034) 0.0026 (0.0041) 0.0024 (0.0032)
## uai 0.0037 (0.0037) 0.0033 (0.0059) 0.0036 (0.0036)
## lto 0.0076* (0.0042) 0.0121* (0.0066) 0.0006 (0.0045)
## ivr -0.0061 (0.0056) -0.0025 (0.0061) -0.0103* (0.0055)
## ln_gdppc 1.3515*** (0.4753) 1.1367*** (0.4130) 1.3505*** (0.3487) 1.1367*** (0.4195)
## ln_population 1.1200*** (0.0807) 0.4700 (0.7908) 1.0617*** (0.1149) 0.4700 (0.8030)
## internet_users_pct 0.0404*** (0.0114) 0.0123*** (0.0046) 0.0120** (0.0047) 0.0123*** (0.0047)
## private_credit_pct_gdp -0.0040* (0.0023) -0.0022 (0.0022) -0.0024 (0.0020) -0.0022 (0.0022)
## trade_open_pct_gdp 0.0014 (0.0017) 0.0017 (0.0028) 0.0013 (0.0024) 0.0017 (0.0029)
## secondary_enrollment_pct 0.0045 (0.0056) 0.0049* (0.0029) 0.0047* (0.0028) 0.0049* (0.0029)
## iq_avg_m -0.6472 (0.6548)
## ln_gdppc_m -0.0954 (0.4749)
## ln_population_m 0.7555 (0.8000)
## internet_users_pct_m 0.0958*** (0.0352)
## private_credit_pct_gdp_m -0.0053 (0.0034)
## trade_open_pct_gdp_m 0.0010 (0.0034)
## secondary_enrollment_pct_m 0.0032 (0.0072)
## Observations (N) 481 481 481 481
## Countries 37 37 37 37
# 5B. Kết quả chi tiết mô hình CRE đầy đủ với sai số chuẩn cụm (HC1) và CR2 Satterthwaite
main_cluster <- get_coefs(cre_full, "cluster")
main_cr2 <- get_coefs(cre_full, "cr2")
cat("\n================ BẢNG 5B: MÔ HÌNH CHÍNH CRE - SAI SỐ CHUẨN CỤM QUỐC GIA (HC1) ================\n")
##
## ================ BẢNG 5B: MÔ HÌNH CHÍNH CRE - SAI SỐ CHUẨN CỤM QUỐC GIA (HC1) ================
print(main_cluster[!grepl("^factor\\(year\\)", main_cluster$term), ], digits = 4, row.names = FALSE)
## term estimate std_error statistic df p_value
## (Intercept) -27.4373410 4.120612 -6.6586 NA 8.112e-11
## iq_avg 0.0306525 0.280092 0.1094 NA 9.129e-01
## pdi -0.0035597 0.006098 -0.5838 NA 5.597e-01
## idv -0.0118636 0.005211 -2.2768 NA 2.327e-02
## mas 0.0023906 0.003204 0.7461 NA 4.560e-01
## uai 0.0035679 0.003627 0.9838 NA 3.257e-01
## lto 0.0006345 0.004484 0.1415 NA 8.875e-01
## ivr -0.0102893 0.005498 -1.8715 NA 6.192e-02
## ln_gdppc 1.1366774 0.419451 2.7099 NA 6.988e-03
## ln_population 0.4699575 0.803016 0.5852 NA 5.587e-01
## internet_users_pct 0.0122864 0.004652 2.6410 NA 8.556e-03
## private_credit_pct_gdp -0.0022078 0.002188 -1.0089 NA 3.136e-01
## trade_open_pct_gdp 0.0016669 0.002866 0.5815 NA 5.612e-01
## secondary_enrollment_pct 0.0048980 0.002948 1.6613 NA 9.736e-02
## iq_avg_m -0.6472420 0.654794 -0.9885 NA 3.235e-01
## ln_gdppc_m -0.0953957 0.474862 -0.2009 NA 8.409e-01
## ln_population_m 0.7554798 0.800030 0.9443 NA 3.455e-01
## internet_users_pct_m 0.0957848 0.035190 2.7219 NA 6.743e-03
## private_credit_pct_gdp_m -0.0053033 0.003422 -1.5496 NA 1.219e-01
## trade_open_pct_gdp_m 0.0010400 0.003372 0.3085 NA 7.579e-01
## secondary_enrollment_pct_m 0.0032321 0.007171 0.4507 NA 6.524e-01
cat("\n================ BẢNG 5C: MÔ HÌNH CHÍNH CRE - SAI SỐ CHUẨN CR2 SATTERTHWAITE ================\n")
##
## ================ BẢNG 5C: MÔ HÌNH CHÍNH CRE - SAI SỐ CHUẨN CR2 SATTERTHWAITE ================
cat("(Khuyến nghị: Dùng cột p_value của bảng CR2 này để kết luận cho các biến văn hoá)\n")
## (Khuyến nghị: Dùng cột p_value của bảng CR2 này để kết luận cho các biến văn hoá)
print(main_cr2[!grepl("^factor\\(year\\)", main_cr2$term), ], digits = 4, row.names = FALSE)
## term estimate std_error statistic df p_value
## (Intercept) -27.4373410 5.639475 -4.86523 8.863 0.0009297
## iq_avg 0.0306525 0.310009 0.09888 16.427 0.9224331
## pdi -0.0035597 0.008050 -0.44219 7.951 0.6701334
## idv -0.0118636 0.006881 -1.72406 8.503 0.1207518
## mas 0.0023906 0.004027 0.59358 7.902 0.5693905
## uai 0.0035679 0.004689 0.76090 11.534 0.4619919
## lto 0.0006345 0.005962 0.10643 10.672 0.9172124
## ivr -0.0102893 0.007376 -1.39492 10.445 0.1919921
## ln_gdppc 1.1366774 0.479583 2.37014 12.501 0.0346270
## ln_population 0.4699575 0.905670 0.51891 8.409 0.6172031
## internet_users_pct 0.0122864 0.005142 2.38952 7.103 0.0476983
## private_credit_pct_gdp -0.0022078 0.002460 -0.89729 12.783 0.3861497
## trade_open_pct_gdp 0.0016669 0.003362 0.49580 7.779 0.6337353
## secondary_enrollment_pct 0.0048980 0.003129 1.56523 7.200 0.1603237
## iq_avg_m -0.6472420 0.866641 -0.74684 12.003 0.4695370
## ln_gdppc_m -0.0953957 0.581798 -0.16397 14.262 0.8720581
## ln_population_m 0.7554798 0.896454 0.84274 9.267 0.4205899
## internet_users_pct_m 0.0957848 0.047249 2.02723 9.423 0.0718532
## private_credit_pct_gdp_m -0.0053033 0.004331 -1.22454 10.489 0.2475390
## trade_open_pct_gdp_m 0.0010400 0.003927 0.26482 11.682 0.7957639
## secondary_enrollment_pct_m 0.0032321 0.009112 0.35472 7.017 0.7332199
# 5C. [FIX 4] Tách bạch tác động Trong nước (Within) và Giữa các nước (Between) của Thể chế (H1)
institution_effects <- function(m, inst = "iq_avg") {
if (is.null(m)) return(NULL)
b <- stats::coef(m)
im <- paste0(inst, "_m")
if (!all(c(inst, im) %in% names(b))) return(NULL)
V1 <- vcov_country(m)
V2 <- as.matrix(clubSandwich::vcovCR(m, cluster = "individual", type = "CR2"))
one <- function(w, label) {
est <- sum(w * b[names(w)])
se1 <- sqrt(as.numeric(t(w) %*% V1[names(w), names(w)] %*% w))
se2 <- sqrt(as.numeric(t(w) %*% V2[names(w), names(w)] %*% w))
data.frame(
effect = label, estimate = est, clustered_se = se1,
clustered_p = 2 * stats::pnorm(-abs(est / se1)),
CR2_se = se2, CR2_p_normal = 2 * stats::pnorm(-abs(est / se2))
)
}
rbind(
one(stats::setNames(1, inst), "Within-country (Tác động thay đổi theo thời gian = FE)"),
one(stats::setNames(c(1, 1), c(inst, im)), "Between-country (Chênh lệch mức thể chế giữa các nước)")
)
}
inst_effects <- safe_call(institution_effects(cre_full), "Within/between institution effect")
cat("\n================ BẢNG 5D [FIX 4]: TÁC ĐỘNG TRONG NƯỚC VS GIỮA CÁC NƯỚC CỦA THỂ CHẾ (H1) ================\n")
##
## ================ BẢNG 5D [FIX 4]: TÁC ĐỘNG TRONG NƯỚC VS GIỮA CÁC NƯỚC CỦA THỂ CHẾ (H1) ================
print(inst_effects, digits = 4, row.names = FALSE)
## effect estimate clustered_se clustered_p CR2_se CR2_p_normal
## Within-country (Tác động thay đổi theo thời gian = FE) 0.03065 0.2801 0.9129 0.3100 0.9212
## Between-country (Chênh lệch mức thể chế giữa các nước) -0.61659 0.6357 0.3321 0.8499 0.4682
# ==============================================================================
# BƯỚC 6: KIỂM ĐỊNH LỰA CHỌN MÔ HÌNH, KHUYẾT TẬT & ĐỘ NHẠY SAI SỐ CHUẨN ====
# -> Kiểm định F, Breusch-Pagan LM, Hausman, Wald cho các biến trung bình Mundlak,
# tự tương quan (Breusch-Godfrey), phụ thuộc chéo (Pesaran CD), phương sai
# thay đổi (Breusch-Pagan) và so sánh 4 phương pháp tính sai số chuẩn.
# -> Cách chạy: Bôi đen toàn bộ Bước 6 và nhấn Ctrl + Enter.
# ==============================================================================
pool_common <- safe_call(fit_noncre(df, "pooling", with_culture = FALSE), "Pooled common")
re_common <- safe_call(fit_noncre(df, "random", with_culture = FALSE), "RE common")
post_tests <- list()
post_tests[["FE_vs_Pooled_F"]] <- safe_call(plm::pFtest(fe_model, pool_common), "F test: FE vs pooled")
post_tests[["RE_vs_Pooled_LM"]] <- safe_call(plm::plmtest(pool_common, effect = "individual", type = "bp"), "Breusch-Pagan LM")
post_tests[["FE_vs_RE_Hausman"]] <- safe_call(plm::phtest(fe_model, re_common), "Hausman test")
post_tests[["FE_Breusch_Godfrey_serial"]] <- safe_call(plm::pbgtest(fe_model), "Breusch-Godfrey serial correlation")
post_tests[["FE_Pesaran_CD"]] <- safe_call(plm::pcdtest(fe_model, test = "cd"), "Pesaran CD")
pooled_lm_for_bp <- stats::lm(make_formula(culture = character(), add_mundlak = FALSE), data = df)
post_tests[["Pooled_Breusch_Pagan_heteroskedasticity"]] <- safe_call(
lmtest::bptest(pooled_lm_for_bp), "Breusch-Pagan heteroskedasticity"
)
# Kiểm định Wald vững theo cụm quốc gia cho 7 biến trung bình Mundlak (_m = 0)
wald_mundlak <- function(m, means) {
if (is.null(m)) return(NULL)
b <- stats::coef(m)
if (!all(means %in% names(b))) {
stop("Không ước lượng được đầy đủ các biến trung bình Mundlak: ",
paste(setdiff(means, names(b)), collapse = ", "))
}
V <- vcov_country(m)
b0 <- b[means]
V0 <- V[means, means, drop = FALSE]
stat <- as.numeric(crossprod(b0, solve(V0, b0)))
list(statistic = stat, df = length(means),
p_value = stats::pchisq(stat, df = length(means), lower.tail = FALSE))
}
mean_terms <- paste0(c("iq_avg", controls_full), "_m")
mundlak_test <- safe_call(wald_mundlak(cre_full, mean_terms), "Robust joint Wald test of Mundlak means")
post_test_table <- do.call(rbind, lapply(names(post_tests), function(label) {
test <- post_tests[[label]]
if (is.null(test)) return(data.frame(
test = label, statistic = NA_real_, df = NA_character_, p_value = NA_real_,
method = "TEST FAILED", stringsAsFactors = FALSE
))
data.frame(
test = label, statistic = unname(as.numeric(test$statistic[1])),
df = paste(test$parameter, collapse = ", "),
p_value = test$p.value, method = test$method
)
}))
if (!is.null(mundlak_test)) {
post_test_table <- rbind(post_test_table, data.frame(
test = "Mundlak_means_joint_cluster_robust_Wald",
statistic = mundlak_test$statistic,
df = as.character(mundlak_test$df), p_value = mundlak_test$p_value,
method = "Country-clustered HC1 Wald chi-square (7 biến trung bình Mundlak đồng thời = 0)"
))
}
cat("\n================ BẢNG 6A: KIỂM ĐỊNH CHỌN MÔ HÌNH & KHUYẾT TẬT ================\n")
##
## ================ BẢNG 6A: KIỂM ĐỊNH CHỌN MÔ HÌNH & KHUYẾT TẬT ================
print(post_test_table, digits = 4, row.names = FALSE)
## test statistic df p_value
## FE_vs_Pooled_F 103.537 36, 425 1.481e-186
## RE_vs_Pooled_LM 1895.745 1 0.000e+00
## FE_vs_RE_Hausman 1.432 19 1.000e+00
## FE_Breusch_Godfrey_serial 205.649 13 9.437e-37
## FE_Pesaran_CD -1.537 1.243e-01
## Pooled_Breusch_Pagan_heteroskedasticity 126.313 19 7.245e-18
## Mundlak_means_joint_cluster_robust_Wald 10.662 7 1.541e-01
## method
## F test for individual effects
## Lagrange Multiplier Test - (Breusch-Pagan)
## Hausman Test
## Breusch-Godfrey/Wooldridge test for serial correlation in panel models
## Pesaran CD test for cross-sectional dependence in panels
## studentized Breusch-Pagan test
## Country-clustered HC1 Wald chi-square (7 biến trung bình Mundlak đồng thời = 0)
cat("\n[Lưu ý giải thích]:\n",
"- Kiểm định Hausman cổ điển không vững khi có tự tương quan và phương sai thay đổi -> chỉ dùng tham khảo.\n",
"- Hãy dựa vào kiểm định 'Mundlak_means_joint_cluster_robust_Wald' (đã vững theo cụm quốc gia).\n",
"- Tự tương quan (Breusch-Godfrey) và phương sai thay đổi (Breusch-Pagan) đều có p < 0.001,\n",
" do đó việc sử dụng sai số chuẩn cụm theo quốc gia (HC1 / CR2) là hoàn toàn bắt buộc.\n", sep = "")
##
## [Lưu ý giải thích]:
## - Kiểm định Hausman cổ điển không vững khi có tự tương quan và phương sai thay đổi -> chỉ dùng tham khảo.
## - Hãy dựa vào kiểm định 'Mundlak_means_joint_cluster_robust_Wald' (đã vững theo cụm quốc gia).
## - Tự tương quan (Breusch-Godfrey) và phương sai thay đổi (Breusch-Pagan) đều có p < 0.001,
## do đó việc sử dụng sai số chuẩn cụm theo quốc gia (HC1 / CR2) là hoàn toàn bắt buộc.
# 6B. Độ nhạy của sai số chuẩn (HC1 cluster, CR2, Driscoll-Kraay, Double clustering)
se_sensitivity <- do.call(rbind, lapply(
list(FE = fe_model, CRE = cre_full), function(m) {
if (is.null(m)) return(NULL)
do.call(rbind, lapply(c("cluster", "cr2", "dk", "double"), function(inf) {
collect_coefs(m, model_label = if (identical(m, fe_model)) "FE" else "CRE",
inference = inf, keep = c("iq_avg", culture_vars))
}))
}
))
## Warning in sqrt(diag(se)): NaNs produced
## Warning in sqrt(diag(se)): NaNs produced
cat("\n================ BẢNG 6B: ĐỘ NHẠY SAI SỐ CHUẨN (CLUSTER HC1 / CR2 / DK / DOUBLE) ================\n")
##
## ================ BẢNG 6B: ĐỘ NHẠY SAI SỐ CHUẨN (CLUSTER HC1 / CR2 / DK / DOUBLE) ================
cat("(Lưu ý: Không dùng dk và double cho biến văn hoá bất biến theo thời gian khi T = 13)\n")
## (Lưu ý: Không dùng dk và double cho biến văn hoá bất biến theo thời gian khi T = 13)
print(se_sensitivity, digits = 4, row.names = FALSE)
## model inference N countries term estimate std_error statistic df p_value
## FE cluster 481 37 iq_avg 0.0306525 0.2758152 0.11113 NA 9.116e-01
## FE cr2 481 37 iq_avg 0.0306525 0.3100090 0.09888 16.427 9.224e-01
## FE dk 481 37 iq_avg 0.0306525 0.0797103 0.38455 NA 7.008e-01
## FE double 481 37 iq_avg 0.0306525 0.2549579 0.12023 NA 9.044e-01
## CRE cluster 481 37 iq_avg 0.0306525 0.2800917 0.10944 NA 9.129e-01
## CRE cluster 481 37 pdi -0.0035597 0.0060980 -0.58375 NA 5.597e-01
## CRE cluster 481 37 idv -0.0118636 0.0052106 -2.27682 NA 2.327e-02
## CRE cluster 481 37 mas 0.0023906 0.0032043 0.74606 NA 4.560e-01
## CRE cluster 481 37 uai 0.0035679 0.0036266 0.98383 NA 3.257e-01
## CRE cluster 481 37 lto 0.0006345 0.0044844 0.14150 NA 8.875e-01
## CRE cluster 481 37 ivr -0.0102893 0.0054978 -1.87154 NA 6.192e-02
## CRE cr2 481 37 iq_avg 0.0306525 0.3100090 0.09888 16.427 9.224e-01
## CRE cr2 481 37 pdi -0.0035597 0.0080502 -0.44219 7.951 6.701e-01
## CRE cr2 481 37 idv -0.0118636 0.0068812 -1.72406 8.503 1.208e-01
## CRE cr2 481 37 mas 0.0023906 0.0040274 0.59358 7.902 5.694e-01
## CRE cr2 481 37 uai 0.0035679 0.0046891 0.76090 11.534 4.620e-01
## CRE cr2 481 37 lto 0.0006345 0.0059617 0.10643 10.672 9.172e-01
## CRE cr2 481 37 ivr -0.0102893 0.0073762 -1.39492 10.445 1.920e-01
## CRE dk 481 37 iq_avg 0.0306525 0.1053164 0.29105 NA 7.711e-01
## CRE dk 481 37 pdi -0.0035597 0.0073165 -0.48654 NA 6.268e-01
## CRE dk 481 37 idv -0.0118636 0.0097528 -1.21643 NA 2.245e-01
## CRE dk 481 37 mas 0.0023906 0.0086819 0.27535 NA 7.832e-01
## CRE dk 481 37 uai 0.0035679 0.0063302 0.56364 NA 5.733e-01
## CRE dk 481 37 lto 0.0006345 0.0047186 0.13447 NA 8.931e-01
## CRE dk 481 37 ivr -0.0102893 0.0093857 -1.09627 NA 2.735e-01
## CRE double 481 37 iq_avg 0.0306525 0.2630111 0.11654 NA 9.073e-01
## CRE double 481 37 pdi -0.0035597 0.0003154 -11.28586 NA 3.616e-26
## CRE double 481 37 idv -0.0118636 NaN NaN NA NaN
## CRE double 481 37 mas 0.0023906 0.0043954 0.54389 NA 5.868e-01
## CRE double 481 37 uai 0.0035679 0.0012291 2.90283 NA 3.881e-03
## CRE double 481 37 lto 0.0006345 0.0011663 0.54404 NA 5.867e-01
## CRE double 481 37 ivr -0.0102893 0.0052397 -1.96371 NA 5.018e-02
# ==============================================================================
# BƯỚC 7: KIỂM ĐỊNH GIẢ THUYẾT VĂN HOÁ H2a - H2f (TỪNG CHIỀU RIÊNG LẺ) ====
# -> Chạy 6 mô hình CRE riêng cho từng chiều văn hoá Hofstede (PDI, IDV, MAS,
# UAI, LTO, IVR) và kiểm định với p-value CR2 cùng hiệu chỉnh đa kiểm định Holm [FIX 3].
# -> Cách chạy: Bôi đen toàn bộ Bước 7 và nhấn Ctrl + Enter.
# ==============================================================================
culture_models <- lapply(culture_vars, function(v) {
safe_call(fit_cre(df, culture = v), paste0("Culture dimension: ", v))
})
names(culture_models) <- culture_vars
# 7A. Bảng so sánh 6 mô hình văn hoá riêng lẻ
print_model_table(
culture_models,
heading = "BẢNG 7A: KẾT QUẢ 6 MÔ HÌNH CRE RIÊNG CHO TỪNG CHIỀU VĂN HOÁ HOFSTEDE (H2a - H2f)",
keep = c("(Intercept)", "iq_avg", culture_vars, controls_full)
)
##
## ==============================================================================
## BẢNG 7A: KẾT QUẢ 6 MÔ HÌNH CRE RIÊNG CHO TỪNG CHIỀU VĂN HOÁ HOFSTEDE (H2a - H2f)
## Ghi chú: Hệ số ước lượng kèm sai số chuẩn (CLUSTER) trong ngoặc đơn.
## Mức ý nghĩa: * p < 0.10; ** p < 0.05; *** p < 0.01. Đã kiểm soát hiệu ứng năm.
## ==============================================================================
## Variable pdi idv mas uai
## (Intercept) -24.4982*** (4.0969) -25.8216*** (4.2083) -24.1770*** (4.3964) -24.8924*** (4.0227)
## iq_avg 0.0307 (0.2785) 0.0307 (0.2785) 0.0307 (0.2785) 0.0307 (0.2785)
## pdi 0.0021 (0.0056)
## idv -0.0127** (0.0054)
## mas 0.0010 (0.0041)
## uai 0.0088* (0.0045)
## lto
## ivr
## ln_gdppc 1.1367*** (0.4171) 1.1367*** (0.4171) 1.1367*** (0.4171) 1.1367*** (0.4171)
## ln_population 0.4700 (0.7986) 0.4700 (0.7986) 0.4700 (0.7986) 0.4700 (0.7986)
## internet_users_pct 0.0123*** (0.0046) 0.0123*** (0.0046) 0.0123*** (0.0046) 0.0123*** (0.0046)
## private_credit_pct_gdp -0.0022 (0.0022) -0.0022 (0.0022) -0.0022 (0.0022) -0.0022 (0.0022)
## trade_open_pct_gdp 0.0017 (0.0029) 0.0017 (0.0029) 0.0017 (0.0029) 0.0017 (0.0029)
## secondary_enrollment_pct 0.0049* (0.0029) 0.0049* (0.0029) 0.0049* (0.0029) 0.0049* (0.0029)
## Observations (N) 481 481 481 481
## Countries 37 37 37 37
## lto ivr
## -24.5428*** (4.1640) -26.8947*** (4.3423)
## 0.0307 (0.2785) 0.0307 (0.2785)
##
##
##
##
## 0.0097*** (0.0037)
## -0.0126*** (0.0044)
## 1.1367*** (0.4171) 1.1367*** (0.4171)
## 0.4700 (0.7986) 0.4700 (0.7986)
## 0.0123*** (0.0046) 0.0123*** (0.0046)
## -0.0022 (0.0022) -0.0022 (0.0022)
## 0.0017 (0.0029) 0.0017 (0.0029)
## 0.0049* (0.0029) 0.0049* (0.0029)
## 481 481
## 37 37
# 7B. Bảng tổng hợp kiểm định giả thuyết H2a - H2f kèm hiệu chỉnh Holm [FIX 3]
expected_sign <- c(pdi = -1, idv = 1, mas = -1, uai = -1, lto = 1, ivr = 1)
hypothesis_label <- setNames(paste0("H2", letters[1:6]), culture_vars)
culture_results <- do.call(rbind, lapply(culture_vars, function(v) {
m <- culture_models[[v]]
cluster <- get_coefs(m, "cluster")
cr2 <- get_coefs(m, "cr2")
a <- cluster[cluster$term == v, , drop = FALSE]
b <- cr2[cr2$term == v, , drop = FALSE]
p2 <- if (nrow(b)) b$p_value[1] else NA_real_
support <- if (!is.finite(p2)) "CR2 unavailable"
else if (p2 < .05 && sign(a$estimate[1]) == expected_sign[[v]])
"Đúng dấu & CR2 p < 0.05 (trước Holm)"
else "Không đạt ý nghĩa ở mức CR2 p < 0.05"
data.frame(
hypothesis = hypothesis_label[[v]],
variable = v,
expected_direction = ifelse(expected_sign[[v]] > 0, "+", "-"),
beta = a$estimate[1],
clustered_se = a$std_error[1],
clustered_p = a$p_value[1],
CR2_se = if (nrow(b)) b$std_error[1] else NA_real_,
CR2_df = if (nrow(b)) b$df[1] else NA_real_,
CR2_p = p2,
preliminary_result = support
)
}))
culture_results$CR2_p_holm <- stats::p.adjust(culture_results$CR2_p, method = "holm")
culture_results$result_after_holm <- ifelse(
is.finite(culture_results$CR2_p_holm) & culture_results$CR2_p_holm < .05 &
((culture_results$beta > 0 & culture_results$expected_direction == "+") |
(culture_results$beta < 0 & culture_results$expected_direction == "-")),
"Ủng hộ sau hiệu chỉnh Holm", "Không ủng hộ sau hiệu chỉnh Holm"
)
cat("\n================ BẢNG 7B [FIX 3]: TỔNG HỢP KIỂM ĐỊNH GIẢ THUYẾT H2a - H2f ================\n")
##
## ================ BẢNG 7B [FIX 3]: TỔNG HỢP KIỂM ĐỊNH GIẢ THUYẾT H2a - H2f ================
print(culture_results, digits = 4, row.names = FALSE)
## hypothesis variable expected_direction beta clustered_se clustered_p CR2_se CR2_df CR2_p
## H2a pdi - 0.002090 0.005588 0.708495 0.006467 10.463 0.75288
## H2b idv + -0.012688 0.005355 0.018232 0.006343 11.419 0.06981
## H2c mas - 0.001040 0.004082 0.798940 0.004599 9.416 0.82585
## H2d uai - 0.008756 0.004538 0.054269 0.005553 13.385 0.13816
## H2e lto + 0.009712 0.003654 0.008148 0.004285 12.398 0.04204
## H2f ivr + -0.012631 0.004403 0.004317 0.005288 11.962 0.03427
## preliminary_result CR2_p_holm result_after_holm
## Không đạt ý nghĩa ở mức CR2 p < 0.05 1.0000 Không ủng hộ sau hiệu chỉnh Holm
## Không đạt ý nghĩa ở mức CR2 p < 0.05 0.2793 Không ủng hộ sau hiệu chỉnh Holm
## Không đạt ý nghĩa ở mức CR2 p < 0.05 1.0000 Không ủng hộ sau hiệu chỉnh Holm
## Không đạt ý nghĩa ở mức CR2 p < 0.05 0.4145 Không ủng hộ sau hiệu chỉnh Holm
## Đúng dấu & CR2 p < 0.05 (trước Holm) 0.2102 Không ủng hộ sau hiệu chỉnh Holm
## Không đạt ý nghĩa ở mức CR2 p < 0.05 0.2056 Không ủng hộ sau hiệu chỉnh Holm
# ==============================================================================
# BƯỚC 8: PHÂN TÍCH SÁU THÀNH PHẦN THỂ CHẾ WGI (VA, PV, GE, RQ, RL, CC) ====
# -> Thay chỉ số tổng hợp iq_avg bằng từng cấu phần WGI để xem khía cạnh thể chế
# nào có vai trò nổi bật nhất (tính lại trung bình Mundlak cho từng thành phần).
# -> Cách chạy: Bôi đen toàn bộ Bước 8 và nhấn Ctrl + Enter.
# ==============================================================================
wgi_models <- lapply(wgi_vars, function(v) {
safe_call(fit_cre(df, institution = v, culture = culture_vars),
paste0("WGI component: ", v))
})
names(wgi_models) <- wgi_vars
print_model_table(
wgi_models,
heading = "BẢNG 8A: KẾT QUẢ CRE KHI THAY IQ_AVG BẰNG TỪNG THÀNH PHẦN WGI (VA, PV, GE, RQ, RL, CC)",
keep = c(wgi_vars, culture_vars, controls_full)
)
##
## ==============================================================================
## BẢNG 8A: KẾT QUẢ CRE KHI THAY IQ_AVG BẰNG TỪNG THÀNH PHẦN WGI (VA, PV, GE, RQ, RL, CC)
## Ghi chú: Hệ số ước lượng kèm sai số chuẩn (CLUSTER) trong ngoặc đơn.
## Mức ý nghĩa: * p < 0.10; ** p < 0.05; *** p < 0.01. Đã kiểm soát hiệu ứng năm.
## ==============================================================================
## Variable va pv ge rq rl
## va -0.2628 (0.1821)
## pv 0.0720 (0.1277)
## ge 0.0439 (0.1632)
## rq 0.3087* (0.1772)
## rl 0.0608 (0.1998)
## cc
## pdi -0.0027 (0.0051) -0.0015 (0.0052) -0.0026 (0.0058) -0.0035 (0.0059) -0.0039 (0.0068)
## idv -0.0065* (0.0035) -0.0100* (0.0052) -0.0131** (0.0064) -0.0123** (0.0055) -0.0119** (0.0053)
## mas 0.0023 (0.0033) 0.0000 (0.0044) 0.0017 (0.0035) 0.0026 (0.0033) 0.0020 (0.0033)
## uai 0.0092* (0.0049) 0.0068 (0.0045) 0.0031 (0.0034) 0.0025 (0.0035) 0.0047 (0.0040)
## lto -0.0004 (0.0041) 0.0017 (0.0045) 0.0019 (0.0046) 0.0003 (0.0045) 0.0008 (0.0045)
## ivr -0.0080* (0.0045) 0.0006 (0.0074) -0.0078 (0.0051) -0.0088* (0.0050) -0.0097* (0.0058)
## ln_gdppc 1.1645** (0.4569) 1.1258*** (0.4204) 1.1201*** (0.4089) 1.1210*** (0.3941) 1.1240*** (0.4148)
## ln_population 0.5570 (0.7536) 0.4676 (0.7362) 0.4947 (0.8070) 0.7337 (0.7511) 0.5010 (0.7862)
## internet_users_pct 0.0114** (0.0048) 0.0114** (0.0044) 0.0124*** (0.0046) 0.0125*** (0.0045) 0.0124*** (0.0046)
## private_credit_pct_gdp -0.0028 (0.0022) -0.0021 (0.0022) -0.0022 (0.0022) -0.0022 (0.0021) -0.0022 (0.0022)
## trade_open_pct_gdp 0.0020 (0.0028) 0.0016 (0.0029) 0.0017 (0.0028) 0.0013 (0.0026) 0.0017 (0.0029)
## secondary_enrollment_pct 0.0042 (0.0029) 0.0049* (0.0029) 0.0049* (0.0029) 0.0054* (0.0029) 0.0050* (0.0029)
## Observations (N) 481 481 481 481 481
## Countries 37 37 37 37 37
## cc
##
##
##
##
##
## -0.0793 (0.1842)
## -0.0042 (0.0062)
## -0.0120** (0.0054)
## 0.0017 (0.0034)
## 0.0019 (0.0035)
## 0.0006 (0.0046)
## -0.0100* (0.0052)
## 1.1620*** (0.4370)
## 0.3857 (0.7793)
## 0.0122*** (0.0046)
## -0.0022 (0.0022)
## 0.0016 (0.0030)
## 0.0049* (0.0029)
## 481
## 37
wgi_results_cr2 <- do.call(rbind, lapply(wgi_vars, function(v) {
collect_coefs(wgi_models[[v]], paste0("WGI_", v), "cr2", keep = v)
}))
cat("\n================ BẢNG 8B: CHI TIẾT 6 THÀNH PHẦN WGI THEO SAI SỐ CHUẨN CR2 ================\n")
##
## ================ BẢNG 8B: CHI TIẾT 6 THÀNH PHẦN WGI THEO SAI SỐ CHUẨN CR2 ================
print(wgi_results_cr2, digits = 4, row.names = FALSE)
## model inference N countries term estimate std_error statistic df p_value
## WGI_va cr2 481 37 va -0.26277 0.2001 -1.3131 9.568 0.2198
## WGI_pv cr2 481 37 pv 0.07197 0.1371 0.5248 16.730 0.6066
## WGI_ge cr2 481 37 ge 0.04389 0.1809 0.2427 20.313 0.8107
## WGI_rq cr2 481 37 rq 0.30869 0.1916 1.6108 21.279 0.1220
## WGI_rl cr2 481 37 rl 0.06077 0.2220 0.2737 13.086 0.7886
## WGI_cc cr2 481 37 cc -0.07926 0.1980 -0.4003 22.207 0.6927
# ==============================================================================
# BƯỚC 9: HỆ SỐ CHUẨN HOÁ & DIỄN GIẢI Ý NGHĨA KINH TẾ [FIX 5] ====
# -> Tính hệ số chuẩn hoá (Standardized Beta) và mức thay đổi % của doanh số TMĐT
# khi biến độc lập thay đổi 1 độ lệch chuẩn trong nước (within-country SD)
# hoặc 1 đơn vị đo lường thực tế (+1 điểm WGI, +10 điểm Hofstede, +1 điểm %).
# -> Cách chạy: Bôi đen toàn bộ Bước 9 và nhấn Ctrl + Enter.
# ==============================================================================
if (!is.null(cre_full)) {
source_coefs <- get_coefs(cre_full, "cluster")
substantive <- c("iq_avg", culture_vars, controls_full)
source_coefs <- source_coefs[source_coefs$term %in% substantive, , drop = FALSE]
source_coefs$sd_x <- vapply(source_coefs$term, function(v) stats::sd(df[[v]]), numeric(1))
source_coefs$sd_y <- stats::sd(df$ln_ecom)
source_coefs$beta_standardized <- source_coefs$estimate * source_coefs$sd_x / source_coefs$sd_y
# [FIX 5] Độ lệch chuẩn trong từng quốc gia (Within-country SD) cho các biến thay đổi theo thời gian
within_sd <- function(v) stats::sd(df[[v]] - stats::ave(df[[v]], df$iso3))
source_coefs$sd_x_within <- vapply(source_coefs$term, function(v) {
if (v %in% culture_vars) NA_real_ else within_sd(v)
}, numeric(1))
source_coefs$ecom_pct_change_per_within_sd <- 100 *
(exp(source_coefs$estimate * source_coefs$sd_x_within) - 1)
source_coefs$change_in_x <- ifelse(source_coefs$term %in% culture_vars, 10, 1)
source_coefs$ecom_percent_change <- 100 *
(exp(source_coefs$estimate * source_coefs$change_in_x) - 1)
source_coefs$interpretation_unit <- ifelse(
source_coefs$term %in% culture_vars, "+10 điểm Hofstede",
ifelse(source_coefs$term %in% c("ln_gdppc", "ln_population"),
"+1 đơn vị log (Độ co giãn = estimate)",
ifelse(source_coefs$term %in% c("internet_users_pct",
"private_credit_pct_gdp",
"trade_open_pct_gdp",
"secondary_enrollment_pct"),
"+1 điểm phần trăm (%)", "+1 đơn vị WGI"))
)
cat("\n================ BẢNG 9: HỆ SỐ CHUẨN HOÁ & Ý NGHĨA KINH TẾ (MÔ HÌNH CRE CHÍNH) ================\n")
print(source_coefs[, c("term", "estimate", "beta_standardized", "sd_x_within",
"ecom_pct_change_per_within_sd", "ecom_percent_change",
"interpretation_unit")], digits = 4, row.names = FALSE)
}
##
## ================ BẢNG 9: HỆ SỐ CHUẨN HOÁ & Ý NGHĨA KINH TẾ (MÔ HÌNH CRE CHÍNH) ================
## term estimate beta_standardized sd_x_within ecom_pct_change_per_within_sd ecom_percent_change
## iq_avg 0.0306525 0.010392 0.08340 0.256 3.1127
## pdi -0.0035597 -0.038807 NA NA -3.4971
## idv -0.0118636 -0.139254 NA NA -11.1869
## mas 0.0023906 0.026719 NA NA 2.4194
## uai 0.0035679 0.043977 NA NA 3.6323
## lto 0.0006345 0.006961 NA NA 0.6365
## ivr -0.0102893 -0.102649 NA NA -9.7776
## ln_gdppc 1.1366774 0.468370 0.08909 10.658 211.6397
## ln_population 0.4699575 0.394800 0.03363 1.593 59.9926
## internet_users_pct 0.0122864 0.087872 8.72162 11.311 1.2362
## private_credit_pct_gdp -0.0022078 -0.049253 11.40214 -2.486 -0.2205
## trade_open_pct_gdp 0.0016669 0.071003 11.27172 1.897 0.1668
## secondary_enrollment_pct 0.0048980 0.040400 5.18632 2.573 0.4910
## interpretation_unit
## +1 đơn vị WGI
## +10 điểm Hofstede
## +10 điểm Hofstede
## +10 điểm Hofstede
## +10 điểm Hofstede
## +10 điểm Hofstede
## +10 điểm Hofstede
## +1 đơn vị log (Độ co giãn = estimate)
## +1 đơn vị log (Độ co giãn = estimate)
## +1 điểm phần trăm (%)
## +1 điểm phần trăm (%)
## +1 điểm phần trăm (%)
## +1 điểm phần trăm (%)
# ==============================================================================
# BƯỚC 10: KIỂM ĐỊNH ĐỘ BỀN (ROBUSTNESS CHECKS - 13 ĐẶC TẢ) [FIX 2, 6, 8] ====
# -> Chạy 13 đặc tả độ bền (bỏ biến giáo dục, chỉ giữ quan sát gốc, mẫu thận trọng,
# bỏ năm COVID 2020-2021, bỏ 3 nền kinh tế lớn USA-CHN-JPN, TMĐT bình quân đầu
# người, trễ 1 năm của thể chế, giai đoạn 2016-2024, loại ô nội suy/bước nhảy...).
# -> Cách chạy: Bôi đen toàn bộ Bước 10 và nhấn Ctrl + Enter (mất khoảng 10-20 giây).
# ==============================================================================
df_percap <- df |>
dplyr::mutate(ecom_pc_usd = 1e6 * ecom_musd / population_total,
ln_ecom_pc = log(ecom_pc_usd))
df_lag <- df |>
dplyr::group_by(iso3) |>
dplyr::arrange(year, .by_group = TRUE) |>
dplyr::mutate(
iq_avg_l1 = dplyr::if_else(
year - dplyr::lag(year) == 1L, dplyr::lag(iq_avg), NA_real_
)
) |>
dplyr::ungroup() |>
dplyr::filter(!is.na(iq_avg_l1)) |>
as.data.frame()
robust_specifications <- list(
Baseline_full = list(data = df, institution = "iq_avg", controls = controls_full, outcome = "ln_ecom"),
Without_education = list(data = df, institution = "iq_avg",
controls = setdiff(controls_full, "secondary_enrollment_pct"), outcome = "ln_ecom"),
Observed_only = list(data = subset(df, observed_only_sample == 1),
institution = "iq_avg", controls = controls_full, outcome = "ln_ecom"),
Conservative_sample = list(data = subset(df, conservative_sample == 1),
institution = "iq_avg", controls = controls_full, outcome = "ln_ecom"),
No_high_risk = list(data = subset(df, high_risk_cell_count == 0),
institution = "iq_avg", controls = controls_full, outcome = "ln_ecom"),
Without_COVID_2020_2021 = list(data = subset(df, !(year %in% c(2020, 2021))),
institution = "iq_avg", controls = controls_full, outcome = "ln_ecom"),
Without_USA_CHN_JPN = list(data = subset(df, !(iso3 %in% c("USA", "CHN", "JPN"))),
institution = "iq_avg", controls = controls_full, outcome = "ln_ecom"),
Ecom_per_capita = list(data = df_percap, institution = "iq_avg",
controls = setdiff(controls_full, "ln_population"), outcome = "ln_ecom_pc"),
Institution_lag_1y = list(data = df_lag, institution = "iq_avg_l1",
controls = controls_full, outcome = "ln_ecom"),
Period_2016_2024 = list(data = subset(df, year >= 2016), institution = "iq_avg",
controls = controls_full, outcome = "ln_ecom")
)
if ("ecom_imputed" %in% names(df)) {
robust_specifications[["Ecom_observed_only"]] <- list(
data = subset(df, ecom_imputed == 0), institution = "iq_avg",
controls = controls_full, outcome = "ln_ecom"
)
}
if (all(c("ecom_imputed", "ecom_unverified") %in% names(df))) {
robust_specifications[["Ecom_source_values_only"]] <- list(
data = subset(df, ecom_imputed == 0 & ecom_unverified == 0),
institution = "iq_avg", controls = controls_full, outcome = "ln_ecom"
)
}
if ("ecom_jump_country" %in% names(df)) {
robust_specifications[["Without_ecom_jump_countries"]] <- list(
data = subset(df, ecom_jump_country == 0), institution = "iq_avg",
controls = controls_full, outcome = "ln_ecom"
)
}
robust_joint_models <- list()
robust_details <- list()
robust_samples <- list()
for (label in names(robust_specifications)) {
sp <- robust_specifications[[label]]
robust_samples[[label]] <- data.frame(
specification = label,
sample_rows = nrow(sp$data),
sample_countries = dplyr::n_distinct(sp$data$iso3),
sample_years = dplyr::n_distinct(sp$data$year)
)
joint <- safe_call(fit_cre(sp$data, institution = sp$institution,
culture = culture_vars, controls = sp$controls,
y = sp$outcome), paste0(label, " [joint 6 culture]"))
robust_joint_models[[label]] <- joint
for (inf in c("cluster", "cr2")) {
robust_details[[paste(label, "joint", inf, sep = "_")]] <- collect_coefs(
joint, paste0(label, "_joint"), inf, keep = c(sp$institution, culture_vars)
)
}
# Đối chiếu hệ số thể chế với FE trên cùng mẫu con
fe_same <- safe_call({
need <- c("iso3", "year", sp$outcome, sp$institution, sp$controls)
dd <- sp$data[stats::complete.cases(sp$data[need]), , drop = FALSE]
plm::plm(make_formula(sp$outcome, sp$institution, character(), sp$controls),
data = plm::pdata.frame(dd, index = c("iso3", "year"), drop.index = FALSE),
model = "within", effect = "individual")
}, paste0(label, " [FE check]"))
for (inf in c("cluster", "cr2")) {
robust_details[[paste(label, "FE", inf, sep = "_")]] <- collect_coefs(
fe_same, paste0(label, "_FE_check"), inf, keep = sp$institution
)
}
}
## [Lưu ý FIX 2] Swamy-Arora không khả dụng trên mẫu không cân bằng; tự động chuyển sang Wallace-Hussain.
## [Lưu ý FIX 2] Swamy-Arora không khả dụng trên mẫu không cân bằng; tự động chuyển sang Wallace-Hussain.
## [Lưu ý FIX 2] Swamy-Arora không khả dụng trên mẫu không cân bằng; tự động chuyển sang Wallace-Hussain.
## [Lưu ý FIX 2] Swamy-Arora không khả dụng trên mẫu không cân bằng; tự động chuyển sang Wallace-Hussain.
## [Lưu ý FIX 2] Swamy-Arora không khả dụng trên mẫu không cân bằng; tự động chuyển sang Wallace-Hussain.
robust_samples_df <- do.call(rbind, robust_samples)
cat("\n================ BẢNG 10A: CỠ MẪU CỦA CÁC ĐẶC TẢ KIỂM ĐỊNH ĐỘ BỀN ================\n")
##
## ================ BẢNG 10A: CỠ MẪU CỦA CÁC ĐẶC TẢ KIỂM ĐỊNH ĐỘ BỀN ================
print(robust_samples_df, row.names = FALSE)
## specification sample_rows sample_countries sample_years
## Baseline_full 481 37 13
## Without_education 481 37 13
## Observed_only 420 36 13
## Conservative_sample 449 36 13
## No_high_risk 449 36 13
## Without_COVID_2020_2021 407 37 11
## Without_USA_CHN_JPN 442 34 13
## Ecom_per_capita 481 37 13
## Institution_lag_1y 444 37 12
## Period_2016_2024 333 37 9
## Ecom_observed_only 444 37 13
## Ecom_source_values_only 374 36 13
## Without_ecom_jump_countries 338 26 13
# Hiển thị bảng kết quả các mô hình độ bền trên Console
print_model_table(
robust_joint_models,
heading = "BẢNG 10B: KẾT QUẢ KIỂM ĐỊNH ĐỘ BỀN (MÔ HÌNH ĐẦY ĐỦ 6 CHIỀU VĂN HOÁ - HC1)",
keep = c("iq_avg", "iq_avg_l1", culture_vars)
)
##
## ==============================================================================
## BẢNG 10B: KẾT QUẢ KIỂM ĐỊNH ĐỘ BỀN (MÔ HÌNH ĐẦY ĐỦ 6 CHIỀU VĂN HOÁ - HC1)
## Ghi chú: Hệ số ước lượng kèm sai số chuẩn (CLUSTER) trong ngoặc đơn.
## Mức ý nghĩa: * p < 0.10; ** p < 0.05; *** p < 0.01. Đã kiểm soát hiệu ứng năm.
## ==============================================================================
## Variable Baseline_full Without_education Observed_only Conservative_sample No_high_risk
## iq_avg 0.0307 (0.2801) 0.0060 (0.2805) 0.0195 (0.2787) -0.0645 (0.2778) -0.0645 (0.2778)
## iq_avg_l1
## pdi -0.0036 (0.0061) -0.0028 (0.0056) -0.0030 (0.0068) -0.0027 (0.0065) -0.0027 (0.0065)
## idv -0.0119** (0.0052) -0.0094** (0.0043) -0.0099** (0.0048) -0.0101* (0.0052) -0.0101* (0.0052)
## mas 0.0024 (0.0032) 0.0011 (0.0031) 0.0018 (0.0033) 0.0021 (0.0032) 0.0021 (0.0032)
## uai 0.0036 (0.0036) 0.0049 (0.0039) 0.0024 (0.0038) 0.0012 (0.0037) 0.0012 (0.0037)
## lto 0.0006 (0.0045) 0.0010 (0.0044) 0.0046 (0.0051) 0.0043 (0.0049) 0.0043 (0.0049)
## ivr -0.0103* (0.0055) -0.0091* (0.0048) -0.0080 (0.0057) -0.0090 (0.0058) -0.0090 (0.0058)
## Observations (N) 481 481 420 449 449
## Countries 37 37 36 36 36
## Without_COVID_2020_2021 Without_USA_CHN_JPN Ecom_per_capita Institution_lag_1y Period_2016_2024
## 0.0318 (0.2891) -0.0945 (0.2842) 0.0656 (0.2591) -0.1278 (0.2426)
## 0.0023 (0.2729)
## -0.0035 (0.0065) -0.0030 (0.0064) -0.0009 (0.0063) -0.0034 (0.0061) -0.0024 (0.0053)
## -0.0118** (0.0054) -0.0120** (0.0054) -0.0094* (0.0049) -0.0115** (0.0050) -0.0109** (0.0050)
## 0.0026 (0.0033) 0.0022 (0.0040) 0.0054* (0.0030) 0.0023 (0.0031) 0.0029 (0.0031)
## 0.0035 (0.0037) 0.0056 (0.0044) 0.0000 (0.0040) 0.0029 (0.0037) 0.0025 (0.0036)
## -0.0002 (0.0047) 0.0005 (0.0051) 0.0045 (0.0044) 0.0010 (0.0045) 0.0043 (0.0044)
## -0.0104* (0.0059) -0.0068 (0.0056) -0.0087 (0.0057) -0.0107** (0.0054) -0.0095** (0.0047)
## 407 442 481 444 333
## 37 34 37 37 37
## Ecom_observed_only Ecom_source_values_only Without_ecom_jump_countries
## 0.1519 (0.2664) 0.2284 (0.2898) 0.0950 (0.2302)
##
## -0.0052 (0.0061) -0.0021 (0.0057) 0.0025 (0.0049)
## -0.0141*** (0.0053) -0.0094** (0.0043) -0.0143* (0.0079)
## 0.0030 (0.0032) 0.0018 (0.0032) 0.0040 (0.0040)
## 0.0039 (0.0039) 0.0028 (0.0045) 0.0002 (0.0035)
## -0.0011 (0.0046) 0.0042 (0.0051) -0.0015 (0.0043)
## -0.0085 (0.0052) -0.0038 (0.0062) -0.0086 (0.0070)
## 444 374 338
## 37 36 26
# Tổng hợp riêng kết quả kiểm định độ bền theo CR2 để tra cứu nhanh
robust_results <- dplyr::bind_rows(robust_details)
robust_cr2_summary <- subset(robust_results, inference == "cr2" & grepl("_joint$|_FE_check$", model))
cat("\n================ BẢNG 10C: HỆ SỐ THỂ CHẾ & VĂN HOÁ QUA CÁC KIỂM ĐỊNH ĐỘ BỀN (CR2) ================\n")
##
## ================ BẢNG 10C: HỆ SỐ THỂ CHẾ & VĂN HOÁ QUA CÁC KIỂM ĐỊNH ĐỘ BỀN (CR2) ================
print(robust_cr2_summary[, c("model", "N", "countries", "term", "estimate", "std_error", "df", "p_value")],
digits = 4, row.names = FALSE)
## model N countries term estimate std_error df p_value
## Baseline_full_joint 481 37 iq_avg 0.03065253 0.310009 16.427 0.92243
## Baseline_full_joint 481 37 pdi -0.00355973 0.008050 7.951 0.67013
## Baseline_full_joint 481 37 idv -0.01186360 0.006881 8.503 0.12075
## Baseline_full_joint 481 37 mas 0.00239060 0.004027 7.902 0.56939
## Baseline_full_joint 481 37 uai 0.00356793 0.004689 11.534 0.46199
## Baseline_full_joint 481 37 lto 0.00063453 0.005962 10.672 0.91721
## Baseline_full_joint 481 37 ivr -0.01028930 0.007376 10.445 0.19199
## Baseline_full_FE_check 481 37 iq_avg 0.03065253 0.310009 16.427 0.92243
## Without_education_joint 481 37 iq_avg 0.00602201 0.311282 16.403 0.98480
## Without_education_joint 481 37 pdi -0.00279604 0.007189 8.208 0.70721
## Without_education_joint 481 37 idv -0.00940816 0.005575 10.736 0.12027
## Without_education_joint 481 37 mas 0.00109009 0.003770 7.496 0.78030
## Without_education_joint 481 37 uai 0.00485541 0.005024 12.561 0.35203
## Without_education_joint 481 37 lto 0.00098004 0.005689 10.731 0.86643
## Without_education_joint 481 37 ivr -0.00908780 0.006203 10.057 0.17343
## Without_education_FE_check 481 37 iq_avg 0.00602201 0.311282 16.403 0.98480
## Observed_only_joint 420 36 iq_avg 0.01949649 0.312028 15.218 0.95099
## Observed_only_joint 420 36 pdi -0.00304646 0.008721 8.171 0.73568
## Observed_only_joint 420 36 idv -0.00992635 0.006238 8.237 0.14910
## Observed_only_joint 420 36 mas 0.00184294 0.004203 8.000 0.67262
## Observed_only_joint 420 36 uai 0.00235143 0.005030 10.864 0.64937
## Observed_only_joint 420 36 lto 0.00463339 0.007281 9.601 0.53941
## Observed_only_joint 420 36 ivr -0.00802859 0.007547 10.012 0.31237
## Observed_only_FE_check 420 36 iq_avg 0.02024706 0.311394 15.217 0.94900
## Conservative_sample_joint 449 36 iq_avg -0.06453848 0.312056 15.848 0.83879
## Conservative_sample_joint 449 36 pdi -0.00269848 0.008439 8.122 0.75721
## Conservative_sample_joint 449 36 idv -0.01014042 0.006846 8.201 0.17594
## Conservative_sample_joint 449 36 mas 0.00206432 0.004038 7.860 0.62324
## Conservative_sample_joint 449 36 uai 0.00119604 0.004786 10.986 0.80726
## Conservative_sample_joint 449 36 lto 0.00429181 0.006943 9.513 0.55100
## Conservative_sample_joint 449 36 ivr -0.00897837 0.007807 10.014 0.27688
## Conservative_sample_FE_check 449 36 iq_avg -0.06410890 0.311725 15.848 0.83968
## No_high_risk_joint 449 36 iq_avg -0.06453848 0.312056 15.848 0.83879
## No_high_risk_joint 449 36 pdi -0.00269848 0.008439 8.122 0.75721
## No_high_risk_joint 449 36 idv -0.01014042 0.006846 8.201 0.17594
## No_high_risk_joint 449 36 mas 0.00206432 0.004038 7.860 0.62324
## No_high_risk_joint 449 36 uai 0.00119604 0.004786 10.986 0.80726
## No_high_risk_joint 449 36 lto 0.00429181 0.006943 9.513 0.55100
## No_high_risk_joint 449 36 ivr -0.00897837 0.007807 10.014 0.27688
## No_high_risk_FE_check 449 36 iq_avg -0.06410890 0.311725 15.848 0.83968
## Without_COVID_2020_2021_joint 407 37 iq_avg 0.03178328 0.317640 17.154 0.92146
## Without_COVID_2020_2021_joint 407 37 pdi -0.00351390 0.008482 8.002 0.68956
## Without_COVID_2020_2021_joint 407 37 idv -0.01182007 0.007044 8.649 0.12902
## Without_COVID_2020_2021_joint 407 37 mas 0.00262953 0.004207 7.802 0.54982
## Without_COVID_2020_2021_joint 407 37 uai 0.00350914 0.004826 11.693 0.48144
## Without_COVID_2020_2021_joint 407 37 lto -0.00022171 0.006356 11.056 0.97279
## Without_COVID_2020_2021_joint 407 37 ivr -0.01037652 0.007913 10.309 0.21819
## Without_COVID_2020_2021_FE_check 407 37 iq_avg 0.03178328 0.317640 17.154 0.92146
## Without_USA_CHN_JPN_joint 442 34 iq_avg -0.09448836 0.315163 14.948 0.76845
## Without_USA_CHN_JPN_joint 442 34 pdi -0.00301010 0.008625 7.966 0.73613
## Without_USA_CHN_JPN_joint 442 34 idv -0.01203595 0.007360 8.324 0.13916
## Without_USA_CHN_JPN_joint 442 34 mas 0.00218291 0.005168 7.323 0.68486
## Without_USA_CHN_JPN_joint 442 34 uai 0.00563633 0.006015 8.976 0.37327
## Without_USA_CHN_JPN_joint 442 34 lto 0.00054095 0.006642 8.583 0.93696
## Without_USA_CHN_JPN_joint 442 34 ivr -0.00675191 0.007577 8.051 0.39873
## Without_USA_CHN_JPN_FE_check 442 34 iq_avg -0.09448836 0.315163 14.948 0.76845
## Ecom_per_capita_joint 481 37 iq_avg 0.06563699 0.284263 16.885 0.82017
## Ecom_per_capita_joint 481 37 pdi -0.00086644 0.008293 7.653 0.91948
## Ecom_per_capita_joint 481 37 idv -0.00944629 0.006489 9.106 0.17909
## Ecom_per_capita_joint 481 37 mas 0.00539844 0.003764 10.455 0.18069
## Ecom_per_capita_joint 481 37 uai 0.00000329 0.005044 8.251 0.99950
## Ecom_per_capita_joint 481 37 lto 0.00446703 0.005532 11.233 0.43617
## Ecom_per_capita_joint 481 37 ivr -0.00872354 0.007398 10.883 0.26345
## Ecom_per_capita_FE_check 481 37 iq_avg 0.06563699 0.284263 16.885 0.82017
## Institution_lag_1y_joint 444 37 iq_avg_l1 0.00232010 0.299788 14.980 0.99393
## Institution_lag_1y_joint 444 37 pdi -0.00340672 0.007973 7.868 0.68064
## Institution_lag_1y_joint 444 37 idv -0.01145468 0.006585 8.546 0.11775
## Institution_lag_1y_joint 444 37 mas 0.00230501 0.003950 7.888 0.57586
## Institution_lag_1y_joint 444 37 uai 0.00292602 0.004695 11.216 0.54560
## Institution_lag_1y_joint 444 37 lto 0.00103084 0.005939 10.582 0.86547
## Institution_lag_1y_joint 444 37 ivr -0.01074806 0.007165 10.495 0.16304
## Institution_lag_1y_FE_check 444 37 iq_avg_l1 0.00232010 0.299788 14.980 0.99393
## Period_2016_2024_joint 333 37 iq_avg -0.12782641 0.277302 14.146 0.65183
## Period_2016_2024_joint 333 37 pdi -0.00235110 0.006942 7.346 0.74432
## Period_2016_2024_joint 333 37 idv -0.01094492 0.006495 8.119 0.12991
## Period_2016_2024_joint 333 37 mas 0.00286322 0.003879 8.118 0.48122
## Period_2016_2024_joint 333 37 uai 0.00248421 0.004595 10.531 0.60001
## Period_2016_2024_joint 333 37 lto 0.00425871 0.005688 9.833 0.47155
## Period_2016_2024_joint 333 37 ivr -0.00949201 0.006047 10.925 0.14498
## Period_2016_2024_FE_check 333 37 iq_avg -0.12782641 0.277302 14.146 0.65183
## Ecom_observed_only_joint 444 37 iq_avg 0.15193771 0.288388 17.406 0.60495
## Ecom_observed_only_joint 444 37 pdi -0.00522797 0.007949 7.994 0.52924
## Ecom_observed_only_joint 444 37 idv -0.01406013 0.007028 8.782 0.07728
## Ecom_observed_only_joint 444 37 mas 0.00304687 0.003991 7.889 0.46741
## Ecom_observed_only_joint 444 37 uai 0.00392693 0.005031 11.809 0.45039
## Ecom_observed_only_joint 444 37 lto -0.00111318 0.006172 10.333 0.86035
## Ecom_observed_only_joint 444 37 ivr -0.00847149 0.006924 10.319 0.24835
## Ecom_observed_only_FE_check 444 37 iq_avg 0.15116091 0.288558 17.406 0.60699
## Ecom_source_values_only_joint 374 36 iq_avg 0.22835102 0.327851 14.447 0.49717
## Ecom_source_values_only_joint 374 36 pdi -0.00212213 0.007542 7.318 0.78621
## Ecom_source_values_only_joint 374 36 idv -0.00937591 0.005553 8.778 0.12645
## Ecom_source_values_only_joint 374 36 mas 0.00179775 0.003965 7.520 0.66304
## Ecom_source_values_only_joint 374 36 uai 0.00284782 0.005627 11.206 0.62259
## Ecom_source_values_only_joint 374 36 lto 0.00416364 0.006593 11.031 0.54058
## Ecom_source_values_only_joint 374 36 ivr -0.00383110 0.008102 10.826 0.64570
## Ecom_source_values_only_FE_check 374 36 iq_avg 0.22926980 0.328030 14.447 0.49570
## Without_ecom_jump_countries_joint 338 26 iq_avg 0.09500651 0.263265 11.555 0.72470
## Without_ecom_jump_countries_joint 338 26 pdi 0.00246785 0.007012 4.503 0.74075
## Without_ecom_jump_countries_joint 338 26 idv -0.01425277 0.009617 8.312 0.17520
## Without_ecom_jump_countries_joint 338 26 mas 0.00401523 0.005640 5.451 0.50583
## Without_ecom_jump_countries_joint 338 26 uai 0.00023326 0.004714 6.308 0.96207
## Without_ecom_jump_countries_joint 338 26 lto -0.00148212 0.005493 7.671 0.79443
## Without_ecom_jump_countries_joint 338 26 ivr -0.00862477 0.009777 6.326 0.40997
## Without_ecom_jump_countries_FE_check 338 26 iq_avg 0.09500651 0.263265 11.555 0.72470
# ==============================================================================
# BƯỚC 11: PHÂN MẪU CÁC QUỐC GIA TƯƠNG ĐỒNG VỚI VIỆT NAM (TUỲ CHỌN) ====
# -> Việt Nam không nằm trong 37 nước của panel. Nếu bạn đã có điểm WGI trung bình
# (iq_avg_m giai đoạn 2012-2024) và 6 điểm Hofstede của Việt Nam, hãy điền
# vào dòng `vietnam_ref` bên dưới (thay các chữ NA_real_ bằng con số thực tế)
# hoặc đặt file `vietnam_reference.csv` cùng thư mục rồi chạy Bước 11.
# -> Cách chạy: Bôi đen toàn bộ Bước 11 và nhấn Ctrl + Enter.
# ==============================================================================
# Bạn có thể điền trực tiếp 7 con số của Việt Nam vào đây (nếu không dùng file CSV):
vietnam_ref <- data.frame(
iso3 = "VNM",
iq_avg_m = NA_real_, # Điền trung bình WGI của Việt Nam (ví dụ: -0.42)
pdi = NA_real_, # Điền PDI của Việt Nam (ví dụ: 70)
idv = NA_real_, # Điền IDV của Việt Nam (ví dụ: 20)
mas = NA_real_, # Điền MAS của Việt Nam (ví dụ: 40)
uai = NA_real_, # Điền UAI của Việt Nam (ví dụ: 30)
lto = NA_real_, # Điền LTO của Việt Nam (ví dụ: 57)
ivr = NA_real_ # Điền IVR của Việt Nam (ví dụ: 35)
)
reference_path <- file.path(base_dir, "vietnam_reference.csv")
if (file.exists(reference_path)) {
log_note("Đang đọc thông số tham chiếu Việt Nam từ file vietnam_reference.csv...")
vietnam_ref <- safe_call(readr::read_csv(reference_path, show_col_types = FALSE),
"Vietnam reference CSV")
}
required_ref <- c("iq_avg_m", culture_vars)
valid_ref <- !is.null(vietnam_ref) && nrow(vietnam_ref) == 1L &&
all(required_ref %in% names(vietnam_ref)) &&
all(vapply(vietnam_ref[required_ref], function(z) is.numeric(z) &&
length(z) == 1L && is.finite(z), logical(1)))
if (!valid_ref) {
cat("\n[BƯỚC 11 - BỎ QUA]: Chưa có đủ 7 thông số tham chiếu của Việt Nam (iq_avg_m, pdi, idv, mas, uai, lto, ivr).\n")
cat("-> Để chạy Bước 11, bạn chỉ cần thay các giá trị NA_real_ trong biến `vietnam_ref` ở đầu Bước 11\n")
cat(" bằng số liệu WGI & Hofstede của Việt Nam rồi bôi đen Bước 11 nhấn Ctrl + Enter!\n")
} else {
country_data <- df |>
dplyr::group_by(iso3, country) |>
dplyr::summarise(
iq_avg_m = mean(iq_avg),
dplyr::across(dplyr::all_of(culture_vars), dplyr::first),
.groups = "drop"
)
country_data <- country_data[country_data$iso3 != "VNM", , drop = FALSE]
attribute_matrix <- as.matrix(rbind(country_data[required_ref], vietnam_ref[required_ref]))
storage.mode(attribute_matrix) <- "double"
attr_mean <- colMeans(attribute_matrix)
attr_sd <- apply(attribute_matrix, 2, sd)
z <- scale(attribute_matrix, center = attr_mean, scale = attr_sd)
ref_z <- z[nrow(z), , drop = FALSE]
distance <- sqrt(rowSums((z[-nrow(z), , drop = FALSE] -
matrix(as.numeric(ref_z), nrow = nrow(z) - 1L,
ncol = ncol(z), byrow = TRUE))^2))
country_data$distance_to_VNM <- round(as.numeric(distance), 4)
country_data <- country_data[order(country_data$distance_to_VNM, country_data$iso3), , drop = FALSE]
country_data$similarity_rank <- seq_len(nrow(country_data))
cat("\n================ BẢNG 11A: XẾP HẠNG MỨC ĐỘ TƯƠNG ĐỒNG VỚI VIỆT NAM ================\n")
print(as.data.frame(country_data), digits = 4, row.names = FALSE)
matched_results <- list()
for (k in c(10L, 12L, 15L)) {
eligible <- head(country_data$iso3, k)
matched <- df[df$iso3 %in% eligible, , drop = FALSE]
for (cv in culture_vars) {
m <- safe_call(fit_cre(matched, culture = cv),
paste0("Vietnam match ", k, " countries / ", cv))
for (inf in c("cluster", "cr2")) {
matched_results[[paste(k, cv, inf, sep = "_")]] <- collect_coefs(
m, paste0("VNM_nearest", k, "_", cv), inf, keep = c("iq_avg", cv)
)
}
}
}
all_matched <- dplyr::bind_rows(matched_results)
cat("\n================ BẢNG 11B: KẾT QUẢ CRE TRÊN PHÂN MẪU GẦN VIỆT NAM (10, 12, 15 NƯỚC) ================\n")
print(subset(all_matched, inference == "cr2"), digits = 4, row.names = FALSE)
}
##
## [BƯỚC 11 - BỎ QUA]: Chưa có đủ 7 thông số tham chiếu của Việt Nam (iq_avg_m, pdi, idv, mas, uai, lto, ivr).
## -> Để chạy Bước 11, bạn chỉ cần thay các giá trị NA_real_ trong biến `vietnam_ref` ở đầu Bước 11
## bằng số liệu WGI & Hofstede của Việt Nam rồi bôi đen Bước 11 nhấn Ctrl + Enter!