Photo by Helloquence on Unsplash

Workshop: A Practical Guide for Meta-Analyses in R

Overview

About this file

This blogpost is an adaption from my 90 minute workshop on all basic steps in conducting a meta-analysis in R. You can find the data set and code here.

The audience of the original workshop were graduate psychology students at the University of Vienna with a focus on research methods. Most statistical and methodical elements were already discussed, yet the implementation in R was new to all of them. And programming in general was new to most of them. I therefore introduced R very shortly, before explaining the actual meta-analytical analyses more briefly. I added some explanations and references to make the file more self-explanatory for independent study. Therefore, this post also contains materials and references on how to install R, R Studio, and some basic concepts of programming.

This blog version of my workshop can guide you through an entire meta-analysis of correlational data of a fake data frame. It covers effect size estimation, data visualization, moderator analyses for categorical and metric moderators, investigation of publication bias, sensitivity analysis, and recommendations for further reading.

Both the entire workshop and this blogpost are based on From pre-registration to publication: a non-technical primer for conducting a meta-analysis to synthesize correlational data by Daniel Quintana. You can find the original paper here and I highly recommend reading it! You can find a video of the author walking you through the entire code here.


Why R?

Reproducible Workflow

Using R enables us to easily share code between colleagues and make changes in our analysis transparent. We don´t have to simply trust the results because we can see how these results were obtained. This makes the entire scientific enterprise more trustworthy.

Most meta-analysis are NOT reproducible, to avoid this we need to adopt open data and open code practices.

For further information see recommendations by Lakens et al., 2014.

Almost anything is possible

R gets its power by its ability to load packages that enable new functionality. The way we want to do our meta-analysis would´t be possible without R!


Getting R up and running

Installing

Both R and R Studio are important and should be downloaded.

R: download via the cloud mirror which selects a good mirror automatically https://clouds.r-project.org.

R Studio: Integrated development environment (IDE) for R, http://www.rstudio.com/download.

I highly recommend this guide to become familiar with R Markdown.


R Basics

“To understand computations in R, two slogans are helpful:

  • Everything that exists is an object.

  • Everything that happens is a function call."

— John Chambers


A function call “does” something, stands before parentheses, and looks like this:

function()

function calls work on objects.

Objects can be almost anything, like numbers, matrices or entire data frames.

function(object1, object2, object3)

For example:

sum(1, 2, 3) # The function "sum()" works on the objects = numbers 1, 2, and 3
[1] 6


Missing values can produce problems, so it is important to tell R what it should do with them. Such specifications, telling R what it should do, are called arguments

mean(c(1, 2, 3, NA)) 
[1] NA
mean(c(1, 2, 3, NA), na.rm = TRUE) # the argument na.rm = TRUE tells R that all NA should be removed
[1] 2

That is it, you do not need to understand anything more on programming in order to apply all functions to conduct a meta-analysis in R!


Performing the Meta-analyses

Install and load packages

First of all, you need to install the packages needed to perform the meta-analysis.

#install.packages(c("metafor", "tidyverse", "metaviz", "puniform", "magrittr", "knitr", "readxl")) # remove the # if you have not installed those packages.

library(tidyverse)  # for data cleaning
── Attaching packages ─────────────
✓ ggplot2 3.2.1     ✓ purrr   0.3.3
✓ tibble  2.1.3     ✓ dplyr   0.8.4
✓ tidyr   1.0.2     ✓ stringr 1.4.0
✓ readr   1.3.1     ✓ forcats 0.4.0
── Conflicts ──────────────────────
x dplyr::filter() masks stats::filter()
x dplyr::lag()    masks stats::lag()
library(metafor)    # for meta-analysis
Loading required package: Matrix

Attaching package: 'Matrix'
The following objects are masked from 'package:tidyr':

    expand, pack, unpack
Loading 'metafor' package (version 2.1-0). For an overview 
and introduction to the package please type: help(metafor).
library(metaviz)    # for visualization of meta-analysis
library(puniform)   # for p-uniform analysis
Loading 'puniform' package (Version 0.2.1)
library(magrittr)   # for the pipe symbol %>%

Attaching package: 'magrittr'
The following object is masked from 'package:purrr':

    set_names
The following object is masked from 'package:tidyr':

    extract
library(knitr)      # to easily use rmarkdown
library(readxl)     # to easily read excel files


Loading data using code

You can load the data with the following code. It is very important to save the data in the exact same folder as the RMD document! Otherwise R cannot find it and won´t be able to load the data frame.

data <- read_excel("data.xlsx") 


Descriptives

To look at the data frame data you just type its name.

data
# A tibble: 20 x 9
   study_id authors    year short_ref   age    study_n r      metric_mod cat_mod
      <dbl> <chr>     <dbl> <chr>       <chr>    <dbl> <chr>  <chr>      <chr>  
 1        1 Mitchell…  2013 Mitchell e… 27.64…     649 0.304… 0.3072197… A      
 2        2 Busch et…  2002 Busch et a… 10.09…     604 0.229… 0.3289834… B      
 3        3 el-Kader…  2017 el-Kader e… 16.69…     179 0.205… 0.8773893… C      
 4        4 Garcia e…  2012 Garcia et … 26.68…     773 0.509… 0.4216361… A      
 5        5 Avalos R…  2017 Avalos Riv… 26.27…     799 0.426… 0.3038734… B      
 6        6 Benjamin…  2001 Benjamin e… 29.30…     890 0.384… 0.6374664… C      
 7        7 Jack et …  2016 Jack et al… 14.93…    1034 0.374… 0.3989912… A      
 8        8 Leon et …  2000 Leon et al… 24.83…     276 0.226… 0.9315439… B      
 9        9 Shrestha…  2009 Shrestha e… 24.91…     368 0.342… 0.3800404… C      
10       10 Smith et…  2016 Smith et a… 15.46…     864 0.391… 0.3610906… A      
11       11 Sterrett…  2011 Sterrett e… 19.29…     123 0.358… 0.5447850… B      
12       12 Heidrich…  2007 Heidrich e… 31.73…     761 0.393… 0.2687553… C      
13       13 Pham et …  2004 Pham et al… 22.26…     902 0.272… 0.5844837… A      
14       14 Meganck …  2004 Meganck et… 23.75…     753 0.522… 0.2350489… B      
15       15 Schuster…  1998 Schuster e… 18.64…    1034 0.387… 0.5282168… C      
16       16 Allen et…  2009 Allen et a… 23.38…     964 0.468… 0.3927904… A      
17       17 el-Kalil…  2017 el-Kalil e… 25.45…     787 0.360… 0.4376787… B      
18       18 al-Zia e…  2010 al-Zia et … 26.11…     275 0.272… 0.8112219… C      
19       19 Garcia e…  2007 Garcia et … 32.72…     187 0.348… 0.4103933… A      
20       20 Melton e…  2000 Melton et … 10.26…     897 0.374… 0.5642247… B      

This is the data frame we just loaded. You can see the variable names in the first row, each row is a study. We have data on the authors, year, age, sample size, correlations as well as measures on moderators like quality and group.

Another way to look at the data is glimpse():

glimpse(data)
Observations: 20
Variables: 9
$ study_id   <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17,…
$ authors    <chr> "Mitchell et al.", "Busch et al.", "el-Kader et al.", "Gar…
$ year       <dbl> 2013, 2002, 2017, 2012, 2017, 2001, 2016, 2000, 2009, 2016…
$ short_ref  <chr> "Mitchell et al., 2013", "Busch et al., 2002", "el-Kader e…
$ age        <chr> "27.64255416795089", "10.098209693014788", "16.69867223923…
$ study_n    <dbl> 649, 604, 179, 773, 799, 890, 1034, 276, 368, 864, 123, 76…
$ r          <chr> "0.30491540177964377", "0.22938155641276098", "0.205558246…
$ metric_mod <chr> "0.30721970419723954", "0.3289834977800418", "0.8773893884…
$ cat_mod    <chr> "A", "B", "C", "A", "B", "C", "A", "B", "C", "A", "B", "C"…

But something weird happened: age, r and metric_mod are for some reason characters!

This tells us something about structure of the data, if we have variables consisting of numerical values/metric (dbl), characters (chr) or factors/categorical (fct).

We need to turn them into numeric values like this:

data$age <- as.numeric(data$age)

data$r <- as.numeric(data$r)

data$metric_mod <- as.numeric(data$metric_mod)


Lets turn our categorical moderator cat_mod from a character string to a real categorical variable like this:

data$cat_mod <- as.factor(data$cat_mod)


We can summarize our data to get an overview with summary()

summary(data)
    study_id       authors               year       short_ref        
 Min.   : 1.00   Length:20          Min.   :1998   Length:20         
 1st Qu.: 5.75   Class :character   1st Qu.:2004   Class :character  
 Median :10.50   Mode  :character   Median :2009   Mode  :character  
 Mean   :10.50                      Mean   :2008                     
 3rd Qu.:15.25                      3rd Qu.:2014                     
 Max.   :20.00                      Max.   :2017                     
      age           study_n             r            metric_mod     cat_mod
 Min.   :10.10   Min.   : 123.0   Min.   :0.2056   Min.   :0.2350   A:7    
 1st Qu.:18.16   1st Qu.: 345.0   1st Qu.:0.2968   1st Qu.:0.3531   B:7    
 Median :24.30   Median : 767.0   Median :0.3674   Median :0.4160   C:6    
 Mean   :22.53   Mean   : 656.0   Mean   :0.3578   Mean   :0.4863          
 3rd Qu.:26.38   3rd Qu.: 891.8   3rd Qu.:0.3924   3rd Qu.:0.5693          
 Max.   :32.73   Max.   :1034.0   Max.   :0.5229   Max.   :0.9315          

Here you can see the range, mean, median, and missing values for each variable.


Calculate the effect size

Transform r to Z

For synthesizing correlational data, it is important to transform Pearson’s r into Fisher’s z scale, as Pearson’s r is not normally distributed and calculate the corresponding sample variances.

data_z <- escalc(measure = "ZCOR",  # function escalc() calulates z from r and sample size
                 ri = r, 
                 ni = study_n, 
                 data = data)

Let´s have a look at the new data frame data_z with two new columns, yi containing the z-transformed values as well as vi, containing the corresponding sampling variances.

data_z
   study_id              authors year                  short_ref      age 
1         1      Mitchell et al. 2013      Mitchell et al., 2013 27.64255 
2         2         Busch et al. 2002         Busch et al., 2002 10.09821 
3         3      el-Kader et al. 2017      el-Kader et al., 2017 16.69867 
4         4        Garcia et al. 2012        Garcia et al., 2012 26.68570 
5         5 Avalos Rivera et al. 2017 Avalos Rivera et al., 2017 26.27226 
6         6      Benjamin et al. 2001      Benjamin et al., 2001 29.30365 
7         7          Jack et al. 2016          Jack et al., 2016 14.93451 
8         8          Leon et al. 2000          Leon et al., 2000 24.83702 
9         9      Shrestha et al. 2009      Shrestha et al., 2009 24.91594 
10       10         Smith et al. 2016         Smith et al., 2016 15.46601 
11       11      Sterrett et al. 2011      Sterrett et al., 2011 19.29807 
12       12      Heidrich et al. 2007      Heidrich et al., 2007 31.73846 
13       13          Pham et al. 2004          Pham et al., 2004 22.26795 
14       14       Meganck et al. 2004       Meganck et al., 2004 23.75831 
15       15      Schuster et al. 1998      Schuster et al., 1998 18.64512 
16       16         Allen et al. 2009         Allen et al., 2009 23.38977 
17       17      el-Kalil et al. 2017      el-Kalil et al., 2017 25.45586 
18       18        al-Zia et al. 2010        al-Zia et al., 2010 26.11898 
19       19        Garcia et al. 2007        Garcia et al., 2007 32.72533 
20       20        Melton et al. 2000        Melton et al., 2000 10.26528 
   study_n         r metric_mod cat_mod     yi     vi 
1      649 0.3049154  0.3072197       A 0.3149 0.0015 
2      604 0.2293816  0.3289835       B 0.2335 0.0017 
3      179 0.2055582  0.8773894       C 0.2085 0.0057 
4      773 0.5097461  0.4216361       A 0.5624 0.0013 
5      799 0.4262435  0.3038734       B 0.4553 0.0013 
6      890 0.3848341  0.6374664       C 0.4057 0.0011 
7     1034 0.3743632  0.3989913       A 0.3935 0.0010 
8      276 0.2269614  0.9315440       B 0.2310 0.0037 
9      368 0.3428623  0.3800405       C 0.3573 0.0027 
10     864 0.3919133  0.3610907       A 0.4141 0.0012 
11     123 0.3580845  0.5447851       B 0.3747 0.0083 
12     761 0.3938829  0.2687553       C 0.4164 0.0013 
13     902 0.2726068  0.5844837       A 0.2797 0.0011 
14     753 0.5229400  0.2350489       B 0.5804 0.0013 
15    1034 0.3874099  0.5282169       C 0.4087 0.0010 
16     964 0.4680389  0.3927904       A 0.5076 0.0010 
17     787 0.3604838  0.4376788       B 0.3774 0.0013 
18     275 0.2722842  0.8112219       C 0.2793 0.0037 
19     187 0.3482871  0.4103933       A 0.3635 0.0054 
20     897 0.3742383  0.5642247       B 0.3933 0.0011 

Removing missing values

In case some effect sizes are missing, it is necessary to exclude these studies like this:

data_z <- data_z %>% 
  drop_na(yi)

Random-effects model

Choosing which model and which estimator you should choose would outgrow the time frame, yet in our case the random effects model is more appropriate.

rem <- rma(yi, vi,                  # use yi (effect size) and vi (sampling variance)
           data = data_z,           # use this data frame
           slab = data_z$short_ref, # use short ref for study names
           method = "REML")         # REML method for estimation
rem

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

tau^2 (estimated amount of total heterogeneity): 0.0083 (SE = 0.0033)
tau (square root of estimated tau^2 value):      0.0911
I^2 (total heterogeneity / total variability):   84.27%
H^2 (total variability / sampling variability):  6.36

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

Model Results:

estimate      se     zval    pval   ci.lb   ci.ub 
  0.3844  0.0227  16.8978  <.0001  0.3398  0.4290  *** 

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

The estimate provides the estimated model coefficient (i.e., the summary effect size) with standard error(se). In this case, the summary effect size would be 0.3843944.

The z-value is the corresponding test statistic, pval is the corresponding p-value, ci.lb the the lower bound of the confidence interval and ci.ub the upper bound of the confidence interval.

Transform z back to r

predict(rem, digits = 3, transf = transf.ztor)

  pred ci.lb ci.ub cr.lb cr.ub 
 0.367 0.327 0.404 0.198 0.514 

These two lines display the transformation of Fisher’s z back to Pearson’s r (“pred”), and the 95% CI for r (ci.lb and ci.ub) for reporting the meta-analysis.

CIs for heterogeneity (I2).

Calculate confidence interval for the amount of heterogeneity (I2).

confint(rem) 

       estimate   ci.lb   ci.ub 
tau^2    0.0083  0.0040  0.0200 
tau      0.0911  0.0630  0.1413 
I^2(%)  84.2735 71.9428 92.8046 
H^2      6.3587  3.5642 13.8978 

Report the results

metafor contains the very helpful function reporter() which produces a complete summary of the results with one line of code!

#reporter(rem) # delete the # to execute the code

Visualization

Forest Plot

For reasons outlined in Finding your way out of the forest without a trail of bread crumbs: development and evaluation of two novel displays of forest plots I prefer meta-analytic data-visualization via metaviz.

viz_forest(rem,                              # rma object from the previous step
           summary_label = "Summary effect", # how the summary effect is called
           study_labels = data_z$short_ref,  # the short reference is used for study labels
           method = "REML",                  # choose estimation method, REML recommended
           xlab = "Correlation Coefficient", # name x-axis
           variant = "rain",                 # choose visualization: rainforest, thickforest...
           annotate_CI = TRUE,               # argument, if 95%CIs should be displayed
           table_headers = "Correlation [95% CI]")

Baujat Plot

Baujat plots can be used to identify sources of effect size heterogeneity, studies in the upper right quadrant are disproportionately influential on heterogeneity and the summary effect size. For further explanation see here.

baujat(rem, symbol = "slab")

Funnel Plot

Funnel plots are great tools for visual inspection of publication bias or other sources of asymmetry. It also contains an Eggers regression line and contours. For further explanations on the interpretation of these plots, have a look at this great post by Daniel Quintana and this description of the metaviz package.

viz_funnel(x = rem, 
           contours = TRUE, # Turn on contours
           egger = TRUE)    # Turn on Eggers regression line
Note: method argument used differs from input object of class rma.uni (metafor)


Moderator analyses

Moderator analyses examine possible influences on effect size heterogeneity, meta-regressions are used for metric moderators, while subgroup analyses are used for categorical analyses.

Meta-Regression

rem_reg<- rma(yi, vi, 
              mods = ~ metric_mod,      # metric moderator of interest
              data = data_z, digits = 3) 
rem_reg

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

tau^2 (estimated amount of residual heterogeneity):     0.006 (SE = 0.003)
tau (square root of estimated tau^2 value):             0.077
I^2 (residual heterogeneity / unaccounted variability): 79.45%
H^2 (unaccounted variability / sampling variability):   4.87
R^2 (amount of heterogeneity accounted for):            28.35%

Test for Residual Heterogeneity:
QE(df = 18) = 87.460, p-val < .001

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

Model Results:

            estimate     se    zval   pval   ci.lb   ci.ub 
intrcpt        0.527  0.055   9.618  <.001   0.420   0.634  *** 
metric_mod    -0.300  0.109  -2.761  0.006  -0.513  -0.087   ** 

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Bubble Plot

bubble_mod <- ggplot(data_z, 
                     aes(metric_mod, r,      # metric moderator on x, r on y axis
                     size = 1 / sqrt(vi))) + # size of bubbles in increased by 1/sqrt
  geom_point(shape = 1) +                    # shape = 1 for hollow dots
  geom_smooth(method = lm,                   # add regression line 
              color = "grey", 
              se=T,                          # display confidence intervall
              size = 0.5, 
              alpha = .2) +  
  ylab(expression("Correlation coefficient"~italic("r"))) + # change y-axis name
  xlab("Metric Moderator") +                 # change x-axis name
  theme_classic() +                          # use theme for aesthetic reasons
  scale_size_continuous(guide = "none")      # removes legend for different point-sizes
  
bubble_mod

Subgroup Analysis

Calculate the subgroup moderator analysis like this:

rem_sub<- rma(yi, vi,           # use yi as effect size and vi as sampling variance
              mods = ~ cat_mod, # use this moderator
              data = data_z,    # use this data frame
              digits = 3)       # show only 3 digits
rem_sub

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

tau^2 (estimated amount of residual heterogeneity):     0.009 (SE = 0.004)
tau (square root of estimated tau^2 value):             0.096
I^2 (residual heterogeneity / unaccounted variability): 85.23%
H^2 (unaccounted variability / sampling variability):   6.77
R^2 (amount of heterogeneity accounted for):            0.00%

Test for Residual Heterogeneity:
QE(df = 17) = 111.105, p-val < .001

Test of Moderators (coefficients 2:3):
QM(df = 2) = 0.756, p-val = 0.685

Model Results:

          estimate     se    zval   pval   ci.lb  ci.ub 
intrcpt      0.407  0.039  10.369  <.001   0.330  0.484  *** 
cat_modB    -0.024  0.056  -0.426  0.670  -0.134  0.086      
cat_modC    -0.051  0.059  -0.869  0.385  -0.166  0.064      

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

intercept is effect size estimate for Group A, cat_modB is the difference from Group A for Group B, and cat_modC is the difference from Group A for Group C.

To obtain the effect size estimate for each group, you should calculate the difference from the intercept.

tribble(
  ~Group,     ~r,
  "Group A",   rem_sub$b[1],
  "Group B",   rem_sub$b[1] + rem_sub$b[2],
  "Group C",   rem_sub$b[1] + rem_sub$b[3]
)
# A tibble: 3 x 2
  Group       r
  <chr>   <dbl>
1 Group A 0.407
2 Group B 0.383
3 Group C 0.356
Subgroup Forest Plot

You can visualize the subgroups with a forest plot like this:

viz_forest(rem_sub,
           study_labels =  data_z$short_ref, # naming the studies
           method = "REML",                  # use reml method
           group = data$cat_mod,             # categorical moderator
           summary_label = c("Group A",      # naming groups
                             "Group B", 
                             "Group C"), 
           xlab = "Correlation",             # naming x-axis
           variant = "rain",                 # type of forest plot
           text_size = 2,                  
           annotate_CI = T)


Publication bias

There are many methods to estimate, identify, and adjust for publication bias. For a recent overview of these techniques read this great paper.

Egger´s regression test

To test for asymmetry of the effect sizes, Eggers regression can be used like this:

regtest(rem, model = "lm", predictor = "sei")

Regression Test for Funnel Plot Asymmetry

model:     weighted regression with multiplicative dispersion
predictor: standard error

test for funnel plot asymmetry: t = -1.7188, df = 18, p = 0.1028

P-Uniform

To estimate the effect size under publication bias, p-uniform analyses can be used. For interpretation have a look at the documentation.

puniform(ri = data_z$r, 
         ni = data_z$study_n, 
         alpha = .05,
         side = "right", 
         method = "P", 
         plot = TRUE)


Method: P

Effect size estimation p-uniform

       est     ci.lb     ci.ub       L.0    pval      ksig
    0.3844    0.3506    0.4107   -7.6554   <.001        20

===

Publication bias test p-uniform

      L.pb      pval
    1.2443    0.1067

===

Fixed-effect meta-analysis

    est.fe     se.fe   zval.fe pval.fe  ci.lb.fe  ci.ub.fe     Qstat   Qpval
    0.4011    0.0088     45.84   <.001     0.384    0.4183  112.7473   <.001

Sensitivity Analysis

Leave one out sensitivity analysis

By iteratively removing one study at a time and recalculating the summary effect size, studies with a large influence can be identified.

leave1out(rem, digits = 3)

                           estimate    se   zval  pval ci.lb ci.ub       Q 
Mitchell et al., 2013         0.388 0.024 16.412 0.000 0.342 0.434 107.697 
Busch et al., 2002            0.393 0.022 17.700 0.000 0.350 0.437  95.051 
el-Kader et al., 2017         0.391 0.022 17.399 0.000 0.347 0.435 106.129 
Garcia et al., 2012           0.375 0.022 17.289 0.000 0.333 0.418  91.471 
Avalos Rivera et al., 2017    0.380 0.024 16.002 0.000 0.334 0.427 110.261 
Benjamin et al., 2001         0.383 0.024 15.894 0.000 0.336 0.430 112.727 
Jack et al., 2016             0.384 0.024 15.895 0.000 0.336 0.431 112.682 
Leon et al., 2000             0.392 0.023 17.351 0.000 0.347 0.436 104.675 
Shrestha et al., 2009         0.385 0.024 16.160 0.000 0.339 0.432 112.027 
Smith et al., 2016            0.382 0.024 15.896 0.000 0.335 0.430 112.593 
Sterrett et al., 2011         0.385 0.024 16.352 0.000 0.338 0.431 112.663 
Heidrich et al., 2007         0.382 0.024 15.911 0.000 0.335 0.429 112.560 
Pham et al., 2004             0.391 0.023 16.891 0.000 0.345 0.436  98.505 
Meganck et al., 2004          0.374 0.021 17.703 0.000 0.333 0.416  87.183 
Schuster et al., 1998         0.383 0.024 15.878 0.000 0.335 0.430 112.682 
Allen et al., 2009            0.377 0.023 16.426 0.000 0.332 0.422 100.999 
el-Kalil et al., 2017         0.384 0.024 15.968 0.000 0.337 0.432 112.279 
al-Zia et al., 2010           0.389 0.023 16.762 0.000 0.344 0.435 108.626 
Garcia et al., 2007           0.385 0.024 16.273 0.000 0.339 0.431 112.483 
Melton et al., 2000           0.384 0.024 15.909 0.000 0.336 0.431 112.689 
                              Qp  tau2     I2    H2 
Mitchell et al., 2013      0.000 0.009 84.642 6.511 
Busch et al., 2002         0.000 0.007 82.605 5.749 
el-Kader et al., 2017      0.000 0.008 83.816 6.179 
Garcia et al., 2012        0.000 0.007 81.474 5.398 
Avalos Rivera et al., 2017 0.000 0.009 84.602 6.494 
Benjamin et al., 2001      0.000 0.009 84.927 6.634 
Jack et al., 2016          0.000 0.009 84.819 6.587 
Leon et al., 2000          0.000 0.008 83.703 6.136 
Shrestha et al., 2009      0.000 0.009 85.300 6.803 
Smith et al., 2016         0.000 0.009 84.917 6.630 
Sterrett et al., 2011      0.000 0.009 85.373 6.837 
Heidrich et al., 2007      0.000 0.009 85.001 6.667 
Pham et al., 2004          0.000 0.008 83.597 6.097 
Meganck et al., 2004       0.000 0.006 80.508 5.130 
Schuster et al., 1998      0.000 0.009 84.781 6.571 
Allen et al., 2009         0.000 0.008 83.293 5.986 
el-Kalil et al., 2017      0.000 0.009 85.035 6.682 
al-Zia et al., 2010        0.000 0.008 84.623 6.503 
Garcia et al., 2007        0.000 0.009 85.387 6.843 
Melton et al., 2000        0.000 0.009 84.947 6.643 
dat_leave1out <- leave1out(rem, digits = 3)

str(dat_leave1out)
List of 14
 $ estimate: num [1:20] 0.388 0.393 0.391 0.375 0.38 ...
 $ se      : num [1:20] 0.0236 0.0222 0.0225 0.0217 0.0238 ...
 $ zval    : num [1:20] 16.4 17.7 17.4 17.3 16 ...
 $ pval    : num [1:20] 1.56e-60 4.19e-70 8.39e-68 5.71e-67 1.23e-57 ...
 $ ci.lb   : num [1:20] 0.342 0.35 0.347 0.333 0.334 ...
 $ ci.ub   : num [1:20] 0.434 0.437 0.435 0.418 0.427 ...
 $ Q       : num [1:20] 107.7 95.1 106.1 91.5 110.3 ...
 $ Qp      : num [1:20] 8.43e-15 1.77e-12 1.65e-14 7.85e-12 2.81e-15 ...
 $ tau2    : num [1:20] 0.00854 0.00733 0.00771 0.00688 0.00862 ...
 $ I2      : num [1:20] 84.6 82.6 83.8 81.5 84.6 ...
 $ H2      : num [1:20] 6.51 5.75 6.18 5.4 6.49 ...
 $ slab    : chr [1:20] "Mitchell et al., 2013" "Busch et al., 2002" "el-Kader et al., 2017" "Garcia et al., 2012" ...
 $ digits  : Named num [1:9] 3 3 3 3 3 3 3 3 3
  ..- attr(*, "names")= chr [1:9] "est" "se" "test" "pval" ...
 $ transf  : logi FALSE
 - attr(*, "class")= chr "list.rma"

In our case, the effect sizes range from 0.37 to 0.39. Which would indicate no extreme influence of single studies.

Influence analysis

Additional influence analyse can be conducted like this:

inf_neuro_x_ad<- influence(rem)

plot(inf_neuro_x_ad, layout = c(4, 2))


Further Reading and Resources

Daniel Quintana, whose above mentioned paper inspired this workshop, is also is the co-host of the Everything hertz podcast, in which he discusses and interviews people on open science, meta-research, meta analysis, and related topics. I can especially recommend this episode in which he defends the idea of meta-analysis to his cohost James Heathers.

The creator of the package metafor, Wolfgang Viechtbauer walks you through multiple examples for meta-analyses on his website. Most meta-analytic analyses are described there. In particular, I used this example in this file.

Another helpful resource is the book Doing Meta-Analysis in R which covers more theory without being impractical.

If you want to dive further into programming with R or need to look something up, I can highly recommend this free introductory book. It covers all the basics on how to clean, visualize, and report data analyses.

Also, printing out cheat sheets and keeping them next to your computer is a great way to learn R, for example for the basics.

I can also recommend Datacamp courses which cover most topics related to programming in R.

On a site note, the most important resource for troubleshooting is Stackoverflow. If you are facing a problem, chances are very high that somebody else had that same problem and its solution can be found there!