This is the software workflow for our Introduction to Meta-Analysis workshop. The slides for the presentation accompanying this workflow can be found here.
In this workflow we focus on practical data analysis by presenting the steps to run a meta-analysis using the metafor package in R. The R code used to create output is also included.
If you use this workflow, please don’t forget to acknowledge us. You can use the following sentence
The authors acknowledge the Statistical Workshops and Workflows provided by the Sydney Informatics Hub, a Core Research Facility of the University of Sydney
This data is taken from the book “Methods of Meta-analysis in Medical Research” by A. J. Sutton (originally extracted from the article: Smith, George Davey, Fujian Song, and Trevor A Sheldon. “Cholesterol Lowering and Mortality: The Importance of Considering Initial Level of Risk.” British Medical Journal 306.6889 (1993): 1367–. Print.)
The data comprises 34 studies assessing the effectiveness of cholesterol treatments on mortality in randomised controlled trials.
Before we delve into the content, it is necessary to load certain libraries. If you want to use them, please install them beforehand with install.packages("library_name").
# Run once
# install.packages("metafor")
library(metafor)Loading required package: MatrixLoading required package: metadatLoading required package: numDeriv
Loading the 'metafor' package (version 4.8-0). For an
introduction to the package please type: help(metafor)Load the data. Note: check file is in your working directory.
metadata <- read.csv("Meta_Sutton_Smith.csv", header = T)
head(metadata) Sutton.et.al author year g1.r1 g1.r0 g1.n g0.r1 g0.r0 g0.n
1 1 author 1 1992 28 176 204 51 151 202
2 2 author 2 1962 70 215 285 38 109 147
3 3 author 3 1963 37 119 156 40 79 119
4 4 author 4 1981 2 86 88 3 27 30
5 5 author 5 1969 0 30 30 3 30 33
6 6 author 6 1988 61 218 279 82 194 276This function creates variables for yi which is the effect size, and vi which is the variance of yi. These are added to the metadata object dataframe.
Tell the function what effect size to calculate using measure=“OR”,“RR”,“RD”, etc.
We can only choose a measure that matches our data - i.e. it has a binary outcome.
meta.es <-
escalc(measure = "OR", g1.r1, g1.r0, g0.r1, g0.r0, data = metadata)
head(meta.es)
Sutton.et.al author year g1.r1 g1.r0 g1.n g0.r1 g0.r0 g0.n yi vi
1 1 author 1 1992 28 176 204 51 151 202 -0.7528 0.0676
2 2 author 2 1962 70 215 285 38 109 147 -0.0684 0.0544
3 3 author 3 1963 37 119 156 40 79 119 -0.4876 0.0731
4 4 author 4 1981 2 86 88 3 27 30 -1.5640 0.8820
5 5 author 5 1969 0 30 30 3 30 33 -1.9459 2.3513
6 6 author 6 1988 61 218 279 82 194 276 -0.4125 0.0383 We use the escalc function from metafor to create yi which is the log OR effect size and vi which is the variance of yi.
When we view the model object meta.es, the effect size and variance sit alongside our inputs.
meta.re <- rma(yi, vi, data = meta.es) #REML estimation by default
# Use the method argument to change estimation to Der Simonian and Laird, etc.
# for example: meta.reDL <-rma(yi, vi, data = meta.es, method = "DL")
# Now convert summary estimate log(OR) to OR
meta.re.or <- exp(c(meta.re$b, meta.re$ci.lb, meta.re$ci.ub))
summary(meta.re)
Random-Effects Model (k = 34; tau^2 estimator: REML)
logLik deviance AIC BIC AICc
-22.9313 45.8625 49.8625 52.8556 50.2625
tau^2 (estimated amount of total heterogeneity): 0.0478 (SE = 0.0269)
tau (square root of estimated tau^2 value): 0.2186
I^2 (total heterogeneity / total variability): 53.87%
H^2 (total variability / sampling variability): 2.17
Test for Heterogeneity:
Q(df = 33) = 89.1065, p-val < .0001
Model Results:
estimate se zval pval ci.lb ci.ub
-0.1178 0.0628 -1.8750 0.0608 -0.2410 0.0053 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Print OR results
cat("Odds Ratio (OR):", round(meta.re.or[1], 3), "\n","95% CI: [", round(meta.re.or[2], 3), ",", round(meta.re.or[3], 3), "]\n")Odds Ratio (OR): 0.889
95% CI: [ 0.786 , 1.005 ]# For a fixed effects model, use method="FE"
meta.fe <- rma(yi, vi, data = meta.es, method = "FE")
summary(meta.fe)
Fixed-Effects Model (k = 34)
logLik deviance AIC BIC AICc
-41.3091 89.1065 84.6181 86.1445 84.7431
I^2 (total heterogeneity / total variability): 62.97%
H^2 (total variability / sampling variability): 2.70
Test for Heterogeneity:
Q(df = 33) = 89.1065, p-val < .0001
Model Results:
estimate se zval pval ci.lb ci.ub
-0.1691 0.0324 -5.2258 <.0001 -0.2325 -0.1057 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Now convert summary estimate log(OR) to OR
meta.fe.or <- exp(c(meta.fe$b, meta.fe$ci.lb, meta.fe$ci.ub))
cat("Odds Ratio (OR):", round(meta.fe.or[1], 3), "95% CI: [", round(meta.fe.or[2], 3), ",", round(meta.fe.or[3], 3), "]\n")Odds Ratio (OR): 0.844 95% CI: [ 0.793 , 0.9 ]Random-Effects model has
# Create Forest Plots
# First plot is for Random-Effects model
forest(
meta.re,
header = TRUE,
atransf = exp, #display OR on x-axis log scale
at = log(c(0.01, 0.1, 1, 10,50)),
xlim = c(-8, 8),
xlab = "Random Effects Model",
mlab = "Random Effects OR summary estimate",
showweights = TRUE,
pch = 15,
border = "red", #colour of summary diamond border
col = "red", #colour of summary diamond fill
main = "The effect of cholesterol lowering interventions \n on mortality \n Sutton-Smith data"
)
op <- par(cex = 0.75, font = 2)
text(c(-2, 2), 36, c("Favours Treatment", "Favours Control"))
par(op)
# Second plot for Fixed-Effect model
forest(
meta.fe,
header = TRUE,
atransf = exp, #display OR on x-axis log scale
at = log(c(0.01, 0.1, 1, 10, 50)),
xlim = c(-8, 8),
xlab = "Fixed-Effect Model",
mlab = "Fixed-Effect OR summary estimate",
showweights = TRUE,
pch = 15,
border = "red", #colour of summary diamond border
col = "red", #colour of summary diamond fill
main = "The effect of cholesterol lowering interventions \n on mortality \n Sutton-Smith data"
)
op <- par(cex = 0.75, font = 2)
text(c(-2, 2), 36, c("Favours Treatment", "Favours Control"))
par(op)Refer to the output from the rma() function for both models.
No additional code required.
# Basic Funnel Plot
funnel(meta.re, main = "Random Effects Model")
# We can use trimfill(), ranktest(), regtest(), hc() as tests for bias
# Use the trimfill() function to see where "missing" studies might beCheck for asymmetry * No of studies > 10 * Studies are of various sizes and precision * Plot appears to be symmetrical (can use random effects model without problem) * One more check
meta.tf <- trimfill(meta.re)
funnel(meta.tf, main = "Trim and Fill funnel plot \n for RE model")
Note: 3 imputed studies (open circles)
# Other diagnostic plots like Radial (Galbraith) plot are available
# aslo labbe(), qqnorm(), baujat(), gosh()
# lets look at the Radial plots
radial(meta.re, main = "Random Effects model")
radial(meta.fe, main = "Fixed Effects model")
# lets look at some Q-Q plots
qqnorm(meta.re, main = "Random Effects Model")
qqnorm(meta.fe, main = "Fixed Effects Model")
# Can also produce influence statistics
infl <- influence(meta.re)
plot(infl, plotdfb = TRUE)
The original study looked at a sub-group analysis using the factor: Coronary heart disease risk group – high, medium, low
Sub-group analysis found significant differences between these groups, with only the high risk CHD (treatment) group showing a benefit in terms of total mortality.
A meta-regression was also carried out showing a similar trend.