Overview

  • describe and illustrate the ‘standard workflow’ of a meta-analysis conducted with the metafor package (Viechtbauer, 2010)
  • discuss some alternative modeling approaches
  • illustrate a multivariate meta-analysis

 

Setup

# install the metafor package
install.packages("metafor")

# load the metafor package
library(metafor)

# for the example data, we also need the bayesmeta package
install.packages("bayesmeta")

 

Standard Workflow

  • define the goal(s) of the meta-analysis and inclusion/exclusion criteria
  • find relevant studies that have examined the phenomenon of interest
  • quantify the results in terms of an effect size measure
  • quantify the precision of the estimates in terms of their variances
  • meta-analyze the estimates (equal/fixed/random-effects model)

 

Example

  • meta-analysis on the effectiveness of tiotropium for the management of COPD (Karner et al., 2014)
  • study participants were randomly assigned to either a treatment (tiotropium) or a control (placebo) group
  • collected information on the number of patients in each group with exacerbations and the number of deaths during follow-up

 

Effect Size Calculation

  • escalc() allows for the calculation of 70+ different types of effect size measures / outcome measures
  • example: log odds ratio for ≥1 exacerbations in the treatment group compared to the control group
# copy the data to 'dat'
dat <- get(data(KarnerEtAl2014, package="bayesmeta"))

# examine the dataset (for a selection of variables)
dat[c(1,10:11,17:18)]
#             study tiotropium.total tiotropium.exa placebo.total placebo.exa
# 1    Bateman2010a             1989            685          2002         842
# 2    Bateman2010b             1337            495           653         288
# 3        Beeh2006             1236            180           403          80
# 4    Brusasco2003              402            129           400         156
# 5    Casaburi2002              550            198           371         156
# 6        Chan2007              608            268           305         125
# 7      Cooper2010              260            112           259         102
# 8     Covelli2005              100              9            96          12
# 9      Dusser2006              500            213           510         272
# 10    Freeman2007              200             19           195          35
# 11  Johansson2008              107              2           117           4
# 12  Magnussen2008              228             13           244          26
# 13      Moita2008              147              6           164           6
# 14    NCT00144326              123             11           127          12
# 15 Niewoehner2005              914            255           915         296
# 16     Powrie2007               69             30            73          47
# 17        Sun2007               30              0            30           2
# 18    Tashkin2008             2987           2001          3006        2049
# 19     Tonnel2008              266            101           288         130
# 20   Trooster2011              238             11           219          24
# 21  Verkindre2006               46             10            54           8
# 22    Voshaar2008              360             43           181          21

# tiotropium.total - total number of patients in the treatment group
# tiotropium.exa   - patients with ≥1 exacerbations in the treatment group
# placebo.total    - total number of patients in the placebo group
# placebo.exa      - patients with ≥1 exacerbations in the placebo group

# calculate log odds ratios and corresponding sampling variances
dat <- escalc(measure="OR", ai=tiotropium.exa, n1i=tiotropium.total,
                            ci=placebo.exa,    n2i=placebo.total,
                            slab=study, data=dat)

# examine the dataset
cbind(dat[1], "..."="...", dat[28:29])
#             study ...      yi     vi 
# 1    Bateman2010a ... -0.3234 0.0043 
# 2    Bateman2010b ... -0.2943 0.0094 
# 3        Beeh2006 ... -0.3737 0.0221 
# 4    Brusasco2003 ... -0.3023 0.0219 
# 5    Casaburi2002 ... -0.2546 0.0190 
# 6        Chan2007 ...  0.1267 0.0202 
# 7      Cooper2010 ...  0.1526 0.0319 
# 8     Covelli2005 ... -0.3677 0.2173 
# 9      Dusser2006 ... -0.4317 0.0161 
# 10    Freeman2007 ... -0.7342 0.0930 
# 11  Johansson2008 ... -0.6197 0.7684 
# 12  Magnussen2008 ... -0.6793 0.1246 
# 13      Moita2008 ...  0.1138 0.3468 
# 14    NCT00144326 ... -0.0606 0.1919 
# 15 Niewoehner2005 ... -0.2117 0.0104 
# 16     Powrie2007 ... -0.8544 0.1187 
# 17        Sun2007 ... -1.6773 2.4679 
# 18    Tashkin2008 ... -0.0536 0.0030 
# 19     Tonnel2008 ... -0.2958 0.0300 
# 20   Trooster2011 ... -0.9321 0.1421 
# 21  Verkindre2006 ...  0.4683 0.2745 
# 22    Voshaar2008 ...  0.0329 0.0803

 

  • yi – log odds ratios (negative values indicate lower odds of exacerbations in the treatment group)
  • vi – corresponding sampling variances

 

Random-Effects Model

  • rma() is one of the main modeling functions for fitting equal-, fixed-, random-, and mixed-effects meta-regression models
  • for random/mixed-effects models, fits ‘normal-normal models’
  • uses restricted maximum likelihood estimation (REML) by default
  • lots of alternative estimators for \(\tau^2\)
# fit a random-effects model (with K&H method)
res <- rma(yi, vi, data=dat, test="knha")
res
# Random-Effects Model (k = 22; tau^2 estimator: REML)
# 
# tau^2 (estimated amount of total heterogeneity): 0.0204 (SE = 0.0150)
# tau (square root of estimated tau^2 value):      0.1427
# I^2 (total heterogeneity / total variability):   48.61%
# H^2 (total variability / sampling variability):  1.95
# 
# Test for Heterogeneity:
# Q(df = 21) = 42.4592, p-val = 0.0037
# 
# Model Results:
# 
# estimate      se     tval  df    pval    ci.lb    ci.ub      
#  -0.2451  0.0544  -4.5076  21  0.0002  -0.3582  -0.1320  *** 
# 
# ---
# Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# back-transform results to the odds ratio scale
predict(res, transf=exp, digits=2)
#  pred ci.lb ci.ub pi.lb pi.ub 
#  0.78  0.70  0.88  0.57  1.08

# forest plot
par(mar=c(5,4,2,2))
forest(res, xlim=c(-5,3.8), atransf=exp, at=log(2^(-4:3)),
       addpred=TRUE, digits=c(2L,4L), header=TRUE, top=2, shade=TRUE)

  • the look of forest plots is very customizable, but it can take some work to replicate the look of other software (see examples here)
# risk-of-bias ratings extracted from Figure 2 from Karner et al (2014)
rob <- read.table(header=TRUE, text = "
rb.a rb.b rb.c rb.d rb.e rb.f
'+' '+' '+' '+' '+' '+'
'+' '+' '+' '+' '-' '+'
'+' '+' '+' '+' '?' '+'
'+' '+' '+' '+' '?' '+'
'+' '+' '+' '+' '?' '+'
'+' '+' '+' '+' '?' '+'
'+' '+' '+' '+' '-' '+'
'+' '+' '+' '+' '?' '?'
'+' '+' '+' '+' '?' '+'
'+' '+' '+' '+' '?' '+'
'+' '+' '+' '+' '+' '+'
'+' '+' '+' '+' '+' '+'
'+' '+' '+' '+' '+' '+'
'+' '+' '+' '+' '+' '+'
'+' '+' '+' '+' '?' '+'
'+' '+' '+' '+' '?' '+'
'+' '?' '+' '?' '+' '+'
'+' '+' '+' '+' '-' '+'
'+' '+' '+' '+' '?' '+'
'+' '+' '+' '+' '+' '+'
'+' '+' '+' '+' '?' '+'
'+' '+' '+' '+' '+' '+'")

# turn the risk of bias items into factors with levels +, -, and ?
rob <- lapply(rob, factor, levels=c("+", "-", "?"))

# estimated average odds ratio (and 95% CI/PI)
pred <- predict(res, transf=exp, digits=2)

# need the rounded estimate and CI bounds further below
pred <- fmtx(c(pred$pred, pred$ci.lb, pred$ci.ub), 2)

# total number of studies
k <- nrow(dat)

# get the weights and format them as will be used in the forest plot
weights <- paste0(fmtx(weights(res), 1), "%")

# adjust the margins
par(mar=c(10.6,0,2.3,1.3), mgp=c(3,0.4,0), tcl=-0.3)

# forest plot with extra annotations
sav <- forest(res, atransf=exp, at=log(c(0.1, 0.2, 0.5, 1, 2, 5, 10)), xlim=c(-15,5),
       xlab="", efac=c(0,1,1), textpos=c(-15,-2.5), lty=c(1,1,0), refline=NA,
       ilab=cbind(tiotropium.exa, tiotropium.total, placebo.exa, placebo.total, weights),
       ilab.xpos=c(-10.8,-9.6,-8.1,-6.9,-5.2), ilab.pos=2,
       cex=0.79, header=c("Study","95% CI"), mlab="")

# add the horizontal line at the top
segments(sav$xlim[1], k+1, sav$xlim[2], k+1)

# add the vertical reference line at 0
segments(0, -3, 0, k+1)

# now we add a bunch of text; since some of the text falls outside of the
# plot region, we set xpd=NA so nothing gets clipped
par(xpd=NA)

# adjust cex as used in the forest plot and use a bold font
par(cex=sav$cex, font=2)

text(sav$ilab.xpos, k+2, pos=2, c("Events","Total","Events","Total","Weight"))
text(c(mean(sav$ilab.xpos[1:2])+0.2,mean(sav$ilab.xpos[3:4])), k+3, pos=2, c("tiotropium","placebo"))
text(sav$textpos[2], k+3, "Odds ratio", pos=2)
text(0, k+3, "Odds ratio")
text(sav$xlim[2]-0.35, k+3, "Risk of Bias", pos=2)
text(0, k+2, "IV, Random, 95% CI")
text(c(sav$xlim[1],sav$ilab.xpos[c(2,4,5)]), -1, pos=c(4,2,2,2,2),
     c("Total (95% CI)", sum(dat$tiotropium.total), sum(dat$placebo.total), "100.0%"))
text(sav$xlim[1], -7, pos=4, "Risk of bias legend")

# first hide the non-bold summary estimate text and then add it back in bold font
rect(sav$textpos[2], -1.5, sav$ilab.xpos[5], -0.5, col="white", border=NA)
text(sav$textpos[2], -1, paste0(pred[1], " [", pred[2], ",  ", pred[3], "]"), pos=2)

# use a non-bold font for the rest of the text
par(cex=sav$cex, font=1)

# add favours text below the x-axis
text(log(c(1,1)), -4.5, c("Favours tiotropium","Favours placebo"), pos=c(2,4), offset=0.5)

# add text for total events
text(sav$xlim[1], -2, pos=4, "Total events:")
text(sav$ilab.xpos[c(1,3)], -2, c(sum(dat$tiotropium.exa),sum(dat$placebo.exa)), pos=2)

# add text with heterogeneity statistics
text(sav$xlim[1], -3, pos=4, bquote(paste("Heterogeneity: ",
   "Tau"^2, " = ", .(fmtx(res$tau2, 2)), "; ",
   "Chi"^2, " = ", .(fmtx(res$QE, 2)),
   ", df = ", .(res$k - res$p),
   " (", .(fmtp(res$QEp, digits=3, pname="P", add0=TRUE, sep=TRUE, equal=TRUE)), "); ",
   I^2, " = ", .(round(res$I2)), "%")))

# add text for test of overall effect
text(sav$xlim[1], -4, pos=4, bquote(paste("Test for overall effect: Z = ",
     .(fmtx(res$zval, 2)), " (", .(fmtp(res$pval, digits=5, pname="P", add0=TRUE, sep=TRUE, equal=TRUE)), ")")))

# add risk of bias points and symbols
cols <- c("#00cc00", "#cc0000", "#eeee00")
syms <- levels(rob$rb.a)
pos  <- seq(sav$xlim[2]-2.4,sav$xlim[2]-0.5,length=6)
for (i in 1:6) {
   points(rep(pos[i],k), k:1, pch=19, col=cols[rob[[i]]], cex=2.2)
   text(pos[i], k:1, syms[rob[[i]]], font=2)
}
text(pos, k+2, c("A","B","C","D","E","F"), font=2)

# add risk of bias legend
text(sav$xlim[1], -8:-13, pos=4, c(
   "(A) Random sequence generation (selection bias)",
   "(B) Allocation concealment (selection bias)",
   "(C) Blinding of participants and personnel (performance bias)",
   "(D) Blinding of outcome assessment (detection bias)",
   "(E) Incomplete outcome data (attrition bias)",
   "(F) Selective reporting (reporting bias)"))