整合分析(Meta-Analysis)

許清芳,鄭中平

2024-04-27

小題

(1)請用下列方式讀入資料,並繪製森林圖。 dta <- read.csv(“https://raw.githubusercontent.com/ChungPingCheng/R4BS2/main/Data/Deindividuation.csv”, header=T)

(2)請計算平均效果量,並檢驗是否為0。

(3)請繪製漏斗圖,並檢驗是否有出版偏誤。

(4)請檢驗匿名團體的參與者是否在群體中,是否會影響研究中效果的大小。

提示

先找到範例中對應的程式複製到新的 RMD,再往下修改

看一下資料

#整體設定,含載入套件
source("https://raw.githubusercontent.com/ChungPingCheng/R4BS2/main/R4BS_setup.R")
#讀進資料
dta <- read.csv("https://raw.githubusercontent.com/ChungPingCheng/R4BS2/main/Data/Deindividuation.csv", header=T) 
#顯示前六筆資料
dim(dta)
[1] 70  5
head(dta)
Study r N Par_population Par_in_groups
Becker-Haven & Lindskold, 1978 0.07 68 1 1
Becker-Haven & Lindskold, 1979 0.60 68 1 1
Carver, 1974 -0.38 32 1 2
Carver, 1975, Study 1 -0.11 40 1 2
Carver, 1975, Study 2 -0.12 44 1 2
Diener, 1976 -0.37 60 1 1

(1)請用下列方式讀入資料,並繪製森林圖。

描述性統計

#依不同Par_in_groups,看看效果量等的描述性統計
dta |>
  select(-Study) |>
  gtsummary::tbl_summary(by=Par_in_groups,
                statistic = list(all_continuous() ~ "{mean} ({sd})"))

資料圖形化

#效果量直方圖
ggplot(data=dta, 
       aes(r, after_stat(density)))+
  geom_histogram(binwidth = function(x) 2 * IQR(x) / (length(x)^(1/3)),
                 col=1, fill=8)+
  labs(x="研究所得效果(相關係數量)",
       y="密度")

(2)請計算平均效果量,並檢驗是否為0。

整合分析

#利用相關與樣本數計算相關的抽樣變異,並放進原先資料
dta <- metafor::escalc(measure="COR", ri=r, ni=N, data=dta)
head(dta)
Study r N Par_population Par_in_groups yi vi
Becker-Haven & Lindskold, 1978 0.07 68 1 1 0.07 0.0148
Becker-Haven & Lindskold, 1979 0.60 68 1 1 0.60 0.0061
Carver, 1974 -0.38 32 1 2 -0.38 0.0236
Carver, 1975, Study 1 -0.11 40 1 2 -0.11 0.0250
Carver, 1975, Study 2 -0.12 44 1 2 -0.12 0.0226
Diener, 1976 -0.37 60 1 1 -0.37 0.0126
#計算所有研究的平均效果
summary(res <- metafor::rma(yi=yi, vi=vi, data=dta))

Random-Effects Model (k = 70; tau^2 estimator: REML)

  logLik  deviance       AIC       BIC      AICc   
 -7.2204   14.4408   18.4408   22.9090   18.6226   

tau^2 (estimated amount of total heterogeneity): 0.0552 (SE = 0.0121)
tau (square root of estimated tau^2 value):      0.2350
I^2 (total heterogeneity / total variability):   82.94%
H^2 (total variability / sampling variability):  5.86

Test for Heterogeneity:
Q(df = 69) = 394.9507, p-val < .0001

Model Results:

estimate      se    zval    pval   ci.lb   ci.ub 
  0.0985  0.0320  3.0777  0.0021  0.0358  0.1613 

第一小題

森林圖,圖示結果

metaviz::viz_forest(res)

(3)請繪製漏斗圖,並檢驗是否有出版偏誤。

檢查出版偏誤

漏斗圖

#檢查出版偏誤,先繪製 funnelplot
metaviz::viz_funnel(res)

ranktest(res)
Warning in cor.test.default(yi.star, vi, method = "kendall", exact = exact):
無法給連結計算精確 p 值

Rank Correlation Test for Funnel Plot Asymmetry

Kendall's tau = -0.0593, p = 0.4684
regtest(res)

Regression Test for Funnel Plot Asymmetry

Model:     mixed-effects meta-regression model
Predictor: standard error

Test for Funnel Plot Asymmetry: z = -0.6800, p = 0.4965
Limit Estimate (as sei -> 0):   b =  0.1718 (CI: -0.0485, 0.3920)

(4)請檢驗匿名團體的參與者是否在群體中,是否會影響研究中效果的大小。

整合調節分析

dta <- dta |>
  dplyr::mutate(Par_in_groups = factor(Par_in_groups))
set.seed(201408)
dta_13 <- dta |> 
  dplyr::sample_frac(.33) |> 
  arrange(Par_in_groups, r)
#圖11.5
with(dta_13, metafor::forest(yi, vi, slab = Par_in_groups))

dta <- 
  dta |> arrange(Par_in_groups, r)
res_rgn <- rma(yi=yi, vi=vi, mods = ~ Par_in_groups, data=dta)
res_rgn |> summary()

Mixed-Effects Model (k = 70; tau^2 estimator: REML)

  logLik  deviance       AIC       BIC      AICc   
 -6.9130   13.8261   19.8261   26.4846   20.2011   

tau^2 (estimated amount of residual heterogeneity):     0.0547 (SE = 0.0121)
tau (square root of estimated tau^2 value):             0.2339
I^2 (residual heterogeneity / unaccounted variability): 82.75%
H^2 (unaccounted variability / sampling variability):   5.80
R^2 (amount of heterogeneity accounted for):            0.95%

Test for Residual Heterogeneity:
QE(df = 68) = 381.5628, p-val < .0001

Test of Moderators (coefficient 2):
QM(df = 1) = 1.4037, p-val = 0.2361

Model Results:

                estimate      se     zval    pval    ci.lb   ci.ub 
intrcpt           0.1267  0.0398   3.1854  0.0014   0.0488  0.2047 
Par_in_groups2   -0.0789  0.0666  -1.1848  0.2361  -0.2094  0.0516 
#繪製 forest
p1 <- metaviz::viz_forest(res_rgn, group=dta[,'Par_in_groups'],
              xlab = "r", x_limit=c(-.80,1.1), type="study_only")
p2 <- metaviz::viz_forest(res_rgn, group=dta[,'Par_in_groups'],
              xlab = "r", x_limit=c(-.80,1.1),type="summary_only")
p1/p2+plot_layout(heights = c(7, 1))