This is an R Markdown document. Markdown is a simple formatting syntax for authoring HTML, PDF, and MS Word documents. For more details on using R Markdown see http://rmarkdown.rstudio.com.
When you click the Knit button a document will be generated that includes both content as well as the output of any embedded R code chunks within the document. You can embed an R code chunk like this:
# Loading the Dataset
library(readxl)
Coding_Studies_Cleaned <- read_excel("Coding_Studies_Cleaned.xlsx")
#View(Coding_Studies_Cleaned)
#load metafor
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)
#rename the dataset
MetaData <- Coding_Studies_Cleaned
#RME
#extracting effect sizes
dat <- escalc(measure = "SMD", # effect size type
n1i = active_N, # sample size for group 1
m1i = active_rt, # mean for group 1
sd1i = active_SD, # sd for group 1
n2i = sham_N, # sample size for group 2
m2i = sham_rt, # mean for group 2
sd2i = sham_SD, # sd for group 2
data = MetaData, # data set name
var.names = c("g", "varg")) # create names for smd and variance columns
#Random Effects Model
res <- rma(yi = g, # the variable with the ES
vi = varg, # the variable with the var(ES)
data = dat, # data set to use
method = "REML", # method for RE model
slab = paste(dat$author, dat$year))
summary(res)
##
## Random-Effects Model (k = 71; tau^2 estimator: REML)
##
## logLik deviance AIC BIC AICc
## -41.3082 82.6164 86.6164 91.1134 86.7955
##
## tau^2 (estimated amount of total heterogeneity): 0.0691 (SE = 0.0305)
## tau (square root of estimated tau^2 value): 0.2629
## I^2 (total heterogeneity / total variability): 39.06%
## H^2 (total variability / sampling variability): 1.64
##
## Test for Heterogeneity:
## Q(df = 70) = 113.3603, p-val = 0.0008
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## -0.0953 0.0511 -1.8652 0.0622 -0.1954 0.0048 .
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##Including Plots
#Forest Plot
forest.rma(res)
#check for outliers
inf<-influence(res)
plot(inf)
#predicting interval
predict(res)
##
## pred se ci.lb ci.ub pi.lb pi.ub
## -0.0953 0.0511 -0.1954 0.0048 -0.6201 0.4296
#quartiles of g
qrtls <- quantile(dat$g, c(.25, .75))
# 3*IQR fences
fences <- qrtls + 3 * diff(qrtls) * c(-1, 1)
#data frame of fences
fence_dat <- data.frame(qrtl = qrtls, fence = fences)
#Outliers
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
# identify outliers
outliers <-
dat %>%
filter(g < fences[1] | g > fences[2])
# count outliers by analysis
n_outliers <-
outliers %>%
group_by(Study_ID) %>%
summarise(effects = n(), .groups = "drop_last") %>%
summarise(
studies = n(),
effects = sum(effects))
#Distribution plot the effect sizes with fences
library(ggplot2)
ggplot(dat, aes(g)) +
geom_density(fill = "cyan", alpha = 0.2) +
geom_rug(alpha = 0.1) +
geom_vline(data = fence_dat, aes(xintercept = qrtl), linetype = "solid") +
geom_vline(data = fence_dat, aes(xintercept = fence), linetype = "dashed") +
scale_x_continuous(breaks = seq(-6, 6, 2)) +
theme_minimal() +
labs(
x = "Effect size estimate",
y = "",
title = "Empirical distribution of effect size estimates"
)
#calculating the SE
dat <- dat %>%
mutate(SE_mod = sqrt(varg))
#Funnel Plot
ggplot(dat, aes(g, SE_mod, color = factor(Study_ID))) +
geom_hline(yintercept = 0) +
geom_vline(xintercept = 0) +
geom_vline(data = fence_dat, aes(xintercept = fence), linetype = "dashed") +
geom_abline(
intercept = 0,
slope = 1 / qnorm(c(0.005, 0.025, 0.05, 0.95, 0.975, 0.995)),
color = "grey"
) +
scale_x_continuous(breaks = seq(-6, 6, 2)) +
scale_y_reverse(expand = expansion(add = c(0.02, 0))) +
scale_color_discrete(guide = "none") +
geom_point(alpha = 0.2) +
theme_minimal() +
labs(
y = "Modified Standard Error",
x = "Standardized Mean Difference (Hedges' g)",
title = "Funnel plot with guide lines at significance levels of alpha = .01, .05, and .10",
subtitle = "Dashed vertical lines indicate outlier fences"
)
#install.packages("kableExtra")
library(kableExtra)
##
## Attaching package: 'kableExtra'
## The following object is masked from 'package:dplyr':
##
## group_rows
#table of outlier effects
outliers_tbl <- outliers %>%
mutate(SE = sqrt(varg)) %>%
transmute(
Author = author,
Year = year,
g,
SE,
N_active = active_N,
M_active = active_rt,
SD_active = active_SD,
N_sham = sham_N,
M_sham = sham_rt,
SD_sham = sham_SD
) %>%
arrange(desc(abs(g)))
#formatted table
kable(outliers_tbl,
caption = "Effect size estimates identified as outliers") %>%
kable_styling(full_width = FALSE)
| Author | Year | g | SE | N_active | M_active | SD_active | N_sham | M_sham | SD_sham |
|---|---|---|---|---|---|---|---|---|---|
| NA | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| :—— | —-: | –: | –: | ——–: | ——–: | ———: | ——: | ——: | ——-: |
#the output was a table with NAs. Indicating no outliers outside of the fences
#find the outlier
res$slab[29]
## [1] "Wei et al. 2024.3"
#removing the outlier
leave1out(res)
##
## estimate se zval pval ci.lb ci.ub
## -Meng et al. 2024.1 -0.0806 0.0495 -1.6268 0.1038 -0.1777 0.0165
## -Meng et al. 2024.2 -0.0899 0.0514 -1.7514 0.0799 -0.1906 0.0107
## -Meng et al. 2024.3 -0.0927 0.0519 -1.7868 0.0740 -0.1944 0.0090
## -Sipolins 2016 -0.0917 0.0523 -1.7529 0.0796 -0.1942 0.0108
## -Meng et al. 2023.1 -0.0822 0.0500 -1.6420 0.1006 -0.1802 0.0159
## -Meng et al. 2023.2 -0.0928 0.0519 -1.7893 0.0736 -0.1945 0.0089
## -Meng et al. 2023.3 -0.0937 0.0516 -1.8146 0.0696 -0.1949 0.0075
## -Alsultan et al. 2020.1 -0.0961 0.0514 -1.8683 0.0617 -0.1968 0.0047
## -Alsultan et al. 2020.2 -0.0977 0.0515 -1.8982 0.0577 -0.1985 0.0032
## -Alsultan et al. 2020.3 -0.0961 0.0515 -1.8679 0.0618 -0.1970 0.0047
## -Lu et al. 2020.1 -0.0870 0.0511 -1.7011 0.0889 -0.1872 0.0132
## -Lu et al. 2020.2 -0.0889 0.0514 -1.7275 0.0841 -0.1897 0.0120
## -Lu et al. 2020.3 -0.0934 0.0519 -1.8014 0.0716 -0.1950 0.0082
## -Lu et al. 2020.4 -0.1012 0.0515 -1.9653 0.0494 -0.2022 -0.0003
## -Lu et al. 2020.5 -0.0873 0.0512 -1.7057 0.0881 -0.1877 0.0130
## -Lu et al. 2020.6 -0.0839 0.0505 -1.6633 0.0963 -0.1828 0.0150
## -Lu et al. 2020.7 -0.1031 0.0512 -2.0116 0.0443 -0.2035 -0.0026
## -Lu et al. 2020.8 -0.1031 0.0512 -2.0131 0.0441 -0.2035 -0.0027
## -Lu et al. 2020.9 -0.0970 0.0519 -1.8699 0.0615 -0.1986 0.0047
## -Lu et al. 2020.10 -0.0944 0.0519 -1.8201 0.0687 -0.1961 0.0073
## -Lu et al. 2020.11 -0.0884 0.0514 -1.7213 0.0852 -0.1891 0.0123
## -Lu et al. 2020.12 -0.0919 0.0518 -1.7747 0.0759 -0.1933 0.0096
## -Lu et al. 2020.13 -0.0935 0.0519 -1.8024 0.0715 -0.1951 0.0082
## -Lu et al. 2020.14 -0.0880 0.0513 -1.7153 0.0863 -0.1886 0.0126
## -Lu et al. 2020.15 -0.0994 0.0517 -1.9230 0.0545 -0.2008 0.0019
## -Lu et al. 2020.16 -0.0913 0.0517 -1.7649 0.0776 -0.1926 0.0101
## -Wei et al. 2024.1 -0.1105 0.0493 -2.2400 0.0251 -0.2071 -0.0138
## -Wei et al. 2024.2 -0.1042 0.0510 -2.0411 0.0412 -0.2042 -0.0041
## -Wei et al. 2024.3 -0.0772 0.0481 -1.6025 0.1090 -0.1715 0.0172
## -Schmidt et al. 2024.1 -0.0934 0.0517 -1.8059 0.0709 -0.1948 0.0080
## -Schmidt et al. 2024.2 -0.0937 0.0517 -1.8108 0.0702 -0.1951 0.0077
## -Schmidt et al. 2024.3 -0.0937 0.0518 -1.8095 0.0704 -0.1951 0.0078
## -Schmidt et al. 2024.4 -0.0937 0.0517 -1.8108 0.0702 -0.1951 0.0077
## -Schmidt et al. 2024.5 -0.0965 0.0518 -1.8646 0.0622 -0.1979 0.0049
## -Schmidt et al. 2024.6 -0.0939 0.0518 -1.8151 0.0695 -0.1954 0.0075
## -Schmidt et al. 2024.7 -0.1014 0.0514 -1.9733 0.0485 -0.2020 -0.0007
## -Schmidt et al. 2024.8 -0.0987 0.0516 -1.9108 0.0560 -0.1999 0.0025
## -Schmidt et al. 2024.9 -0.0994 0.0516 -1.9253 0.0542 -0.2005 0.0018
## -Schmidt et al. 2024.10 -0.1028 0.0512 -2.0091 0.0445 -0.2030 -0.0025
## -Schmidt et al. 2024.11 -0.0991 0.0516 -1.9208 0.0548 -0.2003 0.0020
## -Schmidt et al. 2024.12 -0.1021 0.0513 -1.9928 0.0463 -0.2026 -0.0017
## -Schmidt et al. 2024.13 -0.0962 0.0518 -1.8585 0.0631 -0.1977 0.0052
## -Schmidt et al. 2024.14 -0.0956 0.0518 -1.8463 0.0649 -0.1970 0.0059
## -Schmidt et al. 2024.15 -0.1006 0.0515 -1.9545 0.0506 -0.2015 0.0003
## -Schmidt et al. 2024.16 -0.0974 0.0517 -1.8821 0.0598 -0.1987 0.0040
## -Schmidt et al. 2024.17 -0.1007 0.0515 -1.9563 0.0504 -0.2016 0.0002
## -Schmidt et al. 2024.18 -0.1013 0.0514 -1.9706 0.0488 -0.2020 -0.0005
## -Khanmohammadi et al. 2025.1 -0.0946 0.0519 -1.8223 0.0684 -0.1963 0.0071
## -Khanmohammadi et al. 2025.2 -0.0932 0.0519 -1.7966 0.0724 -0.1948 0.0085
## -Khanmohammadi et al. 2025.3 -0.0955 0.0519 -1.8401 0.0657 -0.1973 0.0062
## -Fujiyama et al. 2025.1 -0.1037 0.0510 -2.0341 0.0419 -0.2036 -0.0038
## -Fujiyama et al. 2025.2 -0.0998 0.0515 -1.9362 0.0528 -0.2008 0.0012
## -Bashir et al. 2020.1 -0.0890 0.0513 -1.7343 0.0829 -0.1895 0.0116
## -Bashir et al. 2020.2 -0.0913 0.0516 -1.7709 0.0766 -0.1924 0.0098
## -Friehs et al. 2021 -0.0900 0.0517 -1.7425 0.0814 -0.1913 0.0112
## -McDermott et al. 2019.1 -0.0951 0.0518 -1.8368 0.0662 -0.1967 0.0064
## -McDermott et al. 2019.2 -0.0951 0.0518 -1.8361 0.0663 -0.1966 0.0064
## -McDermott et al. 2019.3 -0.1064 0.0504 -2.1102 0.0348 -0.2052 -0.0076
## -McDermott et al. 2019.4 -0.1060 0.0505 -2.0966 0.0360 -0.2050 -0.0069
## -Miler et al. 2018.1 -0.0873 0.0511 -1.7098 0.0873 -0.1874 0.0128
## -Miler et al. 2018.2 -0.0973 0.0517 -1.8814 0.0599 -0.1987 0.0041
## -Miler et al. 2018.3 -0.0956 0.0518 -1.8470 0.0648 -0.1971 0.0059
## -Brambilla et al. 2021.1 -0.0903 0.0515 -1.7526 0.0797 -0.1913 0.0107
## -Brambilla et al. 2021.2 -0.0947 0.0518 -1.8274 0.0676 -0.1962 0.0069
## -Hausman 2024 -0.1029 0.0520 -1.9814 0.0476 -0.2048 -0.0011
## -Diedrich et al. 2025.1 -0.0963 0.0522 -1.8438 0.0652 -0.1986 0.0061
## -Diedrich et al. 2025.2 -0.0998 0.0520 -1.9196 0.0549 -0.2018 0.0021
## -Diedrich et al. 2025.3 -0.0961 0.0522 -1.8401 0.0658 -0.1984 0.0063
## -Souki et al. 2025.1 -0.0939 0.0517 -1.8181 0.0690 -0.1951 0.0073
## -Souki et al. 2025.2 -0.0918 0.0515 -1.7815 0.0748 -0.1928 0.0092
## -Souki et al. 2025.3 -0.1006 0.0514 -1.9583 0.0502 -0.2013 0.0001
## Q Qp tau2 I2 H2
## -Meng et al. 2024.1 104.3628 0.0039 0.0562 34.2565 1.5211
## -Meng et al. 2024.2 111.9392 0.0008 0.0691 39.1980 1.6447
## -Meng et al. 2024.3 113.0737 0.0007 0.0718 39.9016 1.6639
## -Sipolins 2016 112.7088 0.0007 0.0731 39.6102 1.6559
## -Meng et al. 2023.1 106.2189 0.0027 0.0595 35.5333 1.5512
## -Meng et al. 2023.2 113.1050 0.0006 0.0718 39.9197 1.6644
## -Meng et al. 2023.3 113.2374 0.0006 0.0711 39.8540 1.6626
## -Alsultan et al. 2020.1 113.3185 0.0006 0.0703 39.7080 1.6586
## -Alsultan et al. 2020.2 113.0298 0.0007 0.0703 39.6531 1.6571
## -Alsultan et al. 2020.3 113.3173 0.0006 0.0705 39.7395 1.6595
## -Lu et al. 2020.1 110.4166 0.0011 0.0669 38.3001 1.6207
## -Lu et al. 2020.2 111.6017 0.0009 0.0690 39.0146 1.6397
## -Lu et al. 2020.3 113.2058 0.0006 0.0718 39.9548 1.6654
## -Lu et al. 2020.4 111.9142 0.0008 0.0694 39.1527 1.6435
## -Lu et al. 2020.5 110.6538 0.0011 0.0673 38.4442 1.6245
## -Lu et al. 2020.6 107.7997 0.0020 0.0624 36.6835 1.5794
## -Lu et al. 2020.7 110.8741 0.0010 0.0676 38.5121 1.6263
## -Lu et al. 2020.8 110.8348 0.0010 0.0675 38.4877 1.6257
## -Lu et al. 2020.9 113.2483 0.0006 0.0718 39.9660 1.6657
## -Lu et al. 2020.10 113.3269 0.0006 0.0720 40.0222 1.6673
## -Lu et al. 2020.11 111.3588 0.0009 0.0686 38.8692 1.6358
## -Lu et al. 2020.12 112.8593 0.0007 0.0712 39.7558 1.6599
## -Lu et al. 2020.13 113.2149 0.0006 0.0718 39.9600 1.6656
## -Lu et al. 2020.14 111.1047 0.0010 0.0681 38.7165 1.6318
## -Lu et al. 2020.15 112.6628 0.0007 0.0708 39.6091 1.6559
## -Lu et al. 2020.16 112.6707 0.0007 0.0709 39.6462 1.6569
## -Wei et al. 2024.1 103.6051 0.0045 0.0549 33.7576 1.5096
## -Wei et al. 2024.2 110.1157 0.0012 0.0662 38.0323 1.6137
## -Wei et al. 2024.3 98.5927 0.0112 0.0477 30.7192 1.4434
## -Schmidt et al. 2024.1 113.2003 0.0006 0.0714 39.8923 1.6637
## -Schmidt et al. 2024.2 113.2431 0.0006 0.0714 39.9153 1.6643
## -Schmidt et al. 2024.3 113.2373 0.0006 0.0715 39.9202 1.6645
## -Schmidt et al. 2024.4 113.2431 0.0006 0.0714 39.9153 1.6643
## -Schmidt et al. 2024.5 113.2986 0.0006 0.0715 39.9373 1.6649
## -Schmidt et al. 2024.6 113.2745 0.0006 0.0715 39.9321 1.6648
## -Schmidt et al. 2024.7 111.7306 0.0009 0.0689 39.0425 1.6405
## -Schmidt et al. 2024.8 112.8586 0.0007 0.0707 39.6844 1.6579
## -Schmidt et al. 2024.9 112.6453 0.0007 0.0704 39.5673 1.6547
## -Schmidt et al. 2024.10 110.8630 0.0010 0.0674 38.5491 1.6273
## -Schmidt et al. 2024.11 112.7168 0.0007 0.0705 39.6036 1.6557
## -Schmidt et al. 2024.12 111.2762 0.0010 0.0681 38.7840 1.6336
## -Schmidt et al. 2024.13 113.3261 0.0006 0.0715 39.9536 1.6654
## -Schmidt et al. 2024.14 113.3574 0.0006 0.0716 39.9731 1.6659
## -Schmidt et al. 2024.15 112.1263 0.0008 0.0695 39.2676 1.6466
## -Schmidt et al. 2024.16 113.1781 0.0006 0.0713 39.8674 1.6630
## -Schmidt et al. 2024.17 112.0925 0.0008 0.0695 39.2495 1.6461
## -Schmidt et al. 2024.18 111.7901 0.0009 0.0689 39.0763 1.6414
## -Khanmohammadi et al. 2025.1 113.3367 0.0006 0.0721 40.0341 1.6676
## -Khanmohammadi et al. 2025.2 113.1663 0.0006 0.0718 39.9387 1.6650
## -Khanmohammadi et al. 2025.3 113.3590 0.0006 0.0721 40.0435 1.6679
## -Fujiyama et al. 2025.1 110.1780 0.0012 0.0663 38.1595 1.6171
## -Fujiyama et al. 2025.2 112.4663 0.0007 0.0701 39.4610 1.6518
## -Bashir et al. 2020.1 111.4807 0.0009 0.0685 38.9522 1.6381
## -Bashir et al. 2020.2 112.6248 0.0007 0.0703 39.5675 1.6547
## -Friehs et al. 2021 112.2065 0.0008 0.0703 39.3910 1.6499
## -McDermott et al. 2019.1 113.3590 0.0006 0.0717 39.9923 1.6665
## -McDermott et al. 2019.2 113.3583 0.0006 0.0717 39.9921 1.6664
## -McDermott et al. 2019.3 107.8297 0.0019 0.0625 36.7612 1.5813
## -McDermott et al. 2019.4 108.2923 0.0018 0.0632 37.0319 1.5881
## -Miler et al. 2018.1 110.4154 0.0011 0.0669 38.3584 1.6223
## -Miler et al. 2018.2 113.1836 0.0006 0.0713 39.8705 1.6631
## -Miler et al. 2018.3 113.3565 0.0006 0.0716 39.9725 1.6659
## -Brambilla et al. 2021.1 112.2588 0.0008 0.0699 39.3840 1.6497
## -Brambilla et al. 2021.2 113.3409 0.0006 0.0717 39.9857 1.6663
## -Hausman 2024 110.2742 0.0012 0.0701 38.2445 1.6193
## -Diedrich et al. 2025.1 113.3257 0.0006 0.0731 40.0602 1.6683
## -Diedrich et al. 2025.2 112.5383 0.0007 0.0717 39.5797 1.6551
## -Diedrich et al. 2025.3 113.3386 0.0006 0.0732 40.0687 1.6686
## -Souki et al. 2025.1 113.2650 0.0006 0.0711 39.8677 1.6630
## -Souki et al. 2025.2 112.7648 0.0007 0.0704 39.6153 1.6560
## -Souki et al. 2025.3 112.0095 0.0008 0.0692 39.2012 1.6448
#Robust variance estimate
#As some of our studies are correlated together, we need to take the correlation into account. So we're using RVE and we assume r = 0.60.
#loading the package
library(clubSandwich)
## Registered S3 method overwritten by 'clubSandwich':
## method from
## bread.mlm sandwich
##Create variance-covariance matrix based on known variances and assumed correlation
V_list <- impute_covariance_matrix(vi = dat$varg, #known correlation vector
cluster = dat$Study_ID, #Clustering variable
r = 0.60) #assumed correlation
## Warning: `impute_covariance_matrix()` was deprecated in clubSandwich 0.5.11.
## ℹ Please use `metafor::vcalc()` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
##Create a column for Effect Size ID in order to indicate nesting structure
dat$ES_ID <- 1:nrow(dat)
##Build the CHE model
SummaryModel <- rma.mv(yi=g, #effect size
V = V_list, #covariance matrix
random = ~1 | Study_ID/ES_ID, #nesting structure: highest / lowest
test= "t", #use t-tests
data=dat, #define data
method="REML") #estimate variances using REML
summary(SummaryModel)
##
## Multivariate Meta-Analysis Model (k = 71; method: REML)
##
## logLik Deviance AIC BIC AICc
## -40.4319 80.8638 86.8638 93.6093 87.2274
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0000 0.0000 17 no Study_ID
## sigma^2.2 0.1111 0.3333 71 no Study_ID/ES_ID
##
## Test for Heterogeneity:
## Q(df = 70) = 207.0690, p-val < .0001
##
## Model Results:
##
## estimate se tval df pval ci.lb ci.ub
## -0.1176 0.0859 -1.3690 70 0.1754 -0.2888 0.0537
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##Above model-based standard errors and hypothesis tests are likely incorrect, so we conduct hypothesis testing with Robust Variance Estimation using clubsandwich.
SummaryModelcf <- coef_test(SummaryModel,#estimation model above
cluster=dat$Study_ID, #define cluster IDs
vcov = "CR2") #estimation method (CR2 is best)
SummaryModelcf
## Alternative hypothesis: two-sided
## Coef. Estimate SE Null value t-stat d.f. (Satt) p-val (Satt) Sig.
## intrcpt -0.118 0.0671 0 -1.75 14.4 0.101
#Subgroup Moderator, Stimulation Type
#Categorical moderator --> stimulation type
dat$stimulation_method <- as.factor(dat$stimulation_method)
##table to see the representation in each category
table(dat$stimulation_method)
##
## CT + tRNS HD-tDCS Single-Pulse TMS tACS
## 2 23 2 11
## tDCS tDCS_Anodal tDCS_Cathodal
## 11 11 11
##boxplots of effect sizes by the categorical moderator
boxplot(dat$g ~ dat$stimulation_method)
#model for the res with moderator
res_sti<-rma(yi=g, vi=varg, mods = ~ stimulation_method, data = dat)
res_sti
##
## Mixed-Effects Model (k = 71; tau^2 estimator: REML)
##
## tau^2 (estimated amount of residual heterogeneity): 0.0543 (SE = 0.0291)
## tau (square root of estimated tau^2 value): 0.2329
## I^2 (residual heterogeneity / unaccounted variability): 33.18%
## H^2 (unaccounted variability / sampling variability): 1.50
## R^2 (amount of heterogeneity accounted for): 21.48%
##
## Test for Residual Heterogeneity:
## QE(df = 64) = 94.4454, p-val = 0.0080
##
## Test of Moderators (coefficients 2:7):
## QM(df = 6) = 12.8189, p-val = 0.0460
##
## Model Results:
##
## estimate se zval pval ci.lb
## intrcpt -0.2994 0.3001 -0.9979 0.3183 -0.8875
## stimulation_methodHD-tDCS 0.0951 0.3114 0.3052 0.7602 -0.5153
## stimulation_methodSingle-Pulse TMS -0.2167 0.4371 -0.4957 0.6201 -1.0734
## stimulation_methodtACS 0.0929 0.3218 0.2887 0.7728 -0.5378
## stimulation_methodtDCS 0.2535 0.3258 0.7782 0.4364 -0.3850
## stimulation_methodtDCS_Anodal 0.5489 0.3274 1.6765 0.0936 -0.0928
## stimulation_methodtDCS_Cathodal 0.3270 0.3270 1.0000 0.3173 -0.3139
## ci.ub
## intrcpt 0.2887
## stimulation_methodHD-tDCS 0.7054
## stimulation_methodSingle-Pulse TMS 0.6400
## stimulation_methodtACS 0.7236
## stimulation_methodtDCS 0.8920
## stimulation_methodtDCS_Anodal 1.1907 .
## stimulation_methodtDCS_Cathodal 0.9678
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#no intercept
res_sti_z<-rma(yi=g, vi=varg, mods = ~ stimulation_method -1, data = dat)
res_sti_z
##
## Mixed-Effects Model (k = 71; tau^2 estimator: REML)
##
## tau^2 (estimated amount of residual heterogeneity): 0.0543 (SE = 0.0291)
## tau (square root of estimated tau^2 value): 0.2329
## I^2 (residual heterogeneity / unaccounted variability): 33.18%
## H^2 (unaccounted variability / sampling variability): 1.50
##
## Test for Residual Heterogeneity:
## QE(df = 64) = 94.4454, p-val = 0.0080
##
## Test of Moderators (coefficients 1:7):
## QM(df = 7) = 16.6194, p-val = 0.0200
##
## Model Results:
##
## estimate se zval pval ci.lb
## stimulation_methodCT + tRNS -0.2994 0.3001 -0.9979 0.3183 -0.8875
## stimulation_methodHD-tDCS -0.2044 0.0832 -2.4556 0.0141 -0.3675
## stimulation_methodSingle-Pulse TMS -0.5161 0.3179 -1.6238 0.1044 -1.1391
## stimulation_methodtACS -0.2065 0.1163 -1.7763 0.0757 -0.4344
## stimulation_methodtDCS -0.0459 0.1269 -0.3619 0.7175 -0.2945
## stimulation_methodtDCS_Anodal 0.2495 0.1310 1.9042 0.0569 -0.0073
## stimulation_methodtDCS_Cathodal 0.0275 0.1299 0.2120 0.8321 -0.2270
## ci.ub
## stimulation_methodCT + tRNS 0.2887
## stimulation_methodHD-tDCS -0.0412 *
## stimulation_methodSingle-Pulse TMS 0.1069
## stimulation_methodtACS 0.0214 .
## stimulation_methodtDCS 0.2027
## stimulation_methodtDCS_Anodal 0.5064 .
## stimulation_methodtDCS_Cathodal 0.2821
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#testing contrasts
anova(res_sti_z, btt = 2:6)
##
## Test of Moderators (coefficients 2:6):
## QM(df = 5) = 15.5787, p-val = 0.0082
#plot with the subgroup summary effect
forest.default(res_sti_z$b, sei = res_sti_z$se, slab = c("CT & tRNS", "HD-tDCS", "Single-Pulse TMS", "tACS", "tDCS", "tDCS_Anodal", "tDCS_Cathodal"), header = "Forest Plot with Stimulation Type Summary Effects", cex = 1)
#Subgroup Moderator2, Attention Function
#Categorical moderator --> The Attention Function
dat$Attention_Function <- as.factor(dat$Attention_Function)
##table to see the representation in each category
table(dat$Attention_Function)
##
## alerting Attention Flexibility Attention Switching
## 10 6 3
## Attentional Impulsiveness Divided Attention executive control
## 11 6 14
## orienting response inhibition selective attention
## 8 3 10
##boxplots of effect sizes by the categorical moderator
boxplot(dat$g ~ dat$Attention_Function)
#model for the res with moderator
res_ATF<-rma(yi=g, vi=varg, mods = ~ Attention_Function, data = dat)
res_ATF
##
## Mixed-Effects Model (k = 71; tau^2 estimator: REML)
##
## tau^2 (estimated amount of residual heterogeneity): 0.0488 (SE = 0.0285)
## tau (square root of estimated tau^2 value): 0.2210
## I^2 (residual heterogeneity / unaccounted variability): 30.95%
## H^2 (unaccounted variability / sampling variability): 1.45
## R^2 (amount of heterogeneity accounted for): 29.35%
##
## Test for Residual Heterogeneity:
## QE(df = 62) = 89.4233, p-val = 0.0129
##
## Test of Moderators (coefficients 2:9):
## QM(df = 8) = 16.8513, p-val = 0.0317
##
## Model Results:
##
## estimate se zval pval
## intrcpt -0.0787 0.1249 -0.6304 0.5284
## Attention_FunctionAttention Flexibility 0.2402 0.2145 1.1198 0.2628
## Attention_FunctionAttention Switching 0.1679 0.3580 0.4691 0.6390
## Attention_FunctionAttentional Impulsiveness -0.2936 0.1762 -1.6667 0.0956
## Attention_FunctionDivided Attention 0.2784 0.2146 1.2973 0.1945
## Attention_Functionexecutive control -0.2058 0.1658 -1.2414 0.2145
## Attention_Functionorienting 0.2412 0.1865 1.2929 0.1961
## Attention_Functionresponse inhibition -0.0698 0.2556 -0.2732 0.7847
## Attention_Functionselective attention -0.0069 0.1705 -0.0405 0.9677
## ci.lb ci.ub
## intrcpt -0.3235 0.1660
## Attention_FunctionAttention Flexibility -0.1802 0.6605
## Attention_FunctionAttention Switching -0.5337 0.8695
## Attention_FunctionAttentional Impulsiveness -0.6389 0.0517 .
## Attention_FunctionDivided Attention -0.1422 0.6991
## Attention_Functionexecutive control -0.5308 0.1192
## Attention_Functionorienting -0.1244 0.6068
## Attention_Functionresponse inhibition -0.5707 0.4311
## Attention_Functionselective attention -0.3411 0.3273
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#no intercept
res_ATF_z<-rma(yi=g, vi=varg, mods = ~ Attention_Function -1, data = dat)
res_ATF_z
##
## Mixed-Effects Model (k = 71; tau^2 estimator: REML)
##
## tau^2 (estimated amount of residual heterogeneity): 0.0488 (SE = 0.0285)
## tau (square root of estimated tau^2 value): 0.2210
## I^2 (residual heterogeneity / unaccounted variability): 30.95%
## H^2 (unaccounted variability / sampling variability): 1.45
##
## Test for Residual Heterogeneity:
## QE(df = 62) = 89.4233, p-val = 0.0129
##
## Test of Moderators (coefficients 1:9):
## QM(df = 9) = 20.7834, p-val = 0.0136
##
## Model Results:
##
## estimate se zval pval
## Attention_Functionalerting -0.0787 0.1249 -0.6304 0.5284
## Attention_FunctionAttention Flexibility 0.1614 0.1744 0.9259 0.3545
## Attention_FunctionAttention Switching 0.0892 0.3355 0.2659 0.7903
## Attention_FunctionAttentional Impulsiveness -0.3724 0.1243 -2.9968 0.0027
## Attention_FunctionDivided Attention 0.1997 0.1745 1.1441 0.2526
## Attention_Functionexecutive control -0.2846 0.1091 -2.6088 0.0091
## Attention_Functionorienting 0.1624 0.1386 1.1723 0.2411
## Attention_Functionresponse inhibition -0.1486 0.2230 -0.6663 0.5052
## Attention_Functionselective attention -0.0856 0.1161 -0.7377 0.4607
## ci.lb ci.ub
## Attention_Functionalerting -0.3235 0.1660
## Attention_FunctionAttention Flexibility -0.1803 0.5032
## Attention_FunctionAttention Switching -0.5683 0.7467
## Attention_FunctionAttentional Impulsiveness -0.6159 -0.1288 **
## Attention_FunctionDivided Attention -0.1424 0.5418
## Attention_Functionexecutive control -0.4984 -0.0708 **
## Attention_Functionorienting -0.1091 0.4340
## Attention_Functionresponse inhibition -0.5856 0.2885
## Attention_Functionselective attention -0.3132 0.1419
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#testing contrasts
anova(res_ATF_z, btt = c(4,6))
##
## Test of Moderators (coefficients 4,6):
## QM(df = 2) = 15.7865, p-val = 0.0004
#plot with the subgroup summary effect
forest.default(res_ATF_z$b, sei = res_ATF_z$se, slab = c("Alerting", "Flexibility", "Switching", "Impulsiveness", "Divided Attention", "Executive Control", "Orienting", "Response Inhibition", "Selective Attention"), header = "Forest Plot with Attention Function Summary Effects", cex = 1)
#meta_Regression for Age variable
#treating age as numeric
dat$total_age_mean <- as.numeric(dat$total_age_mean)
## Warning: NAs introduced by coercion
#mean_centering
dat$total_age_mean_centered <-dat$total_age_mean- mean(dat$total_age_mean)
#random effects model of Age
res_age<-rma(yi=g, vi=varg, mods = ~ total_age_mean, data = dat)
## Warning: 3 studies with NAs omitted from model fitting.
res_age
##
## Mixed-Effects Model (k = 68; tau^2 estimator: REML)
##
## tau^2 (estimated amount of residual heterogeneity): 0.0604 (SE = 0.0295)
## tau (square root of estimated tau^2 value): 0.2457
## I^2 (residual heterogeneity / unaccounted variability): 36.10%
## H^2 (unaccounted variability / sampling variability): 1.56
## R^2 (amount of heterogeneity accounted for): 16.76%
##
## Test for Residual Heterogeneity:
## QE(df = 66) = 102.4851, p-val = 0.0027
##
## Test of Moderators (coefficient 2):
## QM(df = 1) = 4.6931, p-val = 0.0303
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## intrcpt -0.2935 0.1044 -2.8126 0.0049 -0.4980 -0.0890 **
## total_age_mean 0.0048 0.0022 2.1664 0.0303 0.0005 0.0091 *
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##Draw a bubble plot
size <- 1/dat$varg
size <- size/max(size)
plot(dat$total_age_mean, dat$g,
cex = size*20,
xlab = "Age",
ylab = "Effect size",
main = "Association between NIBS Effects and Age")
abline(res_age)
b <- coef(res_age) # b[1] = intercept, b[2] = slope
abline(a = b[1], b = b[2])
#Subgroup Moderator3, Age categories
#As the age is varied between the two age groups of young adults and old adults only, we categorized it to two separate groups and ran another subgroup analysis on it
#create a grouping variable
dat$age_group <- ifelse(dat$total_age_mean < 30,
"young_adults",
"older_adults")
#Categorical moderator --> The Age groups
dat$age_group <- factor(dat$age_group, levels = c("young_adults", "older_adults"))
#how many studies in each group?
table(dat$age_group)
##
## young_adults older_adults
## 41 27
#summary of each group
aggregate(total_age_mean ~ age_group, data = dat, summary)
## age_group total_age_mean.Min. total_age_mean.1st Qu. total_age_mean.Median
## 1 young_adults 21.15000 21.15000 21.37000
## 2 older_adults 67.70000 67.70000 67.70000
## total_age_mean.Mean total_age_mean.3rd Qu. total_age_mean.Max.
## 1 22.19073 22.90000 26.27000
## 2 68.33741 68.96000 71.25000
##boxplots of effect sizes by the categorical moderator
boxplot(dat$g ~ dat$age_group)
#model for the res with moderator
res_Age2<-rma(yi=g, vi=varg, mods = ~ age_group, data = dat)
## Warning: 3 studies with NAs omitted from model fitting.
res_Age2
##
## Mixed-Effects Model (k = 68; tau^2 estimator: REML)
##
## tau^2 (estimated amount of residual heterogeneity): 0.0610 (SE = 0.0297)
## tau (square root of estimated tau^2 value): 0.2469
## I^2 (residual heterogeneity / unaccounted variability): 36.34%
## H^2 (unaccounted variability / sampling variability): 1.57
## R^2 (amount of heterogeneity accounted for): 15.94%
##
## Test for Residual Heterogeneity:
## QE(df = 66) = 102.8458, p-val = 0.0025
##
## Test of Moderators (coefficient 2):
## QM(df = 1) = 4.4776, p-val = 0.0343
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## intrcpt -0.1854 0.0661 -2.8069 0.0050 -0.3149 -0.0560 **
## age_groupolder_adults 0.2180 0.1030 2.1160 0.0343 0.0161 0.4200 *
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#no intercept
res_Age2_z<-rma(yi=g, vi=varg, mods = ~ age_group -1, data = dat)
## Warning: 3 studies with NAs omitted from model fitting.
res_Age2_z
##
## Mixed-Effects Model (k = 68; tau^2 estimator: REML)
##
## tau^2 (estimated amount of residual heterogeneity): 0.0610 (SE = 0.0297)
## tau (square root of estimated tau^2 value): 0.2469
## I^2 (residual heterogeneity / unaccounted variability): 36.34%
## H^2 (unaccounted variability / sampling variability): 1.57
##
## Test for Residual Heterogeneity:
## QE(df = 66) = 102.8458, p-val = 0.0025
##
## Test of Moderators (coefficients 1:2):
## QM(df = 2) = 8.0484, p-val = 0.0179
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## age_groupyoung_adults -0.1854 0.0661 -2.8069 0.0050 -0.3149 -0.0560 **
## age_groupolder_adults 0.0326 0.0791 0.4122 0.6802 -0.1224 0.1876
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#testing contrasts
anova(res_Age2_z, btt = 1:2)
##
## Test of Moderators (coefficients 1:2):
## QM(df = 2) = 8.0484, p-val = 0.0179
#plot with the subgroup summary effect
forest.default(res_Age2_z$b, sei = res_Age2_z$se, slab = c("Young Adults", "Old Adults"), cex = 1)