R meta analysis (#2 practice)
Sangyoung
2023-07-17
library(tidyverse)
library(metafor)
library(readxl)
library(countrycode)
library(broom)
library(scales)
theme_set(theme_light())
rm(list = ls()) ## to clear global enviro if neededData loading
setwd("d:/r/training/")
dat <-
read_excel("Practice Wellspring.xlsx") %>%
mutate(region = countrycode(country, origin = 'country.name', destination = c('region'))) %>%
mutate(authors = str_remove_all(authors, "[:digit:]")) %>%
mutate(authors = glue::glue("{authors} ({year_pub})")) %>%
mutate(row_id = row_number())ROADS INTERVENTION: OUTCOME = Volume
Exploratory data analysis
dat %>%
mutate(estimate_id = str_remove(id, "_.*")) %>%
count(estimate_id) ## # A tibble: 18 × 2
## estimate_id n
## <chr> <int>
## 1 73628845 1
## 2 73628851 1
## 3 73628853 1
## 4 73628861 1
## 5 73628867 1
## 6 73797268 1
## 7 73806573 1
## 8 74666345 1
## 9 74713077 1
## 10 74716985 1
## 11 74719673 1
## 12 75350923 1
## 13 75413202 1
## 14 75413210 1
## 15 75413213 1
## 16 75431749 1
## 17 77212051 1
## 18 77270650 1- Only one outcome per study
dat %>%
count(inter_subcode) %>%
mutate(pct = n / 18) %>%
mutate(inter_subcode = fct_reorder(inter_subcode, pct)) %>%
ggplot(aes(pct, inter_subcode, fill = inter_subcode)) +
geom_col(show.legend = F) +
geom_label(aes(label = scales::percent(pct)), show.legend = F) +
labs(x = "Percentage",
y = "") +
theme(axis.text.y = element_text(hjust = 0))
- “Arranged or curated offers from buyers”Domestic transport infrastructure” takes 83%
Country and regions (WB category)
dat %>%
count(region, sort = T) %>%
mutate(pct = n / 18) %>%
mutate(region = fct_reorder(region, pct)) %>%
ggplot(aes(pct, region, fill = region)) +
geom_col(show.legend = F) +
geom_label(aes(label = scales::percent(pct)), show.legend = F) +
labs(x = "Percentage",
y = "") +
theme(axis.text.y = element_text(hjust = 0))
- “Sub-Saharan Africa” takes 61%
meta analysis
random-effects model (k = 18)
res_volume <- rma(yi, vi, method="DL", data=dat, slab= authors)
a <- forest.rma(res_volume,
slab=dat$authors,
order = "obs",
showweights = T,
header = T)
annotation <- sprintf("(Q = %.2f, df = %d, p = %.2f; I^2 = %.1f%%)",
res_volume$QE, res_volume$k - res_volume$p,
res_volume$QEp, res_volume$I2)
text(-3, -2, pos = 4, cex = 1, label = annotation)
using {Meta} package
library(meta)## Warning: 패키지 'meta'는 R 버전 4.2.3에서 작성되었습니다## Loading 'meta' package (version 6.5-0).
## Type 'help(meta)' for a brief overview.
## Readers of 'Meta-Analysis with R (Use R!)' should install
## older version of 'meta' package: https://tinyurl.com/dt4y5drsres_volume_meta <- metagen(TE = yi,
seTE = sqrt(vi),
data=dat,
sm = "SMD",
fixed = F,
random = T,
method.tau = "DL",
studlab = authors)
forest.meta(res_volume_meta,
sortvar = TE,
leftcols = c("studlab"),
rightcols = c("w.random","TE", "ci"),
col.study = "grey10",
col.square = "grey10",
weight.study = "common",
plotwidth = "14cm",
# colgap.forest = "-1cm",
squaresize = 1)
Sensitivity analysis
full and without-outlier analysis
Outlier detection
plot(influence(res_volume))
dat %>%
select(authors, row_id, vi, yi)## # A tibble: 18 × 4
## authors row_id vi yi
## <glue> <int> <dbl> <dbl>
## 1 Kingombe and Di Falco (2012) 1 0.00386 -0.109
## 2 Dumas and Játiva (2020) 2 0.00265 -0.0956
## 3 Khandker and Koolwal (2011) 3 0.00311 -0.0602
## 4 Bonilla et al. (2018) 4 0.00156 0.133
## 5 Khandker et al. (2009) 5 0.00372 0.176
## 6 Jirgi et al. (2019) 6 0.0258 0.486
## 7 Thomas et al. (2021) 7 0.00138 -0.0186
## 8 Fuller (2014) 8 0.00541 0.0541
## 9 Garbero and Songsermsawas (2018) 9 0.00222 -0.160
## 10 Lema et al. (2017) 10 0.0592 0.490
## 11 Skelley (2018) 11 0.0181 0.0264
## 12 Ministry of Foreign Affairs of Denmark (2010) 12 0.0509 0.0564
## 13 Wanjala and Muradia (2013) 13 0.0136 0.163
## 14 Sulaiman et al. (2021) 14 0.0243 1.78
## 15 Ring et al. (2018) 15 0.00285 0.106
## 16 Barnett et al. (2018) 16 0.00209 0.128
## 17 Gebresilasse (2018) 17 0.00203 -0.120
## 18 Torero et al. (2016) 18 0.00243 -0.0587- Study 14 (“Sulaiman et al. (2021)”) might be outlier depending on satatics
Analysis without the outlier (“Sulaiman et al. (2021)”) (k = 17)
res_volume_without_outlier <- rma(yi, vi, method="DL",
data= dat %>% filter(row_id != 14), slab = authors)
b <- forest(res_volume_without_outlier,
slab=dat %>% filter(row_id != 14) %>% select(authors) %>% pull,
order = "obs",
showweights = F,
header = T,
addpred= F)
Comparison between with and without outlier
comparision <-
tidy(res_volume, conf.int = T) %>%
mutate(type = "full model") %>%
bind_rows(tidy(res_volume_without_outlier, conf.int = T) %>%
mutate(type = "without-outlier model"))
comparision %>%
mutate(type = str_to_sentence(type)) %>%
mutate(type = fct_reorder(type, estimate)) %>%
ggplot(aes(x = estimate , y = type, color = type)) +
geom_errorbarh(aes(xmin = conf.low,
xmax = conf.high),
height = 0.3,
size = 2, show.legend = F) +
geom_point(size = 5, show.legend = F) +
geom_vline(xintercept = 0, lty = 2, size = 2, , show.legend = F) +
labs(y = "") +
theme(axis.text.y = element_text(size = 20))## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
- Without-outlier model include 0 in the confidence interval, although estimates are positive. It means estimates are sensitive to the outlier in terms of significance
Publication bias
Funnel plot
funnel(res_volume,
label= 1,
offset = -9,
legend = T)
At first glance, the funnel plot does not appear asymmetrical.
“When there is no publication bias, the data points in such a plot should form a roughly symmetrical, upside-down funnel” https://bookdown.org/MathiasHarrer/Doing_Meta_Analysis_in_R/pub-bias.html?
Study 14 (“Sulaiman et al. (2021)”) might be outlier depending on the figure
Asymmetry test (Egger’s Regression Test)
With a p-value lower than 0.05, we fail to reject the hypothesis that the funnel plot is asymmetric. This suggests the presence of publication bias.
” To evaluate the funnel asymmetry, we inspect the size of B0, and if it differs significantly from zero. When this is the case, Egger’s test indicates funnel plot asymmetry.
https://bookdown.org/MathiasHarrer/Doing_Meta_Analysis_in_R/pub-bias.html?#eggers-test
funnel plot without outlier
funnel(res_volume_without_outlier,
label= 1,
offset = -9,
legend = T)
regtest(res_volume_without_outlier, model = "lm")##
## Regression Test for Funnel Plot Asymmetry
##
## Model: weighted regression with multiplicative dispersion
## Predictor: standard error
##
## Test for Funnel Plot Asymmetry: t = 1.3100, df = 15, p = 0.2099
## Limit Estimate (as sei -> 0): b = -0.0810 (CI: -0.2357, 0.0736)- With a p-value greater than 0.05, we can reject the hypothesis that the funnel plot is asymmetric. In other words, there might not be publication bias, unlike the asymmetry analysis with full models.
Additional analysis (searching for moderator)
- There are no statistically significant moderators
res_exp <- rma(yi, vi, mods = ~ exposure, data=dat)## Warning: 2 studies with NAs omitted from model fitting.res_exp##
## Mixed-Effects Model (k = 16; tau^2 estimator: REML)
##
## tau^2 (estimated amount of residual heterogeneity): 0.0115 (SE = 0.0061)
## tau (square root of estimated tau^2 value): 0.1070
## I^2 (residual heterogeneity / unaccounted variability): 77.14%
## H^2 (unaccounted variability / sampling variability): 4.37
## R^2 (amount of heterogeneity accounted for): 0.00%
##
## Test for Residual Heterogeneity:
## QE(df = 14) = 63.2599, p-val < .0001
##
## Test of Moderators (coefficient 2):
## QM(df = 1) = 0.4398, p-val = 0.5072
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## intrcpt -0.0309 0.0744 -0.4151 0.6781 -0.1766 0.1149
## exposure 0.0007 0.0010 0.6631 0.5072 -0.0013 0.0026
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1res_year <- rma(yi, vi, mods = ~ year_pub, data=dat)
res_year##
## Mixed-Effects Model (k = 18; tau^2 estimator: REML)
##
## tau^2 (estimated amount of residual heterogeneity): 0.2653 (SE = 0.1566)
## tau (square root of estimated tau^2 value): 0.5151
## I^2 (residual heterogeneity / unaccounted variability): 98.91%
## H^2 (unaccounted variability / sampling variability): 91.84
## R^2 (amount of heterogeneity accounted for): 0.00%
##
## Test for Residual Heterogeneity:
## QE(df = 6) = 166.7210, p-val < .0001
##
## Test of Moderators (coefficients 2:12):
## QM(df = 11) = 5.4633, p-val = 0.9067
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## intrcpt 0.1758 0.5187 0.3389 0.7347 -0.8408 1.1924
## year_pub2010 -0.1194 0.7650 -0.1560 0.8760 -1.6188 1.3800
## year_pub2011 -0.2360 0.7331 -0.3219 0.7475 -1.6729 1.2008
## year_pub2012 -0.2845 0.7336 -0.3879 0.6981 -1.7224 1.1533
## year_pub2013 -0.0129 0.7402 -0.0174 0.9861 -1.4637 1.4379
## year_pub2014 -0.1217 0.7347 -0.1657 0.8684 -1.5616 1.3182
## year_pub2016 -0.2345 0.7326 -0.3201 0.7489 -1.6704 1.2015
## year_pub2017 0.3146 0.7704 0.4084 0.6830 -1.1954 1.8246
## year_pub2018 -0.1568 0.5604 -0.2799 0.7796 -1.2552 0.9415
## year_pub2019 0.3100 0.7484 0.4142 0.6787 -1.1568 1.7768
## year_pub2020 -0.2713 0.7328 -0.3703 0.7112 -1.7076 1.1649
## year_pub2021 0.6686 0.6386 1.0469 0.2951 -0.5831 1.9203
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1res_mci <- rma(yi, vi, mods = ~ eval_method, data=dat)
res_mci##
## Mixed-Effects Model (k = 18; tau^2 estimator: REML)
##
## tau^2 (estimated amount of residual heterogeneity): 0.1858 (SE = 0.0720)
## tau (square root of estimated tau^2 value): 0.4310
## I^2 (residual heterogeneity / unaccounted variability): 97.97%
## H^2 (unaccounted variability / sampling variability): 49.21
## R^2 (amount of heterogeneity accounted for): 0.00%
##
## Test for Residual Heterogeneity:
## QE(df = 15) = 199.1034, p-val < .0001
##
## Test of Moderators (coefficients 2:3):
## QM(df = 2) = 0.3296, p-val = 0.8481
##
## Model Results:
##
## estimate se zval pval ci.lb
## intrcpt 0.1873 0.1347 1.3903 0.1644 -0.0767
## eval_methodRDD -0.2460 0.4543 -0.5415 0.5882 -1.1363
## eval_methodStatistical matching -0.0638 0.2247 -0.2839 0.7765 -0.5041
## ci.ub
## intrcpt 0.4513
## eval_methodRDD 0.6443
## eval_methodStatistical matching 0.3765
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1res_region <- rma(yi, vi, mods = ~ region, data=dat)
res_region##
## Mixed-Effects Model (k = 18; tau^2 estimator: REML)
##
## tau^2 (estimated amount of residual heterogeneity): 0.1922 (SE = 0.0767)
## tau (square root of estimated tau^2 value): 0.4384
## I^2 (residual heterogeneity / unaccounted variability): 98.02%
## H^2 (unaccounted variability / sampling variability): 50.58
## R^2 (amount of heterogeneity accounted for): 0.00%
##
## Test for Residual Heterogeneity:
## QE(df = 14) = 185.3988, p-val < .0001
##
## Test of Moderators (coefficients 2:4):
## QM(df = 3) = 0.8588, p-val = 0.8354
##
## Model Results:
##
## estimate se zval pval ci.lb
## intrcpt -0.1597 0.4409 -0.3622 0.7172 -1.0239
## regionLatin America & Caribbean 0.3019 0.5196 0.5810 0.5613 -0.7166
## regionSouth Asia 0.2074 0.5111 0.4058 0.6849 -0.7943
## regionSub-Saharan Africa 0.3712 0.4611 0.8050 0.4208 -0.5325
## ci.ub
## intrcpt 0.7045
## regionLatin America & Caribbean 1.3203
## regionSouth Asia 1.2091
## regionSub-Saharan Africa 1.2748
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1res_inter_subcode <- rma(yi, vi, mods = ~ inter_subcode, data=dat)
res_inter_subcode##
## Mixed-Effects Model (k = 18; tau^2 estimator: REML)
##
## tau^2 (estimated amount of residual heterogeneity): 0.1989 (SE = 0.0800)
## tau (square root of estimated tau^2 value): 0.4460
## I^2 (residual heterogeneity / unaccounted variability): 97.86%
## H^2 (unaccounted variability / sampling variability): 46.84
## R^2 (amount of heterogeneity accounted for): 0.00%
##
## Test for Residual Heterogeneity:
## QE(df = 14) = 174.9045, p-val < .0001
##
## Test of Moderators (coefficients 2:4):
## QM(df = 3) = 0.5327, p-val = 0.9117
##
## Model Results:
##
## estimate
## intrcpt 0.1330
## inter_subcodeDomestic transport infrastructure 0.0442
## inter_subcodeDomestic transport infrastructure ; Alternative physical marketplaces -0.2927
## inter_subcodeDomestic transport infrastructure ; Storage deposit systems -0.0050
## se
## intrcpt 0.4477
## inter_subcodeDomestic transport infrastructure 0.4633
## inter_subcodeDomestic transport infrastructure ; Alternative physical marketplaces 0.6337
## inter_subcodeDomestic transport infrastructure ; Storage deposit systems 0.6336
## zval
## intrcpt 0.2970
## inter_subcodeDomestic transport infrastructure 0.0955
## inter_subcodeDomestic transport infrastructure ; Alternative physical marketplaces -0.4619
## inter_subcodeDomestic transport infrastructure ; Storage deposit systems -0.0079
## pval
## intrcpt 0.7665
## inter_subcodeDomestic transport infrastructure 0.9239
## inter_subcodeDomestic transport infrastructure ; Alternative physical marketplaces 0.6442
## inter_subcodeDomestic transport infrastructure ; Storage deposit systems 0.9937
## ci.lb
## intrcpt -0.7445
## inter_subcodeDomestic transport infrastructure -0.8637
## inter_subcodeDomestic transport infrastructure ; Alternative physical marketplaces -1.5347
## inter_subcodeDomestic transport infrastructure ; Storage deposit systems -1.2468
## ci.ub
## intrcpt 1.0105
## inter_subcodeDomestic transport infrastructure 0.9522
## inter_subcodeDomestic transport infrastructure ; Alternative physical marketplaces 0.9493
## inter_subcodeDomestic transport infrastructure ; Storage deposit systems 1.2368
##
## intrcpt
## inter_subcodeDomestic transport infrastructure
## inter_subcodeDomestic transport infrastructure ; Alternative physical marketplaces
## inter_subcodeDomestic transport infrastructure ; Storage deposit systems
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1