#依不同研究群，看看效果量等的描述性統計
#程式報表11.2
dta |>
select(-Study) |>
gtsummary::tbl_summary(by=Group,
statistic = list(all_continuous() ~ "{mean} ({sd})"))
#依不同區域，看看效果量等的描述性統計
#程式報表11.3
dta |>
select(-Study) |>
gtsummary::tbl_summary(by=Region,
statistic = list(all_continuous() ~ "{mean} ({sd})"))
#效果量直方圖
#圖11.2
ggplot(data=dta,
aes(ES, after_stat(density)))+
geom_histogram(binwidth = function(x) 2 * IQR(x) / (length(x)^(1/3)),
col=1, fill=8)+
labs(x="研究所得效果（相關係數量）",
y="密度")
| Study | Region | Group | N | ES | yi | vi |
|---|---|---|---|---|---|---|
| 10 | 1 | 1 | 37 | 0.29 | 0.29 | 0.0233 |
| 20 | 1 | 1 | 83 | 0.33 | 0.33 | 0.0097 |
| 30 | 99 | 1 | 93 | 0.65 | 0.65 | 0.0036 |
| 41 | 1 | 1 | 46 | 0.62 | 0.62 | 0.0084 |
| 42 | 1 | 1 | 55 | 0.26 | 0.26 | 0.0161 |
| 43 | 1 | 1 | 29 | 0.44 | 0.44 | 0.0232 |
Random-Effects Model (k = 277; tau^2 estimator: REML)
logLik deviance AIC BIC AICc
49.8797 -99.7594 -95.7594 -88.5186 -95.7155
tau^2 (estimated amount of total heterogeneity): 0.0337 (SE = 0.0037)
tau (square root of estimated tau^2 value): 0.1836
I^2 (total heterogeneity / total variability): 90.46%
H^2 (total variability / sampling variability): 10.48
Test for Heterogeneity:
Q(df = 276) = 12088.7813, p-val < .0001
Model Results:
estimate se zval pval ci.lb ci.ub
0.3726 0.0125 29.6979 <.0001 0.3480 0.3972 
- 圖中的每一個線段，對應一個研究，其中黑點是效果量，線段則是95%信賴區間
-
圖中虛線標示效果量為0，如果信賴區間壓到這條虛線，代表效果可能為0，不顯著
Warning in cor.test.default(yi.star, vi, method = "kendall", exact = exact):
無法給連結計算精確 p 值
Rank Correlation Test for Funnel Plot Asymmetry
Kendall's tau = 0.3035, p < .0001
Regression Test for Funnel Plot Asymmetry
Model: mixed-effects meta-regression model
Predictor: standard error
Test for Funnel Plot Asymmetry: z = -7.3122, p < .0001
Limit Estimate (as sei -> 0): b = 0.6401 (CI: 0.5656, 0.7147)#後續要分析區域與研究群的調節效果
#先排除區域不明與第四組研究群體資料
dta <- dta |>
dplyr::filter(Region != 99 & Group < 4) |>
dplyr::mutate( Region = factor(Region),
Group = factor(Group))#抽取277筆中的13%，依區域與 ES 排序，繪製森林圖，標示區域
set.seed(201408)
dta_13 <- dta |>
dplyr::sample_frac(.13) |>
arrange(Region, ES)
#檢查區域效果。排除區域不明與第四組研究群體資料
#程式報表11.6
dta <-
dta |> arrange(Region, ES)
res_rgn <- rma(yi=yi, vi=vi, mods = ~ Region, data=dta)
res_rgn |> summary()
Mixed-Effects Model (k = 232; tau^2 estimator: REML)
logLik deviance AIC BIC AICc
41.3077 -82.6154 -70.6154 -50.0657 -70.2336
tau^2 (estimated amount of residual heterogeneity): 0.0337 (SE = 0.0041)
tau (square root of estimated tau^2 value): 0.1836
I^2 (residual heterogeneity / unaccounted variability): 88.12%
H^2 (unaccounted variability / sampling variability): 8.42
R^2 (amount of heterogeneity accounted for): 3.21%
Test for Residual Heterogeneity:
QE(df = 227) = 6206.1429, p-val < .0001
Test of Moderators (coefficients 2:5):
QM(df = 4) = 10.2446, p-val = 0.0365
Model Results:
estimate se zval pval ci.lb ci.ub
intrcpt 0.3990 0.0190 20.9591 <.0001 0.3617 0.4363
Region2 -0.0590 0.0351 -1.6810 0.0928 -0.1279 0.0098
Region3 -0.1334 0.0714 -1.8685 0.0617 -0.2734 0.0065
Region4 -0.1064 0.0398 -2.6742 0.0075 -0.1843 -0.0284
Region5 -0.0356 0.0556 -0.6414 0.5212 -0.1445 0.0733 #繪製 forest，標示區域（不只抽出一部份）
#圖11.6
p1 <- metaviz::viz_forest(res_rgn, group=dta[,'Region'],
xlab = "r", x_limit=c(-.80,1.1), type="study_only")
p2 <- metaviz::viz_forest(res_rgn, group=dta[,'Region'],
xlab = "r", x_limit=c(-.80,1.1),type="summary_only")
p1/p2+plot_layout(heights = c(7, 1))
#檢查出版偏誤，繪製 funnelplot，設定點的形狀
#圖11.7
metaviz::viz_funnel(res_rgn, group=dta[,"Region"],group_legend=T)+
scale_shape_manual(values = c(15:17,0,1)) Scale for shape is already present.
Adding another scale for shape, which will replace the existing scale.
#每組抽取10筆，依 Group 與 ES 排序，繪製森林圖，以 Group 標示
set.seed(2014)
dta1 <- dta[dta$Group==1,][sample(145,10),]
dta2 <- dta[dta$Group==2,][sample(36,10),]
dta3 <- dta[dta$Group==3,][sample(93,10),]
dta30 <- rbind(dta1,dta2,dta3)
dta30 <- dta30 |>
arrange(Group, ES)
#圖11.8
with(dta30, metafor::forest(yi, vi, slab = Group))