R Markdown

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)
Effect size estimates identified as outliers
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)