AI trong phân tích dữ liệu với R - Viện Phương Nam (6-11/1/2026)
Thach Tran
2025-12-30
Ngày 5: Meta-analysis
Việc 1. Đánh giá hiệu quả của can thiệp trên tình trạng cao huyết áp
1.1 Nhập liệu
study = c("ATMH", "HEP", "EWPHE", "HDFP", "MRC-1", "MRC-2",
"SHEP", "STOP", "Sy-Chi", "Sy-Eur")
n1 = c(780, 150, 90, 2427, 3546, 1314, 2365, 137, 1252, 2398)
mean1 = c(152, 190, 177, 152, 157, 182, 170, 195, 171, 174)
sd1 = c(15, 16, 16, 20, 16, 13, 10, 12, 11, 10)
n2 = c(750, 199, 82, 2370, 3445, 1337, 2371, 131, 1139, 2297)
mean2 = c(153, 191, 178, 151, 157, 182, 170, 194, 170, 174)
sd2 = c(16, 18, 15, 19, 16, 13, 9, 11, 11, 10)
df = data.frame(study, n1, mean1, sd1, n2, mean2, sd2)
head(df, 10)## study n1 mean1 sd1 n2 mean2 sd2
## 1 ATMH 780 152 15 750 153 16
## 2 HEP 150 190 16 199 191 18
## 3 EWPHE 90 177 16 82 178 15
## 4 HDFP 2427 152 20 2370 151 19
## 5 MRC-1 3546 157 16 3445 157 16
## 6 MRC-2 1314 182 13 1337 182 13
## 7 SHEP 2365 170 10 2371 170 9
## 8 STOP 137 195 12 131 194 11
## 9 Sy-Chi 1252 171 11 1139 170 11
## 10 Sy-Eur 2398 174 10 2297 174 101.2 Thực hiện phân tích gộp với fixed effects
library(metafor)## Warning: package 'metafor' was built under R version 4.3.2## Loading required package: Matrix## Loading required package: metadat## Warning: package 'metadat' was built under R version 4.3.2## Loading required package: numDeriv##
## Loading the 'metafor' package (version 4.4-0). For an
## introduction to the package please type: help(metafor)res_fe <- rma(measure = "MD", n1i = n1, n2i = n2, m1i = mean1, m2i = mean2, sd1i = sd1, sd2i = sd2, data = df, method = "FE")
summary(res_fe)##
## Fixed-Effects Model (k = 10)
##
## logLik deviance AIC BIC AICc
## -10.0963 9.7349 22.1927 22.4953 22.6927
##
## I^2 (total heterogeneity / total variability): 7.55%
## H^2 (total variability / sampling variability): 1.08
##
## Test for Heterogeneity:
## Q(df = 9) = 9.7349, p-val = 0.3724
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## 0.1411 0.1471 0.9594 0.3374 -0.1472 0.4294
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1- Note: measure = “SMD” (standardised mean difference); “SMDH” (standardised mean difference with heteroscedastic population variances)
1.3 Thực hiện phân tích gộp với random effects
res_re <- rma(measure = "MD", n1i = n1, n2i = n2, m1i = mean1, m2i = mean2, sd1i = sd1, sd2i = sd2, data = df, method = "REML")
summary(res_re)##
## Random-Effects Model (k = 10; tau^2 estimator: REML)
##
## logLik deviance AIC BIC AICc
## -9.9428 19.8855 23.8855 24.2800 25.8855
##
## tau^2 (estimated amount of total heterogeneity): 0 (SE = 0.0822)
## tau (square root of estimated tau^2 value): 0
## I^2 (total heterogeneity / total variability): 0.00%
## H^2 (total variability / sampling variability): 1.00
##
## Test for Heterogeneity:
## Q(df = 9) = 9.7349, p-val = 0.3724
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## 0.1411 0.1471 0.9594 0.3374 -0.1472 0.4294
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 11.4 Forest plot
forest(res_re, slab = df$study, xlab = "Mean difference in blood pressure (mmHg)", mlab = "Random-effects model")
1.5 Publication bias
Funnel plot
funnel(res_re, xlab = "Mean difference in blood pressure (mmHg)", ylab = "Standard error of MD")
Egger’s test
regtest(res_re, model = "lm") ##
## Regression Test for Funnel Plot Asymmetry
##
## Model: weighted regression with multiplicative dispersion
## Predictor: standard error
##
## Test for Funnel Plot Asymmetry: t = -0.0822, df = 8, p = 0.9365
## Limit Estimate (as sei -> 0): b = 0.1640 (CI: -0.5788, 0.9067)1.6 ChatGPT
PROMPT: I have a dataset for meta-analysis of continuous outcome. The dataset named “df” (see the uploaded file) with the following variables: study=name of study; n1=sample size in the treated group; mean1=mean blood pressure in the treated group; sd1=standard deviation in the treated group; n2=sample size in the control group; mean2=mean blood pressure in the control group; sd1=standard deviation in the control group. Please conduct a meta-analysis (using R) using the fixed-effects and random-effects model, forest plot, then evaluate publication bias.
Việc 2. Đánh giá hiệu quả của vaccine BCG
2.1 Đọc dữ liệu dat.bcg
data("dat.bcg")
head(dat.bcg, 13)## trial author year tpos tneg cpos cneg ablat alloc
## 1 1 Aronson 1948 4 119 11 128 44 random
## 2 2 Ferguson & Simes 1949 6 300 29 274 55 random
## 3 3 Rosenthal et al 1960 3 228 11 209 42 random
## 4 4 Hart & Sutherland 1977 62 13536 248 12619 52 random
## 5 5 Frimodt-Moller et al 1973 33 5036 47 5761 13 alternate
## 6 6 Stein & Aronson 1953 180 1361 372 1079 44 alternate
## 7 7 Vandiviere et al 1973 8 2537 10 619 19 random
## 8 8 TPT Madras 1980 505 87886 499 87892 13 random
## 9 9 Coetzee & Berjak 1968 29 7470 45 7232 27 random
## 10 10 Rosenthal et al 1961 17 1699 65 1600 42 systematic
## 11 11 Comstock et al 1974 186 50448 141 27197 18 systematic
## 12 12 Comstock & Webster 1969 5 2493 3 2338 33 systematic
## 13 13 Comstock et al 1976 27 16886 29 17825 33 systematic2.2 Thực hiện phân tích gộp với fixed effects
res_fe2 = rma(measure = "RR", ai = tpos, bi = tneg, ci = cpos, di = cneg, data = dat.bcg, method = "FE")
summary(res_fe2)##
## Fixed-Effects Model (k = 13)
##
## logLik deviance AIC BIC AICc
## -70.2236 152.2330 142.4471 143.0121 142.8108
##
## I^2 (total heterogeneity / total variability): 92.12%
## H^2 (total variability / sampling variability): 12.69
##
## Test for Heterogeneity:
## Q(df = 12) = 152.2330, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## -0.4303 0.0405 -10.6247 <.0001 -0.5097 -0.3509 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 12.3 Thực hiện phân tích gộp với random effects
res_re2 = rma(measure = "RR", ai = tpos, bi = tneg, ci = cpos, di = cneg, data = dat.bcg, method = "REML")
summary(res_re2)##
## Random-Effects Model (k = 13; tau^2 estimator: REML)
##
## logLik deviance AIC BIC AICc
## -12.2024 24.4047 28.4047 29.3746 29.7381
##
## tau^2 (estimated amount of total heterogeneity): 0.3132 (SE = 0.1664)
## tau (square root of estimated tau^2 value): 0.5597
## I^2 (total heterogeneity / total variability): 92.22%
## H^2 (total variability / sampling variability): 12.86
##
## Test for Heterogeneity:
## Q(df = 12) = 152.2330, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## -0.7145 0.1798 -3.9744 <.0001 -1.0669 -0.3622 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 12.4 Forest plot
forest(res_re2, slab = paste(dat.bcg$author, dat.bcg$year), xlab = "Relative Risk (RR)", atransf = exp,
at = log(c(0.2, 0.5, 1, 2, 5)), digits = 2)
1.5 Publication bias
Funnel plot
funnel(res_re2, xlab = "Log(RR)", ylab = "Standard error of Log(RR)")
Egger’s test
regtest(res_re2, model = "lm") ##
## Regression Test for Funnel Plot Asymmetry
##
## Model: weighted regression with multiplicative dispersion
## Predictor: standard error
##
## Test for Funnel Plot Asymmetry: t = -1.4013, df = 11, p = 0.1887
## Limit Estimate (as sei -> 0): b = -0.1909 (CI: -0.6753, 0.2935)2.5 ChatGPT
PROMPT: I perform a meta-analysis of a binary outcome using the dataset ‘dat.bcg’ in the ‘metafor’ R package. The dataset ‘dat.bcg’ includes: author = paper authors; year = publication year; tpos = positive (to tuberculosis) cases in the vaccinated group; tneg = negative cases in the vaccinated group; cpos = positive cases in the control group; cneg = negative cases in the control group. Please conduct a meta-analysis (using R) for the fixed-, random-effects models, forest plot (with the x-axis presenting relative risk and corresponding 95% CI) and evaluate publication bias.