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)"))

 

Publication Bias

  • the presence of publication bias (or more accurately, funnel plot asymmetry or “small-study effects”) and its potential impact on the results can be examined via a variety of methods
# funnel plot
par(mar=c(5,4,2,2))
funnel(res, atransf=exp, at=log(2^(-5:5)), ylim=c(0,2),
       las=1, digits=list(5L,1))

# regression test for funnel plot asymmetry
regtest(res)
# Regression Test for Funnel Plot Asymmetry
# 
# Model:     mixed-effects meta-regression model
# Predictor: standard error
# 
# Test for Funnel Plot Asymmetry: t = -1.1076, df = 20, p = 0.2812
# Limit Estimate (as sei -> 0):   b = -0.1675 (CI: -0.3509, 0.0159)

# rank correlation test for funnel plot asymmetry
ranktest(res)
# Rank Correlation Test for Funnel Plot Asymmetry
# 
# Kendall's tau = -0.0130, p = 0.9556

# trim-and-fill method
taf <- trimfill(res)
taf
# Estimated number of missing studies on the right side: 3 (SE = 3.1378)
# 
# Random-Effects Model (k = 25; tau^2 estimator: REML)
# 
# tau^2 (estimated amount of total heterogeneity): 0.0239 (SE = 0.0165)
# tau (square root of estimated tau^2 value):      0.1546
# I^2 (total heterogeneity / total variability):   49.75%
# H^2 (total variability / sampling variability):  1.99
# 
# Test for Heterogeneity:
# Q(df = 24) = 50.0371, p-val = 0.0014
# 
# Model Results:
# 
# estimate      se     zval    pval    ci.lb    ci.ub      
#  -0.2194  0.0534  -4.1089  <.0001  -0.3241  -0.1148  *** 
# 
# ---
# Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# funnel plot
par(mar=c(5,4,2,2))
funnel(taf, atransf=exp, at=log(2^(-5:5)), ylim=c(0,2),
       las=1, digits=list(5L,1))
legend("topright", inset=.01, pch=c(19,21), pt.bg="gray90",
       legend=c("Observed Studies", "Filled Studies"))

# contour-enhanced funnel plot (zoomed in at the top)
par(mar=c(5,4,2,2))
funnel(res, atransf=exp, at=log(2^(-2:2)), xlim=c(-1,1), ylim=c(0,0.6),
       las=1, digits=list(2L,1), refline=0,
       level=c(90,95), shade=c("white","gray75"), legend=TRUE)

# test of excess significance
tes(res)
# Test of Excess Significance
# 
# Observed Number of Significant Findings: 9 (out of 22)
# Expected Number of Significant Findings: 6.9080
# Observed Number / Expected Number:       1.3028
# 
# Estimated Power of Tests (based on theta = -0.2451)
# 
#    min      q1  median      q3     max 
# 0.0537  0.1042  0.2489  0.4509  0.8251 
# 
# Test of Excess Significance: p = 0.1683 (X^2 = 0.9235, df = 1)
# Limit Estimate (theta_lim):  -0.2187 (where p = 0.1)

# fit step function selection model
sel <- selmodel(res, type="stepfun", alternative="less", steps=c(.025,.05,.5,1))
sel
# Random-Effects Model (k = 22; tau^2 estimator: ML)
# 
# tau^2 (estimated amount of total heterogeneity): 0.0144 (SE = 0.0124)
# tau (square root of estimated tau^2 value):      0.1200
# 
# Test for Heterogeneity:
# LRT(df = 1) = 6.3227, p-val = 0.0119
# 
# Model Results:
# 
# estimate      se     zval    pval    ci.lb   ci.ub    
#  -0.1264  0.0887  -1.4260  0.1539  -0.3002  0.0473    
# 
# Test for Selection Model Parameters:
# LRT(df = 3) = 7.2109, p-val = 0.0655
# 
# Selection Model Results:
# 
#                     k  estimate      se     zval    pval   ci.lb   ci.ub      
# 0     < p <= 0.025  9    1.0000     ---      ---     ---     ---     ---      
# 0.025 < p <= 0.05   3    0.8065  0.5991  -0.3230  0.7467  0.0000  1.9807      
# 0.05  < p <= 0.5    5    0.1468  0.1161  -7.3478  <.0001  0.0000  0.3744  *** 
# 0.5   < p <= 1      5    0.2096  0.2216  -3.5675  0.0004  0.0000  0.6438  *** 
# 
# ---
# Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# plot selection function
plot(sel, ylim=c(0,1), las=1, bty="l")

 

Meta-Regression

  • can include one or multiple quantitative and/or categorical predictors (‘moderators’) in the model
# meta-regression model using the mean FEV1 at baseline as predictor
res <- rma(yi, vi, mods = ~ baseline.fev1, data=dat, test="knha")
res
# Mixed-Effects Model (k = 22; tau^2 estimator: REML)
# 
# tau^2 (estimated amount of residual heterogeneity):     0.0109 (SE = 0.0107)
# tau (square root of estimated tau^2 value):             0.1044
# I^2 (residual heterogeneity / unaccounted variability): 33.77%
# H^2 (unaccounted variability / sampling variability):   1.51
# R^2 (amount of heterogeneity accounted for):            46.40%
# 
# Test for Residual Heterogeneity:
# QE(df = 20) = 29.6562, p-val = 0.0756
# 
# Test of Moderators (coefficient 2):
# F(df1 = 1, df2 = 20) = 11.4396, p-val = 0.0030
# 
# Model Results:
# 
#                estimate      se     tval  df    pval    ci.lb    ci.ub     
# intrcpt          0.6740  0.2728   2.4708  20  0.0226   0.1050   1.2430   * 
# baseline.fev1   -0.7935  0.2346  -3.3822  20  0.0030  -1.2829  -0.3041  ** 
# 
# ---
# Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# scatter/bubbeplot of the results
par(mar=c(5,4,2,2))
tmp <- regplot(res, xlab="Baseline Mean FEV1", bty="l", las=1, digits=1, legend=TRUE)
# label some points
lab <- c(11,17,18,20)
text(tmp$xi[lab], tmp$yi[lab], tmp$ids[lab], pos=c(1,3,4,1), offset=0.6)

# predicted pooled odds ratio at FEV1 equal to 1, 1.5, and 2
predict(res, newmods=c(1, 1.5, 2), transf=exp, digits=2)
#   pred ci.lb ci.ub pi.lb pi.ub 
# 1 0.89  0.79  0.99  0.69  1.13 
# 2 0.60  0.49  0.72  0.45  0.80 
# 3 0.40  0.26  0.61  0.25  0.65

# plot of various influence / outlier diagnostics
plot(influence(res), las=2)

 

Other Methods/Models

  • rma.mh() - meta-analysis via the Mantel-Haenszel method
  • rma.peto() - meta-analysis via Peto’s one-step method
  • rma.glmm() - meta-analysis via generalized linear models
# examine the dataset (outcome: mortality)
dat[c(10,14,17,21)]
#    tiotropium.total tiotropium.deaths placebo.total placebo.deaths 
# 1              1989                52          2002             38 
# 2              1337                34           653             10 
# 3              1236                 2           403              2 
# 4               402                 1           400              5 
# 5               550                 7           371              7 
# 6               608                13           305              2 
# 7               260                 6           259              6 
# 8               100                 0            96              0 
# 9               500                 7           510              8 
# 10              200                 1           195              4 
# 11              107                 0           117              1 
# 12              228                 0           244              2 
# 13              147                 2           164              0 
# 14              123                 0           127              0 
# 15              914                22           915             19 
# 16               69                 1            73              2 
# 17               30                 0            30              0 
# 18             2987               381          3006            411 
# 19              266                 3           288              6 
# 20              238                 0           219              0 
# 21               46                 0            54              0 
# 22              360                 2           181              0

# fit a random-effects model using a 'binomial-normal model'
res <- rma.glmm(measure="OR", data=dat,
                ai=tiotropium.deaths, n1i=tiotropium.total,
                ci=placebo.deaths,    n2i=placebo.total)
# Warning: 5 studies with NAs omitted from model fitting.
# Warning: Some yi/vi values are NA.
res
# Random-Effects Model (k = 17; tau^2 estimator: ML)
# Model Type: Unconditional Model with Fixed Study Effects
# 
# tau^2 (estimated amount of total heterogeneity): 0
# tau (square root of estimated tau^2 value):      0
# I^2 (total heterogeneity / total variability):   0.00%
# H^2 (total variability / sampling variability):  1.00
# 
# Tests for Heterogeneity:
# Wld(df = 16) = 14.6596, p-val = 0.5497
# LRT(df = 16) = 25.4335, p-val = 0.0625
# 
# Model Results:
# 
# estimate      se     zval    pval    ci.lb   ci.ub    
#  -0.0234  0.0653  -0.3584  0.7200  -0.1513  0.1045    
# 
# ---
# Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

predict(res, transf=exp, digits=2)
#  pred ci.lb ci.ub pi.lb pi.ub 
#  0.98  0.86  1.11  0.86  1.11

 

Multivariate Models

  • rma.mv() can be used to fit multilevel and multivariate models
  • (also covers network meta-analyses, phylogenetic meta-analyses, diagnostic test meta-analyses, spatio-temporal models, etc.)
  • let \(y_{ij}\) denote the \(j\)th observed effect size in the \(i\)th study
  • assume \(j\) denotes different types of constructs / outcomes
  • the multivariate random-effects model (e.g., Berkey et al., 1998; Kalaian & Raudenbush, 1996) is given by \[y_{ij} = \mu_j + u_{ij} + \varepsilon_{ij}\] where \(\left[ \begin{array}{c} u_{i1} \\ u_{i2} \\ \vdots \end{array} \right] \sim N\left( \left[ \begin{array}{c} 0 \\ 0 \\ \vdots \end{array} \right], \left[ \begin{array}{ccc} \tau_1^2 & \rho_{12}\tau_1\tau_2 & \ldots \\ \rho_{12}\tau_1\tau_2 & \tau_2^2 & \ldots \\ \vdots & \vdots & \ddots \end{array} \right] \right)\)

     

    and \(\left[ \begin{array}{c} \varepsilon_{i1} \\ \varepsilon_{i2} \\ \vdots \end{array} \right] \sim N\left( \left[ \begin{array}{c} 0 \\ 0 \\ \vdots \end{array} \right], \left[ \begin{array}{ccc} v_{i1} & cov_{i12} & \ldots \\ cov_{i12} & v_{i2} & \ldots \\ \vdots & \vdots & \ddots \end{array} \right] \right)\)

     

  • try this out for outcomes ‘exacerbations’ and ‘severe exacerbations’
# copy the data to 'dat' again
dat <- get(data(KarnerEtAl2014, package="bayesmeta"))

# select only relevant variables (also remove two studies)
dat <- dat[-c(17,20),c(1,10:12,17:19)]

# abbreviate the variable names
names(dat) <- abbreviate(names(dat))

# examine the dataset
dat
#              stdy ttrpm.t ttrpm.x ttrpm.s plcb.t plcb.x plcb.s
# 1    Bateman2010a    1989     685     161   2002    842    198
# 2    Bateman2010b    1337     495      78    653    288     44
# 3        Beeh2006    1236     180      29    403     80      7
# 4    Brusasco2003     402     129      48    400    156     90
# 5    Casaburi2002     550     198      30    371    156     35
# 6        Chan2007     608     268      51    305    125     19
# 7      Cooper2010     260     112      21    259    102     16
# 8     Covelli2005     100       9       2     96     12      1
# 9      Dusser2006     500     213      28    510    272     33
# 10    Freeman2007     200      19       2    195     35      1
# 11  Johansson2008     107       2       0    117      4      0
# 12  Magnussen2008     228      13       4    244     26      0
# 13      Moita2008     147       6       1    164      6      1
# 14    NCT00144326     123      11       2    127     12      0
# 15 Niewoehner2005     914     255      64    915    296     87
# 16     Powrie2007      69      30       2     73     47      3
# 18    Tashkin2008    2987    2001     759   3006   2049    811
# 19     Tonnel2008     266     101      11    288    130      8
# 21  Verkindre2006      46      10       0     54      8      3
# 22    Voshaar2008     360      43       2    181     21      0

# ttrpm.t - total number of patients in the treatment group
# ttrpm.x - patients with ≥1 exacerbations in the treatment group
# ttrpm.s - patients with ≥1 severe exacerbations in the treatment group
# plcb.t  - total number of patients in the placebo group
# plcb.x  - patients with ≥1 exacerbations in the placebo group
# plcb.s  - patients with ≥1 severe exacerbations in the placebo group

# calculate log odds ratios for exacerbations
dat <- escalc(measure="OR", ai=ttrpm.x, n1i=ttrpm.t,
                            ci=plcb.x,  n2i=plcb.t,
                            data=dat, var.names=c("y1i","v1i"))

# calculate log odds ratios for severe exacerbations
dat <- escalc(measure="OR", ai=ttrpm.s, n1i=ttrpm.t,
                            ci=plcb.s,  n2i=plcb.t,
                            data=dat, var.names=c("y2i","v2i"))
  • because y1i and y2i are calculate based on the same sample of subjects, they are not independent
  • because ttrpm.s is a subset of ttrpm.x and plcb.s is a subset of plcb.x, we can calculate their covariances (cov12i) without further information (Trikalinos & Olkin, 2012, see Appendix)
# calculate the covariances
dat$cov12i <- with(dat, ttrpm.t / (ttrpm.x * (ttrpm.t-ttrpm.s)) +
                        plcb.t  / (plcb.x * (plcb.t-plcb.s)))

# plot of the log odds ratios against each other (with 50% confidence ellipses)
# set up the plot
par(mar=c(5,4,2,2))
plot(NA, xlim=c(-1.7,1), ylim=c(-3.6,4), bty="l", xlab="Log(OR) for Exacerbations",
     ylab="Log(OR) for Severe Exacerbations")
# add the confidence ellipses (requires that the 'ellipse' package is loaded)
invisible(sapply(1:nrow(dat), function(i) {
   x <- dat[i,]
   xy <- ellipse(matrix(c(x$v1i,x$cov12i,x$cov12i,x$v2i), nrow=2),
                 centre=c(x$y1i,x$y2i), level=0.5)
   lines(xy[,1],xy[,2], col="gray80")
}))
# add the points
points(dat$y1i, dat$y2i, pch=21, bg="gray60", cex=1.5)

# reshape dataset into long format
dat.long <- data.frame(study   = rep(dat$stdy, each=2),
                       outcome = c("exa","sev"),
                       yi      = c(t(dat[c("y1i","y2i")])),
                       vi      = c(t(dat[c("v1i","v2i")])))

# examine the first six rows of the reshaped dataset
head(dat.long)
#          study outcome         yi          vi
# 1 Bateman2010a     exa -0.3233776 0.004276443
# 2 Bateman2010a     sev -0.2200787 0.012363055
# 3 Bateman2010b     exa -0.2942854 0.009419799
# 4 Bateman2010b     sev -0.1537356 0.037984103
# 5     Beeh2006     exa -0.3736609 0.022098500
# 6     Beeh2006     sev  0.3069067 0.180693654

# construct the 2x2 var-cov matrices of the studies
V <- tapply(dat, dat$stdy, function(x) matrix(c(x$v1i, x$cov12i,
                                                x$cov12i, x$v2i), nrow=2))

# examine the first three elements
V[1:3]
# $Bateman2010a
#             [,1]       [,2]
# [1,] 0.004276443 0.00290643
# [2,] 0.002906430 0.01236305
# 
# $Bateman2010b
#             [,1]       [,2]
# [1,] 0.009419799 0.00586845
# [2,] 0.005868450 0.03798410
# 
# $Beeh2006
#           [,1]      [,2]
# [1,] 0.0220985 0.0184100
# [2,] 0.0184100 0.1806937

# fit a multivariate random-effects model
res <- rma.mv(yi, V, mods = ~ factor(outcome) - 1,
              random = ~ outcome | study, struct="UN",
              data=dat.long)
res
# Multivariate Meta-Analysis Model (k = 40; method: REML)
# 
# Variance Components:
# 
# outer factor: study   (nlvls = 20)
# inner factor: outcome (nlvls = 2)
# 
#             estim    sqrt  k.lvl  fixed  level 
# tau^2.1    0.0187  0.1369     20     no    exa 
# tau^2.2    0.0494  0.2222     20     no    sev 
# 
#      rho.exa  rho.sev    exa  sev 
# exa        1               -   20 
# sev   0.4931        1     no    - 
# 
# Test for Residual Heterogeneity:
# QE(df = 38) = 68.8813, p-val = 0.0016
# 
# Test of Moderators (coefficients 1:2):
# QM(df = 2) = 20.2765, p-val < .0001
# 
# Model Results:
# 
#                     estimate      se     zval    pval    ci.lb    ci.ub      
# factor(outcome)exa   -0.2311  0.0513  -4.5018  <.0001  -0.3317  -0.1305  *** 
# factor(outcome)sev   -0.1825  0.0923  -1.9782  0.0479  -0.3634  -0.0017    * 
# 
# ---
# Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# refit the model with cvvc="varcov" to get the var-cov matrix of the elements in G
res <- update(res, cvvc="varcov")
# regression of the log(OR) for severe exacerbations on the log(OR) for exacerbations
sav <- matreg(R=res$G, cov=TRUE, y=2, x=1, V=res$vvc, means=coef(res))
# set up the plot
par(mar=c(5,4,2,2))
plot(NA, xlim=c(-1.7,1), ylim=c(-3.6,4), bty="l", xlab="Log(OR) for Exacerbations",
     ylab="Log(OR) for Severe Exacerbations")
# add the confidence ellipses for the studies
invisible(sapply(1:nrow(dat), function(i) {
   x <- dat[i,]
   xy <- ellipse(matrix(c(x$v1i,x$cov12i,x$cov12i,x$v2i), nrow=2),
                 centre=c(x$y1i,x$y2i), level=0.5)
   lines(xy[,1],xy[,2], col="gray90")
}))
# add the points for the studies
points(dat$y1i, dat$y2i, pch=21, col="gray80", bg="gray90", cex=1)
# add the regression line
abline(a=coef(sav)[1], b=coef(sav)[2], lwd=3, col="gray60", lty="dashed")
# add the 95% PI ellipsis based on the model
xy <- ellipse(res$G, centre=coef(res), level=0.95)
lines(xy[,1],xy[,2], col="gray30", lwd=3, lty="dotted")
# add the 95% CI ellipsis based on the model
xy <- ellipse(vcov(res), centre=coef(res), level=0.95)
lines(xy[,1],xy[,2], col="gray30", lwd=3)
# add the point for the pooled effects
points(coef(res)[1], coef(res)[2], pch=21, bg="gray40", cex=2)
# add the legend
legend("topright", inset=0.01, lty=c("solid","dotted","dashed"), lwd=3,
       col=c("gray30","gray30","gray60"),
       legend=c("95% Confidence Ellipse", "95% Prediction Interval Ellipse", "Latent Regression Model"))

# test if the mean effect size differs for the two outcomes
anova(res, X=c(1,-1))
# Hypothesis:                                               
# 1: factor(outcome)exa - factor(outcome)sev = 0 
# 
# Results:
#    estimate     se    zval   pval 
# 1:  -0.0486 0.0847 -0.5731 0.5666 
# 
# Test of Hypothesis:
# QM(df = 1) = 0.3285, p-val = 0.5666

 

Conclusion

  • metafor provides a comprehensive collection of functions for conducting meta-analyses in R
  • covers standard methods, but also advanced models
  • for further details, see the package website:
    https://www.metafor-project.org

 

References

Berkey, C. S., Hoaglin, D. C., Antczak-Bouckoms, A., Mosteller, F., & Colditz, G. A. (1998). Meta-analysis of multiple outcomes by regression with random effects. Statistics in Medicine, 17(22), 2537–2550. https://doi.org/10.1002/(sici)1097-0258(19981130)17:22<2537::aid-sim953>3.0.co;2-c
Kalaian, H. A., & Raudenbush, S. W. (1996). A multivariate mixed linear model for meta-analysis. Psychological Methods, 1(3), 227–235. https://doi.org/10.1037/1082-989X.1.3.227
Karner, C., Chong, J., & Poole, P. (2014). Tiotropium versus placebo for chronic obstructive pulmonary disease. Cochrane Database of Systematic Reviews, 7, CD009285. https://doi.org/10.1002/14651858.cd009285.pub3
Trikalinos, T. A., & Olkin, I. (2012). Meta-analysis of effect sizes reported at multiple time points: A multivariate approach. Clinical Trials, 9(5), 610–620. https://doi.org/10.1177/1740774512453218
Viechtbauer, W. (2010). Conducting meta-analyses in R with the metafor package. Journal of Statistical Software, 36(3), 1–48. https://doi.org/10.18637/jss.v036.i03