This code through explores how to conduct a meta-analysis using the robumeta package in R (Fisher & Tipton, 2015).
Specifically, we’ll explain and demonstrate how to aggregate data from multiple studies to calculate the average effect of an intervention on an outcome. We will also show how meta-regression can be used to examine how the effect of an intervention varies based on other study factors, as well as how to create a forest plot using forest.robu (Tipton, 2017).
This topic is valuable because aggregating data from multiple studies creates a larger sample size, and therefore, increases the reliability of estimates of treatment effects. With the help of this code through, the robumeta package in R is a free and accessible way to conduct a meta-analysis.
Specifically, you’ll learn how to 1) calculate an average effect across multiple studies; 2) calculate average effects for subsets of data across multiple studies; 3) use meta-regression to assess how effect sizes vary based on other study factors; and 4) create a forest plot using forest.robu.
Here, we’ll show how to conduct a meta-analysis using a dataset with information from 21 studies that evaluated the impact of buprenorphine treatment for opioid use disorder on quality of life (QoL). This dataset is a simplified version of the dataset from my study that was recently published in Drug and Alcohol Dependence (Golan et al., 2022).
A basic example shows how to calculate an average effect across multiple studies. In this meta-analysis, a study could report outcomes for 1-5 QoL outcome measures (overall QoL, physical QoL, psychological QoL, social QoL, and/or environmental QoL). The data includes an effect size (variable name: “es”) and variance (variable name: “var”) for each outcome measure for each study.
The estimate of the average effect of buprenorphine on QoL is 0.757 (p<.001). Because this number is greater than 0, we can say that there is a statistically significant, positive effect of buprenorphine on QoL.
# Read in the data file.
data <- read.csv("~/Desktop/qolbupdata.csv")
# Preview data file.
head(data)# Calculate the overall mean for all QoL domains.
overallmean <- robu(formula = es ~ 1, # es (effect size) is the outcome
# measure in a meta-analysis
data = data,
studynum = author, # studynum tells R to aggregate by study
var.eff.size = var, # var.eff.size associates a variance with
# each effect size
rho = 0.8, # rho assumes a correlation of 0.8 among
# effect sizes
small = TRUE) # small = TRUE adjusts estimates
# (increases std. errors) due to small
# sample sizes
print(overallmean)## RVE: Correlated Effects Model with Small-Sample Corrections
##
## Model: es ~ 1
##
## Number of studies = 21
## Number of outcomes = 59 (min = 1 , mean = 2.81 , median = 3 , max = 5 )
## Rho = 0.8
## I.sq = 96.15286
## Tau.sq = 0.1697477
##
## Estimate StdErr t-value dfs P(|t|>) 95% CI.L 95% CI.U Sig
## 1 X.Intercept. 0.757 0.154 4.9 19.1 0.0000963 0.434 1.08 ***
## ---
## Signif. codes: < .01 *** < .05 ** < .10 *
## ---
## Note: If df < 4, do not trust the results
We may also be interested in the average effect sizes for different domains of QoL (overall QoL, physical QoL, psychological QoL, social QoL, and/or environmental QoL). The estimate is the mean domain effect size controlling for all other domains. Notice that the domain with the largest average effect is “overall QoL” (est.: 0.930), and the domain with the smallest average effect is “environmental QoL” (est.: 0.438). However, for all 5 domains of QoL, buprenorphine had a statistically significant, positive effect.
domainmeans <- robu(formula = es ~ domaincategory - 1, # es ~ domaincategory - 1
# produces separate
# estimates for each
# domain
data = data,
studynum = author,
var.eff.size = var,
rho = 0.8,
small = TRUE)
print(domainmeans)## RVE: Correlated Effects Model with Small-Sample Corrections
##
## Model: es ~ domaincategory - 1
##
## Number of studies = 21
## Number of outcomes = 59 (min = 1 , mean = 2.81 , median = 3 , max = 5 )
## Rho = 0.8
## I.sq = 96.14179
## Tau.sq = 0.1737823
##
## Estimate StdErr t-value dfs P(|t|>) 95% CI.L
## 1 domaincategoryEnvironmental 0.438 0.131 3.35 5.60 0.01709 0.113
## 2 domaincategoryOverall 0.930 0.309 3.00 6.71 0.02084 0.192
## 3 domaincategoryPhysical 0.642 0.181 3.54 11.69 0.00420 0.246
## 4 domaincategoryPsychological 0.796 0.184 4.33 12.62 0.00088 0.397
## 5 domaincategorySocial 0.680 0.209 3.25 10.32 0.00836 0.216
## 95% CI.U Sig
## 1 0.763 **
## 2 1.668 **
## 3 1.037 ***
## 4 1.195 ***
## 5 1.144 ***
## ---
## Signif. codes: < .01 *** < .05 ** < .10 *
## ---
## Note: If df < 4, do not trust the results
What’s more, it can also be used to assess how effect sizes vary based on other study factors. This is called meta-regression. Below, I show how effect size varied based on the percentage of the sample that is male.
The estimate for “male” is positive (est.: 0.705, p=0.282), which indicates that having a higher percentage of male participants was associated with larger positive effects on QoL, on average. However, the p-value is greater than 0.05, so this finding is not statistically significant.
# meta-regression with male
metaregression <- robu(formula = es ~ male, # es ~ male shows how the average
# effect size depends on percentage
# of participants that are male
data = data,
studynum = author,
var.eff.size = var,
rho = 0.8, small = TRUE)
print(metaregression)## RVE: Correlated Effects Model with Small-Sample Corrections
##
## Model: es ~ male
##
## Number of studies = 21
## Number of outcomes = 59 (min = 1 , mean = 2.81 , median = 3 , max = 5 )
## Rho = 0.8
## I.sq = 95.9774
## Tau.sq = 0.1716802
##
## Estimate StdErr t-value dfs P(|t|>) 95% CI.L 95% CI.U Sig
## 1 X.Intercept. 0.197 0.461 0.428 8.95 0.679 -0.846 1.24
## 2 male 0.705 0.618 1.140 9.50 0.282 -0.682 2.09
## ---
## Signif. codes: < .01 *** < .05 ** < .10 *
## ---
## Note: If df < 4, do not trust the results
Most notably, it’s valuable for creating forest plots to visualize meta-analysis data. I will create a forest plot for the 7 studies that reported environmental QoL as an outcome (because the forest plot would be overwhelming with 60 effect sizes). The resulting forest plot shows the effect size for each individual study (black square), as well as the average effect across the studies (the white diamond at the bottom).
# First, I subsetted the data to only include studies that reported an
# environmental QoL outcome measure.
environmentaldata <- data %>%
filter(domaincategory == "Environmental")
# Next, I created a regression model to calculate the average effect of
# buprenorphine on environmental QoL.
forest <- robu(formula = es ~ 1,
data = environmentaldata,
studynum = author,
var.eff.size = var,
rho = 0.8, small = TRUE)
# Lastly, I created the forest plot.
forest.robu(forest, es.lab = "es", study.lab = "author")
Learn more about robumeta and meta-analysis with the following:
Resource I Cran: “robumeta: Robust Variance Meta-Regression”
Resource II Rdocumentation.org: “robumeta (version 2.0)”
Resource III Elizabeth Tipton: “Meta-Analysis”
This code through references and cites the following sources:
Golan, O.K. et al. (2022). Source I. Systematic review and meta-analysis of changes in quality of life following initiation of buprenorphine for opioid use disorder
Fisher, Z. & Tipton, E. (2015). Source II. Robumeta: An R-package for robust variance estimation in meta-analysis
Tipton, E. (2017). Source III. forest. robu 3