############################################################
## Thay đường dẫn sau bằng đường dẫn có chứa file code book trong máy: "/Volumes/PhD - NTTU/CHUYÊN ĐỀ 1/Meta-analysis/10.10/FILE TRICH XUAT DUONG DAN 12.11/FILE TONG/Codebook_bivariate MA_Confirmatory.xlsx"
# 0. CÀI ĐẶT GÓI (chỉ cần chạy 1 lần nếu chưa cài)
############################################################
# install.packages(c("readxl", "janitor", "dplyr", "stringr", "metafor", "purrr"))
############################################################
# 1. NẠP THƯ VIỆN
############################################################
library(readxl)
library(janitor)
##
## Attaching package: 'janitor'
## The following objects are masked from 'package:stats':
##
## chisq.test, fisher.test
library(dplyr)
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(stringr)
library(metafor)
## Loading required package: Matrix
## Loading required package: metadat
## Loading required package: numDeriv
##
## Loading the 'metafor' package (version 4.8-0). For an
## introduction to the package please type: help(metafor)
library(purrr)
############################################################
# 2. ĐỌC FILE EXCEL + CHUẨN HOÁ TÊN CỘT
############################################################
# Đọc file codebook bivariate MA
codebook_raw <- read_excel(
"/Volumes/PhD - NTTU/CHUYÊN ĐỀ 1/Meta-analysis/10.10/FILE TRICH XUAT DUONG DAN 12.11/FILE TONG/Codebook_bivariate MA_Confirmatory.xlsx",
sheet = "Codebook" # nếu tên sheet khác, sửa lại chỗ này
)
# Chuẩn hoá tên cột: sang chữ thường, bỏ dấu, bỏ khoảng trắng
codebook <- codebook_raw %>%
clean_names()
# Kiểm tra tên cột (tham chiếu)
names(codebook)
## [1] "study_id" "nam" "boi_canh_khu_vuc" "n_tong"
## [5] "n_cap" "x_bien_goc" "y_bien_goc" "x_nhom_sor"
## [9] "y_nhom_sor" "canh_h_number" "r" "fisher_z"
## [13] "se_z" "nguon_r"
# Kỳ vọng: "study_id" "nam" "boi_canh_khu_vuc" "n_tong" "n_cap"
# "x_bien_goc" "y_bien_goc" "x_nhom_sor" "y_nhom_sor"
# "canh_h_number" "r" "fisher_z" "se_z" "nguon_r"
############################################################
# 3. TẠO CỘT canh_h (H1…H12) + CHUẨN HOÁ CÁC CỘT SỐ
############################################################
# Mục tiêu: tạo biến canh_h = "H1","H2",..., từ cột canh_h_number
codebook2 <- codebook %>%
mutate(
# Chuẩn hoá cột canh_h_number về chuỗi gọn
canh_h_number_chr = as.character(canh_h_number),
canh_h_number_chr = trimws(canh_h_number_chr),
# Tạo cột canh_h bằng cách trích mẫu "H" + số (VD: "H3: DA→AU" -> "H3")
canh_h = str_to_upper(str_extract(canh_h_number_chr, "H\\d+")),
# Chuẩn hoá fisher_z, se_z, N về số (nếu có dấu phẩy thập phân thì đổi sang chấm)
fisher_z = as.numeric(gsub(",", ".", as.character(fisher_z))),
se_z = as.numeric(gsub(",", ".", as.character(se_z))),
n_tong = as.numeric(gsub(",", ".", as.character(n_tong))),
n_cap = as.numeric(gsub(",", ".", as.character(n_cap)))
)
# Kiểm tra nhanh các giá trị H#
unique(codebook2$canh_h)
## [1] "H12" "H4" "H11" "H6" "H5" "H3" "H9" "H10" "H8" "H7"
summary(codebook2$fisher_z)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.1758 0.5348 0.6575 0.7046 0.8480 1.5220
############################################################
# 4. TẠO r_num TỪ fisher_z (ƯU TIÊN) HOẶC TỪ CỘT r GỐC
############################################################
# Nếu đã có fisher_z: r_num = tanh(fisher_z)
# Nếu fisher_z bị NA: cố gắng parse r gốc.
codebook3 <- codebook2 %>%
mutate(
# Thử parse r gốc thành số (dùng dấu chấm thay dấu phẩy nếu cần)
r_parsed = as.numeric(gsub(",", ".", as.character(r))),
# Nếu có fisher_z thì suy r từ fisher_z; nếu không, dùng r_parsed
r_num = case_when(
!is.na(fisher_z) ~ tanh(fisher_z),
is.na(fisher_z) ~ r_parsed
)
)
# Xem thử phân bố r_num
summary(codebook3$r_num)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.1740 0.4890 0.5767 0.5867 0.6900 0.9090
############################################################
# 5. LỌC DỮ LIỆU MA + HOÀN THIỆN fisher_z, SE_z
############################################################
# Mục tiêu: tạo dat_ma – tập dữ liệu sạch cho MA bivariate
dat_ma <- codebook3 %>%
# 1. Chỉ giữ các cạnh H1–H12 (bỏ H13 – moderator)
filter(!is.na(canh_h),
canh_h %in% paste0("H", 1:12)) %>%
# 2. Loại case thiếu r_num hoặc r_num không hợp lệ
filter(!is.na(r_num),
r_num > -1, r_num < 1) %>%
# 3. Xác định N_effective: ưu tiên N_cặp, nếu NA thì dùng N_tổng
mutate(
n_effective = if_else(!is.na(n_cap), n_cap, n_tong),
n_effective = if_else(is.na(n_effective), n_tong, n_effective)
) %>%
# 4. Tính Fisher’s z nếu còn trống, dựa trên r_num
mutate(
fisher_z = if_else(
is.na(fisher_z),
0.5 * log((1 + r_num) / (1 - r_num)), # z = 1/2 * ln((1+r)/(1-r))
fisher_z
),
# 5. Tính SE_z nếu còn trống
se_z = if_else(
is.na(se_z),
1 / sqrt(n_effective - 3), # SE_z = 1 / sqrt(N - 3)
se_z
)
) %>%
# 6. Đảm bảo N_effective ≥ 4
filter(n_effective >= 4) %>%
# 7. Đặt lại tên r cho gọn (dùng r_num làm r)
mutate(r = r_num)
# Kiểm tra sau khi lọc
dat_ma %>% count(canh_h)
summary(dat_ma$r)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.1740 0.4890 0.5767 0.5867 0.6900 0.9090
############################################################
# 6. BẢNG 1.1 – TỔNG QUAN DỮ LIỆU THEO GIẢ THUYẾT
############################################################
# k = số hiệu ứng (số hệ số r)
# m = số nghiên cứu (số Study_ID khác nhau)
bang_1_1 <- dat_ma %>%
group_by(canh_h) %>%
summarise(
k = n(), # số hiệu ứng
m = n_distinct(study_id), # số nghiên cứu
r_min = min(r, na.rm = TRUE),
r_med = median(r, na.rm = TRUE),
r_max = max(r, na.rm = TRUE),
du_lieu_MA = k >= 2, # k ≥ 2: đủ điều kiện MA bivariate
du_lieu_TSSEM1 = m >= 3, # m ≥ 3: đủ “đẹp” cho TSSEM-1 (tuỳ bạn chỉnh)
.groups = "drop"
) %>%
arrange(canh_h)
bang_1_1
############################################################
# 7. XÁC ĐỊNH H# ĐỦ ĐIỀU KIỆN MA (k ≥ 2 nghiên cứu/hiệu ứng)
############################################################
# Ở đây đơn giản dùng k (số hiệu ứng) ≥ 2
h_list_ma <- bang_1_1 %>%
filter(du_lieu_MA) %>%
pull(canh_h)
h_list_ma
## [1] "H10" "H11" "H12" "H3" "H4" "H5" "H6" "H8" "H9"
############################################################
# 8. HÀM CHẠY META-ANALYSIS CHO MỘT H#
############################################################
# Lưu ý:
# - k = n_effects (số hiệu ứng)
# - m = n_studies (số nghiên cứu)
run_ma_for_h <- function(df_h) {
n_effects <- nrow(df_h)
n_studies <- dplyr::n_distinct(df_h$study_id)
res <- rma(yi = fisher_z, sei = se_z, method = "REML", data = df_h)
pred <- predict(res, transf = transf.ztor)
tibble(
k = n_effects, # số hiệu ứng
m = n_studies, # số nghiên cứu
r_bar = round(pred$pred, 3),
ci95_lower = round(pred$ci.lb, 3),
ci95_upper = round(pred$ci.ub, 3),
tau2 = round(res$tau2, 3),
I2 = round(res$I2, 1),
Q = round(res$QE, 1),
Q_df = res$k - 1,
Q_p = signif(res$QEp, 3),
pi95_lower = round(pred$cr.lb, 3),
pi95_upper = round(pred$cr.ub, 3),
r_min = round(min(df_h$r, na.rm = TRUE), 3),
r_max = round(max(df_h$r, na.rm = TRUE), 3)
)
}
############################################################
# 9. BẢNG TỔNG HỢP KẾT QUẢ MA CHO CÁC H# ĐỦ ĐIỀU KIỆN
# (TƯƠNG ĐƯƠNG BẢNG 2.T – EDGE DECISION)
############################################################
if (length(h_list_ma) == 0) {
warning("Không có giả thuyết nào có k ≥ 2 sau khi lọc dữ liệu.")
results_all <- tibble()
} else {
results_all <- dat_ma %>%
filter(canh_h %in% h_list_ma) %>%
group_by(canh_h) %>%
group_modify(~ run_ma_for_h(.x)) %>%
ungroup() %>%
arrange(canh_h)
}
results_all
############################################################
# 10. HÀM LẤY CHI TIẾT CHO 1 GIẢ THUYẾT H#
############################################################
get_ma_for_one_H <- function(h_code) {
df_h <- dat_ma %>% filter(canh_h == h_code)
if (nrow(df_h) < 2) {
warning(paste("Giả thuyết", h_code, "có k <", 2, "- không đủ dữ liệu để chạy meta-analysis."))
return(NULL)
}
res <- rma(yi = fisher_z, sei = se_z, method = "REML", data = df_h)
pred <- predict(res, transf = transf.ztor)
list(
model = res, # đối tượng rma()
pred = pred, # r̄, CI, PI trên thang r
r_range = range(df_h$r, na.rm = TRUE), # r_min, r_max quan sát
k = nrow(df_h),
m = dplyr::n_distinct(df_h$study_id)
)
}
############################################################
# 11. HÀM IN KẾT QUẢ "ĐẸP" CHO TỪNG H# (PHỤC VỤ BẢNG 2.2–2.9)
############################################################
print_ma_pretty <- function(h_code, edge_label = "") {
obj <- get_ma_for_one_H(h_code)
if (is.null(obj)) return(invisible(NULL))
res <- obj$model
pred <- obj$pred
rr <- obj$r_range
k <- obj$k
m <- obj$m
rbar <- as.numeric(pred$pred)
ci_l <- as.numeric(pred$ci.lb)
ci_u <- as.numeric(pred$ci.ub)
tau2 <- res$tau2
tau <- sqrt(res$tau2)
I2 <- res$I2
Q <- res$QE
dfQ <- res$k - 1
pQ <- res$QEp
pi_l <- as.numeric(pred$cr.lb)
pi_u <- as.numeric(pred$cr.ub)
cat("Bảng: Kết quả MA cho", h_code,
if (edge_label != "") paste0(" – ", edge_label), "\n", sep = "")
cat(sprintf("%s: k = %d hiệu ứng, m = %d nghiên cứu, r̄ = %.3f (CI 95%%: %.3f – %.3f)\n",
h_code, k, m, rbar, ci_l, ci_u))
cat(sprintf("τ² = %.3f (τ = %.3f)\tI² = %.1f%% (Q(%d) = %.1f, p = %.3g)\n",
tau2, tau, I2, dfQ, Q, pQ))
cat(sprintf("PI 95%%: %.3f – %.3f\t(r quan sát dao động từ %.3f – %.3f)\n\n",
pi_l, pi_u, rr[1], rr[2]))
}
# Ví dụ sử dụng:
# print_ma_pretty("H3", "Thuộc tính điểm đến (DA) → Tính xác thực (AU)")
# print_ma_pretty("H4", "Thuộc tính điểm đến (DA) → Trải nghiệm đáng nhớ (ME)")
# print_ma_pretty("H5", "Trải nghiệm 4E → Tính xác thực (AU)")
# print_ma_pretty("H6", "Trải nghiệm 4E → Trải nghiệm đáng nhớ (ME)")
# print_ma_pretty("H8", "Tính xác thực (AU) → Hài lòng (SAT)")
# print_ma_pretty("H9", "Tính xác thực (AU) → Trung thành (LI)")
# print_ma_pretty("H10", "Trải nghiệm đáng nhớ (ME) → Hài lòng (SAT)")
# print_ma_pretty("H11", "Trải nghiệm đáng nhớ (ME) → Trung thành (LI)")
# print_ma_pretty("H12", "Hài lòng (SAT) → Trung thành (LI)")
print_many_ma <- function(h_codes, labels = NULL) {
if (is.null(labels)) {
labels <- rep("", length(h_codes))
}
if (length(labels) != length(h_codes)) {
stop("Số lượng nhãn (labels) phải bằng số lượng H# (h_codes).")
}
for (i in seq_along(h_codes)) {
print_ma_pretty(h_codes[i], labels[i])
}
}
h_vec <- c("H3","H4","H5","H6","H8","H9","H10","H11","H12")
label_vec <- c(
"Thuộc tính điểm đến (DA) → Tính xác thực (AU)",
"Thuộc tính điểm đến (DA) → Trải nghiệm đáng nhớ (ME)",
"Trải nghiệm 4E → Tính xác thực (AU)",
"Trải nghiệm 4E → Trải nghiệm đáng nhớ (ME)",
"Tính xác thực (AU) → Hài lòng (SAT)",
"Tính xác thực (AU) → Trung thành (LI)",
"Trải nghiệm đáng nhớ (ME) → Hài lòng (SAT)",
"Trải nghiệm đáng nhớ (ME) → Trung thành (LI)",
"Hài lòng (SAT) → Trung thành (LI)"
)
print_many_ma(h_vec, label_vec)
## Bảng: Kết quả MA choH3 – Thuộc tính điểm đến (DA) → Tính xác thực (AU)
## H3: k = 7 hiệu ứng, m = 4 nghiên cứu, r̄ = 0.594 (CI 95%: 0.462 – 0.701)
## τ² = 0.059 (τ = 0.242) I² = 94.6% (Q(6) = 124.2, p = 2.11e-24)
## PI 95%: 0.173 – 0.832 (r quan sát dao động từ 0.184 – 0.721)
##
## Bảng: Kết quả MA choH4 – Thuộc tính điểm đến (DA) → Trải nghiệm đáng nhớ (ME)
## H4: k = 6 hiệu ứng, m = 4 nghiên cứu, r̄ = 0.619 (CI 95%: 0.477 – 0.729)
## τ² = 0.063 (τ = 0.251) I² = 96.8% (Q(5) = 138.3, p = 4.15e-28)
## PI 95%: 0.189 – 0.850 (r quan sát dao động từ 0.339 – 0.773)
##
## Bảng: Kết quả MA choH5 – Trải nghiệm 4E → Tính xác thực (AU)
## H5: k = 8 hiệu ứng, m = 3 nghiên cứu, r̄ = 0.500 (CI 95%: 0.470 – 0.529)
## τ² = 0.000 (τ = 0.001) I² = 0.0% (Q(7) = 8.3, p = 0.304)
## PI 95%: 0.470 – 0.529 (r quan sát dao động từ 0.446 – 0.634)
##
## Bảng: Kết quả MA choH6 – Trải nghiệm 4E → Trải nghiệm đáng nhớ (ME)
## H6: k = 50 hiệu ứng, m = 7 nghiên cứu, r̄ = 0.596 (CI 95%: 0.561 – 0.630)
## τ² = 0.033 (τ = 0.182) I² = 91.1% (Q(49) = 564.9, p = 7.46e-89)
## PI 95%: 0.315 – 0.781 (r quan sát dao động từ 0.407 – 0.835)
##
## Bảng: Kết quả MA choH8 – Tính xác thực (AU) → Hài lòng (SAT)
## H8: k = 2 hiệu ứng, m = 2 nghiên cứu, r̄ = 0.593 (CI 95%: 0.182 – 0.827)
## τ² = 0.127 (τ = 0.356) I² = 98.0% (Q(1) = 50.0, p = 1.56e-12)
## PI 95%: -0.173 – 0.912 (r quan sát dao động từ 0.403 – 0.733)
##
## Bảng: Kết quả MA choH9 – Tính xác thực (AU) → Trung thành (LI)
## H9: k = 11 hiệu ứng, m = 8 nghiên cứu, r̄ = 0.582 (CI 95%: 0.476 – 0.670)
## τ² = 0.056 (τ = 0.237) I² = 93.3% (Q(10) = 141.8, p = 1.77e-25)
## PI 95%: 0.176 – 0.819 (r quan sát dao động từ 0.337 – 0.813)
##
## Bảng: Kết quả MA choH10 – Trải nghiệm đáng nhớ (ME) → Hài lòng (SAT)
## H10: k = 4 hiệu ứng, m = 4 nghiên cứu, r̄ = 0.587 (CI 95%: 0.422 – 0.715)
## τ² = 0.049 (τ = 0.222) I² = 95.4% (Q(3) = 69.8, p = 4.77e-15)
## PI 95%: 0.182 – 0.822 (r quan sát dao động từ 0.451 – 0.741)
##
## Bảng: Kết quả MA choH11 – Trải nghiệm đáng nhớ (ME) → Trung thành (LI)
## H11: k = 11 hiệu ứng, m = 7 nghiên cứu, r̄ = 0.584 (CI 95%: 0.454 – 0.691)
## τ² = 0.090 (τ = 0.299) I² = 97.0% (Q(10) = 351.7, p = 1.75e-69)
## PI 95%: 0.055 – 0.857 (r quan sát dao động từ 0.174 – 0.829)
##
## Bảng: Kết quả MA choH12 – Hài lòng (SAT) → Trung thành (LI)
## H12: k = 22 hiệu ứng, m = 18 nghiên cứu, r̄ = 0.691 (CI 95%: 0.620 – 0.751)
## τ² = 0.087 (τ = 0.295) I² = 97.3% (Q(21) = 894.6, p = 7.33e-176)
## PI 95%: 0.253 – 0.894 (r quan sát dao động từ 0.400 – 0.909)
bang_1_1
results_all