Learning outcomes

This tutorial guides you through an example of a meta-analysis (using the method of Hunter and Schmidt) conducted in R statistical programming environment. By the end of the tutorial, a you should be able to:

  • Convert effect sizes reported in publications to a common, standardized effect size
  • Summarise the basic study information and standardized effect sizes in a table
  • Perform a barebones meta-analysis (Hunter and Schmidt method) with heterogeneity analysis
  • Summarize the results of the meta-analysis as a table and forest plot
  • Assess publication bias using funnel plots and cumulative meta-analysis
  • Assess whether some studies could be influential outliers using leave-one-out meta-analysis
  • Interpret and report the analysis results

Installing R

Please follow the link below for a refresher on installing and getting started with R. It includes an overview of R, and a tour of a graphical user interface for R, called RStudio. Some of you may find RStudio helpful but it isn’t necessary to conduct your statistical analysis in R. You might want to skip the sections on ‘R Markdown’, ‘Getting data into RStudio’ and ‘Tables in R Markdown’.

http://milton-the-cat.rocks/learnr/r/r_getting_started/

Note that in this document we use some command line syntax for stringing functions together more efficiently. This syntax is part of a collection of packages called tidyverse and a key feature of it is pipe operators such as %>%. It is used to send the output of one function directly to the input of the next function. The following is a useful introduction to this functionality of tidyverse:

https://www.datacamp.com/community/tutorials/pipe-r-tutorial

# R is open-source software whose functionality can be extended by installing extra packages often 
# written by academics for the research community to use. The following commands are necessary only 
# once to install packages for the first time. We will install some extra packages to conduct the 
# example meta-analysis into an existing R installation.
if (!require(package="psychmeta")) utils::install.packages("psychmeta", dependencies=TRUE)  
if (!require(package="effectsize")) utils::install.packages("effectsize", dependencies=TRUE)    
if (!require(package="esc")) utils::install.packages("esc", dependencies=TRUE)  
if (!require(package="ggplot2")) utils::install.packages("esc", dependencies=TRUE)  
if (!require(package="tidyverse")) utils::install.packages("tidyverse", dependencies=TRUE)  
if (!require(package="knitr")) utils::install.packages("knitr")
if (!require(package="kableExtra")) utils::install.packages("kableExtra")
if (!require(package="ggrepel")) utils::install.packages("ggrepel")
if (!require(package="devEMF")) utils::install.packages("devEMF")

Loading packages

# Whenever you start/restart R you will need to load the above packages. You can do this with the 
# following lines of code.
library(psychmeta)  
library(effectsize) 
library(esc)
library(ggplot2)
library(tidyverse)
library(knitr)
library(kableExtra)
library(ggrepel)
library(devEMF)

Calculating effect sizes

To perform a meta-analysis, one needs to convert the quantitative data reported in publications to a common statistic, i.e. standardize how the data is represented. Our meta-analysis will focus on a meta-analysis on effect sizes (i.e. statistics that quantify the size of the effect being studied), since this is of primary interest in scientific research. For a good primer on effect sizes I direct you to the following, brief article and simply written document:

http://www.discoveringstatistics.com/docs/psypageffectsize.pdf

There are many standardized effect sizes and many of these may be unfamiliar to you (e.g. Cohen’s d, Odd’s ratio). The fact that these so called standardized effect sizes may be unfamilar to many students and researchers does sort of defeat the object of standardizing! But there is a standardized effect size that is worth it’s salt and that is Pearson’s correlation coefficient (r). The correlation coefficient is the slope of the line relating one variable (e.g. y) to another (e.g. x) in units of their respective standard deviations. In other words, r can be thought of as the slope of the relationship between standardized variables. For example, if r = .50, then, for each increase of 1 SD in x, there is an increase of .50 SD in y. The square of the correlation coefficient, r̅² (known as the coefficient of determination), has an intuitive interpretation as the % of variance explained by the effect of the independent variable (e.g. as opposed to sampling error). Cohen provided widely adoped recommendations as to what constitutes small or large effect sizes: An |r| <.10 is considered a small effect (explaining 0 - 1% of the total variance), an |r| of .10 - .30 is a medium effect (accounting for 1 - 9% of the total variance) and an |r| >.50 (large effect, accounting for >25% of the variance). As well as being a statistic familiar to most students with basic maths training, r can be calculated (and interpreted in) many different ways making it an excellent choice for a standardized effect size (Rodgers and Nicewander, 1988). For a meta-analysis, we will also need the sample size of the experiment (n) used to estimate r in each study so that our uncertainty of r can be evaluated.

Depending on the research area in which you are conducting the meta-analysis, one of the hardest parts of the meta-analysis can be standardizing the statistics from different publications. Below we will illustrate how to do some effect size conversions in R. It is important to select the data or statistics that test the common hypothesis of interest across the studies and that are relevant to your research question. There maybe some differences in the experiments but the influence of other variables can be evaluated as part of your meta-analysis.

N.B. The meta-analysis method we will use is a random effects model proposed by Hunter-Schmidt for the analysis of correlation coefficients. Note that if you insist on conducting a meta-analysis using an alternative standardized effect size (rather than you converting them to correlation coefficients) you will need to use a different R package etc.

# Convert two (independent) sample t-test result to r (using the t-statistic and df) 
t <- 4.56; df <- 48; 
convert_es(es=t, input_es="t", output_es="r", df1=df)
## r values converted from t values
## -----------------------------------------
##      r n_effective  n n1 n2   var_e
## 1 0.55          50 50 NA NA 0.00994
# Convert two (independent) sample t-test result to r (using the t-statistic and  sample sizes, N1 and N2)      
t <- 4.96; N1 <- 38; N2 <- 38; 
convert_es(es=t, input_es="t", output_es="r", n1=N1, n2=N2)
## r values converted from t values
## -----------------------------------------
##     r n_effective  n n1 n2   var_e
## 1 0.5          76 76 38 38 0.00751
# Convert two (independent) sample t-test result to r (using the p-value and sample sizes, N1 and N2)       
p <- 0.085; N1 <- 46; N2 <- 47; 
convert_es(es=p, input_es="p.t", output_es="r", n1=N1, n2=N2, tails=2)
## r values converted from p.t values
## -----------------------------------------
##      r n_effective  n n1 n2  var_e
## 1 0.18          93 93 46 47 0.0102
# Convert summary statistics of two independent samples to r (using the M, and SEM and N)       
# N is simply ths sum of N1 and N2: N = N1 + N2
# In this example, 20 + 23 = 43
M1 <- 20.1; M2 <- 5.13; SEM1 <- 0.9; SEM2 <- 1.1; N1 <- 20; N2 <- 23;   
esc_mean_se(M1, SEM1, N1, M2, SEM2, N2, es.type="r")
## 
## Effect Size Calculation for Meta Analysis
## 
##      Conversion: mean and se to effect size correlation
##     Effect Size:   0.8500
##  Standard Error:   0.2254
##        Variance:   0.0508
##        Lower CI:   0.6720
##        Upper CI:   0.9351
##          Weight:  19.6821
##      Fisher's z:   1.2561
##       Lower CIz:   0.8143
##       Upper CIz:   1.6979
# Convert summary statistics of two independent samples to r (using the M, SD and N)        
# N is simply ths sum of N1 and N2: N = N1 + N2
# In this example, 15 + 18 = 33
M1 <- 16.01; M2 <- 9.45; SD1 <- 4.5; SD2 <- 3.8; N1 <- 15; N2 <- 18; 
esc_mean_sd(M1, SD1, N1, M2, SD2, N2, es.type = "r")
## 
## Effect Size Calculation for Meta Analysis
## 
##      Conversion: mean and sd to effect size correlation
##     Effect Size:   0.6203
##  Standard Error:   0.1956
##        Variance:   0.0383
##        Lower CI:   0.3293
##        Upper CI:   0.8036
##          Weight:  26.1403
##      Fisher's z:   0.7254
##       Lower CIz:   0.3421
##       Upper CIz:   1.1088
# Convert one-way ANOVA test result on two samples to r (using the F-statistic and sample sizes, N1 and N2)     
F <- 34.4; N1 <- 35; N2 <- 35; 
convert_es(es=F, input_es="F", output_es="r", n1=N1, n2=N2)
## r values converted from F values
## -----------------------------------------
##      r n_effective  n n1 n2   var_e
## 1 0.58          70 70 35 35 0.00639
# Convert paired t-test result to r (using the t-statistic and df)
# N is simply 1 more than the degrees of freedom: N = df + 1 
# In this example, 44 + 1 = 45
t <- 6.7; df <- 44; 
t_to_r(t, df, Paired=TRUE)
## r    |       95% CI
## -------------------
## 0.71 | [0.54, 0.81]
# Convert chi-squared test result to r (using the c2 statistic and sample sizes, N1 and N2)     
chisq <- 27.5; N1 <- 31; N2 <- 34; 
convert_es(es=chisq, input_es="chisq", output_es="r", n1=N1, n2=N2)
## r values converted from chisq values
## -----------------------------------------
##      r n_effective  n n1 n2  var_e
## 1 0.65          65 65 31 34 0.0052
# Convert Cohen's d to r    
d <- 1.5; N1 <- 20; N2 <- 21; 
convert_es(es=d, input_es="d", output_es="r", n1=N1, n2=N2)
## r values converted from d values
## -----------------------------------------
##     r n_effective  n n1 n2  var_e
## 1 0.6          41 41 20 21 0.0102
# Convert Odds ratio to r       
OR <- 43.4; N1 <- 22; N2 <- 25; 
convert_es(es=OR, input_es="or", output_es="r", n1=N1, n2=N2)
## r values converted from or values
## -----------------------------------------
##      r n_effective  n n1 n2   var_e
## 1 0.72          47 47 22 25 0.00504
# Convert Log Odds ratio to r       
LOR <- 3.77; N1 <- 22; N2 <- 25; 
convert_es(es=LOR, input_es="lor", output_es="r", n1=N1, n2=N2)
## r values converted from lor values
## -----------------------------------------
##      r n_effective  n n1 n2   var_e
## 1 0.72          47 47 22 25 0.00505

Barebones meta-analysis

The following commands demonstrate how to conduct the meta-analysis from a set of standardized effect sizes. For a good primer on meta-analysis I direct you to the following, brief article and simply written document:

https://www.discoveringstatistics.com/docs/meta.pdf

# We illustrated above how one can convert various statistics to correlation coefficients and total sample sizes
# Having collected that data now input that data so that we can perform the meta-analysis in R 
# Assign vector of effect sizes to a new variable ri
ri <- c(.55, .50, .18, .85, .62, .71, .58, .65, .60, .72)
# Assign vector of total sample sizes to a new variable ni  
ni <- c(50, 76, 93, 43, 33, 45, 70, 65, 41, 47)
# Automatically create vector of character strings to use as study IDs
study_id <- format(1:length(ni), digits=1, trim=FALSE) # Study IDs (could be replaced with vector)
# Create a vector of study Author
study <- c("Barnett (2015)",
           "Barnett et al. (2010)",
           "Dodds et al. (2017)",
           "Schroeder and Kennel (2018)",
           "Heyward et al. (2009)",
           "Kennell (2012)",
           "Kennell et al. (2014)",
           "Sharp et al. (2011)",
           "Silver et al. (2016)",
           "Spencer et al. (2008)");
# Create nice table summarizing standardized effect sizes for the studies that will 
# be used in the meta-analysis
tibble(study_id = study_id, study = study, ri = ri, ni = ni) %>%
  knitr::kable(caption = "Summary of study effect sizes", digits = 2,
               col.names = c("Study ID", "Author, Date", "*r*", "*n*")) %>%
  kableExtra::kable_styling(bootstrap_options = c("striped", "hover"))
Summary of study effect sizes
Study ID Author, Date r n
1 Barnett (2015) 0.55 50
2 Barnett et al. (2010) 0.50 76
3 Dodds et al. (2017) 0.18 93
4 Schroeder and Kennel (2018) 0.85 43
5 Heyward et al. (2009) 0.62 33
6 Kennell (2012) 0.71 45
7 Kennell et al. (2014) 0.58 70
8 Sharp et al. (2011) 0.65 65
9 Silver et al. (2016) 0.60 41
10 Spencer et al. (2008) 0.72 47 
# Run the barebones meta-analysis using the method of Hunter & Schmidt
result <- ma_r(ri, ni, ma_method="bb", 
               control=control_psychmeta(cred_level=0.95),
               sample_id=study_id)
##  **** Running ma_r: Meta-analysis of correlations ****
# Preview the summary of results
summary(result)
## Bare-bones meta-analysis results 
## ---------------------------------------------------------------------- 
##   analysis_id pair_id construct_x construct_y analysis_type  k   N mean_r  sd_r
## 1           1       1           X           Y       Overall 10 563   0.56 0.203
##     se_r sd_res CI_LL_95 CI_UL_95 CR_LL_95 CR_UL_95
## 1 0.0641   0.18    0.415    0.705    0.152    0.969
## 
## 
## Information available in the meta-analysis object includes:
##  - meta_tables   [ access using get_metatab() ]
##  - escalc        [ access using get_escalc() ]
# Perform heterogeneity analysis and append the output to results
result <- heterogeneity(result) 
# Preview the results of the heterogenity analysis
result$heterogeneity
## $`analysis id: 1`
## $`analysis id: 1`$barebones
## 
## Heterogeneity results for r 
## ---------------------------
## 
## Accounted for a total of 20.754% of variance
## 
## Correlation between r values and artifactual perturbations: 0.456
## 
## The reliability of observed effect sizes is: 0.792
## 
## 
## Random effects variance estimates
## ---------------------------------
## Hunter-Schmidt method (with k-correction):
##   sd_res  (tau):   0.180, SE = 0.001, 95% CI = [0.104, 0.358] 
##   var_res (tau^2): 0.033, SE = 0.000, 95% CI = [0.011, 0.128] 
## 
##   Q statistic: 43.365 (df = 9, p = 0.000) 
##   H: 2.195   H^2: 4.818   I^2: 79.246
## 
## DerSimonian-Laird method:
##   sd_res  (tau):   0.182
##   var_res (tau^2): 0.033
## 
##   Q statistic: 43.577
##   H: 2.200   H^2: 4.842   I^2: 79.347
## 
## Outlier-robust method (absolute deviation from mean):
##   sd_res  (tau_r):   0.152
##   var_res (tau_r^2): 0.023
## 
##   Q_r statistic: 14.484
##   H_r: 1.914   H_r^2: 3.662   I_r^2: 0.727
## 
## Outlier-robust method (absolute deviation from median):
##   sd_res  (tau_m):   0.130
##   var_res (tau_m^2): 0.017
## 
##   Q_m statistic: 13.688
##   H_m: 1.716   H_m^2: 2.943   I_m^2: 0.660
## 
## 
## $`analysis id: 1`$individual_correction
## NULL
## 
## $`analysis id: 1`$artifact_distribution
## NULL

N.B. The “Accounted for a total of % of variance” refers to % variance accounted for by sampling error. This is the opposite of what I² tells us.

# For convenience, first assign the heterogeneity results to a new variable
het_result <- result$heterogeneity[[1]][[1]]$HS_method

# Set whether to truncate intervals that exceed +/- 1.0 
#
# Correlation coefficients can only adopt values in the range -1.0 to +1.0. 
# For this reason, many methods use Fisher's z transformation before averaging 
# correlation coefficients or calculating variance (or other derived statistics 
# such as confidence intervals). However, Hunter and Schmidt make arguments for 
# performing such calculations on the original scale of the correlation coefficient.
# This can produce intervals with impossible limits (i.e. that exceed +/- 1.0) 
# since the distribution of the correlation coefficient becomes skewed towards 
# the extreme ends. While this may seem unacceptable, absolute values of the 
# correlation coefficient less than 0.5 can be normally distributed and it's 
# this end of confidence or credibility intervals that is often more of interest. 
# The end of the interval that exceeds plausible values will be very close to 
# +/- 1 and so it is reasonable to truncate the interval in respect of the limits.
truncate_intervals = TRUE
if (truncate_intervals==TRUE) {
  lim = 1.0
}  else {
  lim = Inf
}

# Now let's summarise all the results in a nice table!
result %>%
  get_metatab() %>%
  as_tibble() %>%
  mutate(CI_95 = sprintf("[%.2f, %.2f]", max(CI_LL_95,-lim), min(CI_UL_95,lim))) %>%
  mutate(I2 = max(0,het_result$I_squared)) %>%
  mutate(CR_95 = sprintf("[%.2f, %.2f]", max(CR_LL_95,-lim), min(CR_UL_95,lim))) %>%  
  dplyr::select(c(k, N, mean_r, se_r, CI_95, I2, CR_95)) %>%
  knitr::kable(caption = "Summary of meta-analysis results", digits = 2, align = "c",
               col.names = c("*k*", "*N*", "*r\u0305*", "*SE*", "95% *CI*", sprintf("*I*\u00b2"), "95% *CR*")) %>%
  kableExtra::kable_styling(bootstrap_options = c("striped", "hover"), row_label_position="c")
Summary of meta-analysis results
k N r̅ SE 95% CI I² 95% CR
10 563 0.56 0.06 [0.42, 0.71] 79.25 [0.15, 0.97]

From left to right:

  • k = Number of studies in the meta-analysis
  • N = Total sample size of the meta-analysis
  • r̅ = The weighted mean correlation coefficient (effect size)
  • SE = Standard error of the correlation coefficient
  • 95% CI = 95% confidence interval
  • I² is the % variance due to inconsistencies between the population effect across the studies
  • 95% CR = 95% credibility interval

In the above example, you can see that we meta-analysed 10 studies wth a total sample size of 563. The weighted mean effect size was large (r̅ = 0.56) and indicated that the effect of interest accounted for 32% of the total variance (r̅² = 0.32). The overall effect size could be estimated with high enough precision (SE = 0.06) to exclude zero effect in 95% CI [0.42, 0.71]. This indicates that if this meta-analysis (and all the studies whose data on which it depends) was to be repeated many times, we would expect the average effect r̅ to be greater than .42 (18% of the variance) and less than .72 (50% of the variance) 95% of the times. Based on those statistics, we feel confident that there is a medium-to-large effect of…[what you would write here depends on the aim of your project!]. However, the average effect size may not be quantifying a common effect across the set of studies in this meta-analysis. To examine whether the reported study effects are estimates of a single, common underlying effect, we calculated statistics that quantify heterogeneity between population effect sizes. The I² value was high (>50%) and showed that 79% of the total variance in the reported effects was due to inconsistencies in the underlying population effect across the studies. We quantified the uncertainty about the population effect sizes on the scale of the effect size by calculating credibility intervals (CR). From these, we estimate that 95% of the population effect sizes lie in the interval [0.15, 0.97]. This suggests that some other factor or factors, which differ between the studies, are having a pronounced effect on the reported effects of… We will explore this idea further below.

The analysis method used here is the bare-bones analysis proposed by Hunter and Schmidt. The original method is modified here to improve small sample performance and reduce bias, as is done in the function ma_r in the R package psychmeta. For more information on the calculations and how to interpret the statistics from this analsis method, please see Field and Gillet (2010).

# Create forest plot to summarize the effect sizes of all the studies
result <- plot_forest(result, conf_level = 0.95, y_lab="Study ID")

# Get forest plot handle
h1 <- get_plots(result)[["forest"]][[1]][["unmoderated"]][["barebones"]]

# Edit the default graph formatting
# The size of the circle is proportional to the size of the study
circle_size   <- ni*0.1 # 0.1 is just an arbitrary scale factor to adjust the relative size of the points (observed effect). Adjust this as you see fit
triangle_size <- 6     # 10 is just a reasonable size for the triangle on the graph (overall effect). Adjust this as you see fit
h1$layers[[1]]$aes_params$size <- c(circle_size, triangle_size) 
h1$layers[[1]]$aes_params$colour <- "red" # Set marker colour to red
h1$layers <- h1$layers[4:1]               # reverse order of layers so that the markers are on top

# Display forest plot after reformatting using ggplot
# Useful cheat sheet for ggplot2: https://github.com/rstudio/cheatsheets/raw/master/data-visualization.pdf
h1 + theme(text = element_text(size=16)) -> h1  # Set font size
emf("forest.emf", width=5.5, height=3.5)  # Save plot as vector graphic (enhanced meta file)
h1 # Send graph to file
dev.off() 
h1 # Display graph in here knitted document

This forest plot summarizes the effect sizes (correlation coefficient) of all the studies from the meta-analysis in a single graph along with their 95% confidence intervals. Each circular marker point indicates the observed effect size (correlation coefficient) plotted against the study ID assigned to the original data. The size of the marker is proportional to the sample size of that particular study. The triangular marker point indicates the observed effect size (r̅) calculated in the meta-analysis. The dashed vertical line signifies zero effect.

In the above example, we can see that in most studies (except study 3), there is some consistency in the effect size being reported and the overall effect size is large (r̅ > 0.5). Interestingly, the study that reported the smallest effect (study 3) had the largest sample size.

Assess publication bias

# Create funnel plot
result <- plot_funnel(result)

# Get funnel plot handle
h2 <- get_plots(result)[["funnel"]][[1]][["barebones"]]

# Change the size of the points on the graph.
# The size of the circle is proportional to the size of the study
circle_size   <- ni*0.1   # 0.1 is just an arbitrary scale factor to adjust the relative size of the points
h2$layers[[1]]$aes_params$size <- circle_size
h2$layers[[1]]$aes_params$alpha <- 0.5         # Set alpha to 0.5 (to enhance transparency)
h2$layers[[1]]$aes_params$colour <- "red"      # Set marker colour to red

# Display the funnel plot after reformatting using ggplot
h2 + geom_label_repel(aes_(y = ~yi),label=study_id) +
    theme(text = element_text(size=16)) -> h2
emf("funnel.emf", width=5.5, height=3.5)  # Save plot as vector graphic (enhanced meta file)
h2 # Send graph to file
dev.off() 
h2 # Display graph in here knitted document

The funnel plot (shown above) is a graphical method to detect bias in the set of studies included in the meta-analysis (e.g. reporting/publication bias or poor methodological design). Each marker point indicates the observed effect size (r) plotted against the respective standard error for the effect size observed in that study. The size of the marker is proportional to the sample size of that particular study and the label corresponds to the study ID assigned to the original data. The solid, vertical line represents the overall effect size (r̅) calculated in the meta-analysis. The dashed diagonal lines indicate the limits of the 95% confidence region, and the dotted line indicates the limits of the the 99% confidence region. Study effect sizes falling outside these regions usually represent unexpected deviation from the norm and warrant further investigation. However, you can expect 1/20 studies to fall outside the 95% confidence region, and 1/100 studies to fall outside the 99% confidence region and so this must be considered when interpreting the funnel plot. An unbiased sample would ideally show a cloud of data points that is symmetric around the population effect size and has the shape of a funnel. Note that funnel plots are not particularly reliable if the number of studies in the meta-analysis is fewer than 10 and interpreting them is subjective.

In the above example, the points in the above plot are not particularly symmetric in their distribution around the overall effect size estimate (2 or 3 to the left 7 to the right), with the studies tending to overestimate the effect size. 2 out of 10 (20%) of these studies lie outside the 99% confidence region, well above the frequency we would expect (1%). Together, this suggests that small effect sizes are being observed but are not being reported/published (or at the least, are not included in this meta-analysis).

# Perform sensitivity analysis
result <- sensitivity(result,sort_method = c("inv_var"))

# Edit and display forest plot of cumulative meta-analysis results
h3 <- get_cumulative(result)[[1]][["barebones"]][["plots"]]$mean_plot
h3 + geom_count(size=4, colour="red") +
    xlab("Study ID (sorted by effect size precision)") +
    theme(text = element_text(size=16)) -> h3
emf("cumulative.emf", width=5.5, height=3.5)
h3 # Send graph to file
dev.off() 
h3 # Display graph in here knitted document

The above shows another forest plot this time summarising the results of a cumulative meta-analysis. In a cumulative meta-analysis, we investigate the effect of repeating the analysis procedure, each time adding a new study. The order in which studies are added depends on the objective of the analysis. Here, our interest is in publication bias so we have studies sorted from most to least precise (from top to bottom). The axis on the left indicates which study was further added to the analysis. The overall effect size (r̅) estimated for each analysis is plot as a red circle. The intervals plotted represent the corresponding 95% CI and CR. The solid vertical line represents the overall effect size calculated when all studies are included in the analysis. The dashed vertical line demarcates zero effect. You might suspect publication bias if the effects in the most precise studies (at the top) are small but increase as the less precise studies are added (going down) - this is known as the ‘small-study effect’. However, be aware that publication bias is not the only reason why one might see the small-study effect (e.g. smaller studies might have controlled better for confounding variables and nuisance variability).

The example above appears to show the ‘small-study effect’. Together with the result of the funnel plot, we are inclined to suggest that the overall effect size estimated from our meta-analysis may be artificially inflated due to some publication bias.

Influence of outliers

# Edit and display forest plot of leave-one-out meta-analysis results
h4 <- get_leave1out(result)[[1]]$barebones$plots$mean_plot
h4 + geom_count(size=4, colour="red") +
    xlab("Study ID") +
    theme(text = element_text(size=16)) -> h4
emf("leave1out.emf", width=5.5, height=3.5)  # Save plot as vector graphic (enhanced meta file)
h4 # Send graph to file
dev.off() 
h4 # Display graph in here knitted document

The above shows another forest plot this time summarising the results of leave-one-out meta-analysis. In this procedure, we investigate the effect of repeating the meta-analysis, each time with a different study omitted. (Thus unlike the cumulative meta-analysis, each iteration has k-1 studies, where k is the total number of studies in the original meta analysis). The axis on the left indicates which study was omitted from that analysis. The overall effect size (r̅) estimated for each analysis is plot as a red circle. The intervals plotted represent the corresponding 95% CI and CR. The solid vertical line represents the overall effect size calculated when all studies are included in the analysis. The dashed vertical line demarcates zero effect. Leave-one-out meta-analysis is a useful procedure to evaluate whether any particular study is an influential outlier. If a study is an influential outlier, we can expect the overall effect size and it’s intervals from the meta-analysis to change (more than the other leave-one-out meta-analyses).

Casting our thoughts back to the funnel plot, study 3 had a much smaller effect size well beyond the lower limit of the 99% confidence region. While it is clearly a different effect from the other studies, the extent to which it affects our effect size estimate (r̅) is unclear. In the forest plot of our leave-one-out meta-analysis, we can clearly see that omission of study 3 has a huge impact, not only on the magnitude of r̅ but also on the width of the credibility intervals (CR). This suggests that much of the heterogeneity in population effect sizes arises from study 3. It is worth carefully examining the methods of the studies to see if a reason can be identified. However, given our assessment of publication bias discussed earlier, we need to bear in mind that omitting study 3 could just exacerbate the suspected publication bias.

References

  • Brannick, M.T. & Hall, S.M. (2001) Reducing Bias in the Schmidt-Hunter Meta-analysis. Poster presented at the 16th annual conference of the Society for Industrial and Organizational Psychology, San Diego, April 2001
  • Dahlke J, Wiernik B (2019). “psychmeta: An R Package for Psychometric Meta-Analysis.” Applied Psychological Measurement, 43(5), 415–416.
  • Field (2001) Meta-analysis of correlation coefficients: A Monte Carlo comparison of fixed- and random-effects methods. Psychol Methods. 6(2): 161-180
  • Field (2006) A bluffer’s guide to meta-analysis. PsyPAG Quarterly Issue 60
    https://www.discoveringstatistics.com/docs/meta.pdf
  • Field, A.P. & Wright, D. (2006) A bluffer’s guide to effect sizes. PsyPAG Quarterly Issue 58: 9-23
    http://www.discoveringstatistics.com/docs/psypageffectsize.pdf
  • Field, A.P. & Gillett, R. (2010) How to do a meta-analysis. British journal of mathematical and statistical psychology. Vol.63(3), p.665-694
  • Rodgers, J.L. and Nicewander, W.A. (1988) Thirteen ways to look at the correlation coefficient. The American Statistician. 42(1): 59-66
  • Schmidt, F. L., & Hunter, J. E. (2015). Methods of meta-analysis: Correcting error and bias in research findings (3rd ed.). Sage. doi: 10.4135/9781483398105. Chapter 4.
  • Wiernik, B.M. & Dahlke, J.A. (2020) Obtaining Unbiased Results in Meta-Analysis: The Importance of Correcting for Statistical Artifacts. Advances in Methods and Practices in Psychological Science. 3(1): 94-123
  • Wiernik, B. & Dahlke, J. (2019) psychmeta: An Introductory Tutorial.
    https://cran.r-project.org/web/packages/psychmeta/vignettes/overview.html
  • Wiernik, B. (2021) psychmeta: Meta-analyzing correlations.
    https://cran.r-project.org/web/packages/psychmeta/vignettes/ma_r.html