############################################################
# 0. CÀI / NẠP GÓI
############################################################
# Chỉ cài nếu chưa có:
# install.packages(c("readxl","janitor","dplyr","stringr","metafor","purrr"))
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)
############################################################
# 1. ĐỌC FILE CODEBOOK + CHUẨN HOÁ TÊN CỘT
############################################################
# Sửa lại đường dẫn này cho đúng với máy của anh/chị:
file_path <- "/Volumes/PhD - NTTU/CHUYÊN ĐỀ 1/Meta-analysis/10.10/FILE TRICH XUAT DUONG DAN 12.11/FILE TONG/Codebook_bivariate MA_Confirmatory 2.xlsx"
codebook_raw <- read_excel(
path = file_path,
sheet = "Codebook"
)
# Chuẩn hoá tên cột: sang snake_case, bỏ dấu, bỏ khoảng trắng
codebook <- codebook_raw %>%
clean_names()
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" "bang_hinh_trang"
# 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" "r" "fisher_z" "se_z" "nguon_r"
# "bang_hinh_trang" "ghi_chu"
# Trong các version trước dùng "canh_h_number", nên đổi tên nếu cần:
if (!"canh_h_number" %in% names(codebook) && "canh_h" %in% names(codebook)) {
codebook <- codebook %>%
rename(canh_h_number = canh_h)
}
############################################################
# 2. TẠO CỘT canh_h (H1…H12) + CHUẨN HOÁ CÁC CỘT SỐ
############################################################
codebook2 <- codebook %>%
mutate(
# Chuẩn hoá cột Cạnh_H# sau khi clean_names thành 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 = "H2","H3",... từ mẫu "H" + số
canh_h = str_to_upper(str_extract(canh_h_number_chr, "H\\d+")),
# Chuẩn hoá fisher_z, se_z, n_tong, n_cap về dạng số
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)))
)
unique(codebook2$canh_h)
## [1] "H12" "H4" "H11" "H6" "H9" "H5" "H3" "H10" "H8" "H7" "H2"
summary(codebook2$fisher_z)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.1758 0.5347 0.6780 0.7211 0.8830 1.5700
############################################################
# 3. TẠO r_num TỪ fisher_z (ƯU TIÊN) HOẶC TỪ r GỐC
############################################################
codebook3 <- codebook2 %>%
mutate(
# Parse r gốc thành số (nếu có dùng dấu phẩy)
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
)
)
summary(codebook3$r_num)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.1740 0.4890 0.5902 0.5936 0.7079 0.9170
############################################################
# 4. LỌC DỮ LIỆU CHO MA + HOÀN THIỆN fisher_z, SE_z
############################################################
dat_ma <- codebook3 %>%
# 1) Chỉ giữ các cạnh H1–H12
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_cap, nếu NA thì dùng n_tong
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 đang trống
mutate(
fisher_z = if_else(
is.na(fisher_z),
0.5 * log((1 + r_num) / (1 - r_num)),
fisher_z
),
# 5) Tính SE_z nếu trống
se_z = if_else(
is.na(se_z),
1 / sqrt(n_effective - 3),
se_z
)
) %>%
# 6) Đảm bảo n_effective ≥ 4
filter(n_effective >= 4) %>%
# 7) Đặt lại r = r_num cho gọn
mutate(r = r_num)
dat_ma %>% count(canh_h)
## # A tibble: 11 × 2
## canh_h n
## <chr> <int>
## 1 H10 4
## 2 H11 17
## 3 H12 26
## 4 H2 1
## 5 H3 7
## 6 H4 13
## 7 H5 5
## 8 H6 51
## 9 H7 1
## 10 H8 2
## 11 H9 14
summary(dat_ma$r)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.1740 0.4890 0.5902 0.5936 0.7079 0.9170
############################################################
# 5. BẢNG TỔNG QUAN THEO GIẢ THUYẾT (BẢNG 8)
############################################################
bang_1_1 <- dat_ma %>%
group_by(canh_h) %>%
summarise(
k = n(),
m = n_distinct(study_id),
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, # đủ điều kiện MA song biến
.groups = "drop"
) %>%
arrange(canh_h)
bang_1_1
## # A tibble: 11 × 7
## canh_h k m r_min r_med r_max du_lieu_MA
## <chr> <int> <int> <dbl> <dbl> <dbl> <lgl>
## 1 H10 4 3 0.451 0.572 0.741 TRUE
## 2 H11 17 9 0.174 0.639 0.829 TRUE
## 3 H12 26 20 0.400 0.685 0.917 TRUE
## 4 H2 1 1 0.467 0.467 0.467 FALSE
## 5 H3 7 4 0.184 0.633 0.721 TRUE
## 6 H4 13 5 0.329 0.476 0.773 TRUE
## 7 H5 5 2 0.446 0.477 0.634 TRUE
## 8 H6 51 7 0.407 0.602 0.835 TRUE
## 9 H7 1 1 0.705 0.705 0.705 FALSE
## 10 H8 2 2 0.403 0.568 0.733 TRUE
## 11 H9 14 10 0.268 0.512 0.813 TRUE
# Nếu muốn thêm dòng H1 = 0,0,… cho giống bảng trong luận án:
bang_1_1_with_H1 <- bang_1_1 %>%
bind_rows(
tibble(
canh_h = "H1",
k = 0L,
m = 0L,
r_min = NA_real_,
r_med = NA_real_,
r_max = NA_real_,
du_lieu_MA = FALSE
)
) %>%
arrange(canh_h)
bang_1_1_with_H1
## # A tibble: 12 × 7
## canh_h k m r_min r_med r_max du_lieu_MA
## <chr> <int> <int> <dbl> <dbl> <dbl> <lgl>
## 1 H1 0 0 NA NA NA FALSE
## 2 H10 4 3 0.451 0.572 0.741 TRUE
## 3 H11 17 9 0.174 0.639 0.829 TRUE
## 4 H12 26 20 0.400 0.685 0.917 TRUE
## 5 H2 1 1 0.467 0.467 0.467 FALSE
## 6 H3 7 4 0.184 0.633 0.721 TRUE
## 7 H4 13 5 0.329 0.476 0.773 TRUE
## 8 H5 5 2 0.446 0.477 0.634 TRUE
## 9 H6 51 7 0.407 0.602 0.835 TRUE
## 10 H7 1 1 0.705 0.705 0.705 FALSE
## 11 H8 2 2 0.403 0.568 0.733 TRUE
## 12 H9 14 10 0.268 0.512 0.813 TRUE
############################################################
# 6. CHỌN H# ĐỦ DỮ LIỆU ĐỂ CHẠY MA
############################################################
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"
# Kỳ vọng: H2–H12 (H1 không đủ dữ liệu nên không có trong danh sách)
############################################################
# 7. HÀM CHẠY META-ANALYSIS CHO 1 GIẢ THUYẾT
############################################################
run_ma_for_h <- function(df_h) {
n_effects <- nrow(df_h)
n_studies <- 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,
m = n_studies,
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)
)
}
############################################################
# 8. BẢNG KẾT QUẢ MA TỔNG HỢP (TƯƠNG ĐƯƠNG results_all)
############################################################
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
## # A tibble: 9 × 15
## canh_h k m r_bar ci95_lower ci95_upper tau2 I2 Q Q_df
## <chr> <int> <int> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 H10 4 3 0.599 0.451 0.715 0.042 95.3 69.9 3
## 2 H11 17 9 0.618 0.528 0.695 0.078 96.6 463. 16
## 3 H12 26 20 0.71 0.639 0.769 0.113 97.9 1155. 25
## 4 H3 7 4 0.594 0.462 0.701 0.059 94.6 124. 6
## 5 H4 13 5 0.554 0.444 0.648 0.072 97.2 406. 12
## 6 H5 5 2 0.502 0.448 0.553 0.003 47.2 8.2 4
## 7 H6 51 7 0.612 0.578 0.643 0.032 91.1 571. 50
## 8 H8 2 2 0.593 0.182 0.827 0.127 98 50 1
## 9 H9 14 10 0.554 0.46 0.636 0.054 93.7 202. 13
## # ℹ 5 more variables: Q_p <dbl>, pi95_lower <dbl>, pi95_upper <dbl>,
## # r_min <dbl>, r_max <dbl>
############################################################
# 9. HÀM XEM CHI TIẾT CHO 1 GIẢ THUYẾT (phục vụ viết luận)
############################################################
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 MA."))
return(NULL)
}
res <- rma(yi = fisher_z, sei = se_z,
method = "REML", data = df_h)
pred <- predict(res, transf = transf.ztor)
list(
model = res,
pred = pred,
r_range = range(df_h$r, na.rm = TRUE),
k = nrow(df_h),
m = n_distinct(df_h$study_id)
)
}
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ụ:
# print_ma_pretty("H3", "Thuộc tính điểm đến (DA) → Tính xác thực (AU)")
# print_ma_pretty("H10", "Trải nghiệm đáng nhớ (ME) → Hài lòng (SAT)")