This vignette provides up-to-date commands for the analyses in “How to perform a meta-analysis with R: a practical tutorial”, Evid Based Ment Health (Balduzzi, Rücker, and Schwarzer 2019).

Install R packages

install.packages(c("meta", "metasens"))

Make R packages available

library(meta)
#> Loading required package: metadat
#> Loading 'meta' package (version 8.0-2).
#> Type 'help(meta)' for a brief overview.

library(metasens)

Note, a similar message would be printed for R package metasens. However, this vignette does not actually load metasens as it might not be installed in addition to meta.

Default settings for R session

Print results with two significant digits and use Paule-Mandel estimator for between-study variance.

settings.meta(digits = 2, method.tau = "PM")

Note, in the publication, argument ‘method.tau’ was used in R function metabin(). Here, we set the Paule-Mandel method as the default for any meta-analysis conducted in the current R session.

Import the dataset

joy = read.csv("Joy2006.txt")
# Add new variable: miss
joy$miss = ifelse((joy$drop.h + joy$drop.p) == 0, 
  "Without missing data", "With missing data")
head(joy)
#>      author year resp.h fail.h drop.h resp.p fail.p drop.p                 miss
#> 1 Arvanitis 1997     25     25      2     18     33      0    With missing data
#> 2   Beasley 1996     29     18     22     20     14     34    With missing data
#> 3  Bechelli 1983     12     17      1      2     28      1    With missing data
#> 4   Borison 1992      3      9      0      0     12      0 Without missing data
#> 5 Chouinard 1993     10     11      0      3     19      0 Without missing data
#> 6    Durost 1964     11      8      0      1     14      0 Without missing data
str(joy)
#> 'data.frame':    17 obs. of  9 variables:
#>  $ author: chr  "Arvanitis" "Beasley" "Bechelli" "Borison" ...
#>  $ year  : int  1997 1996 1983 1992 1993 1964 1962 1974 1994 1982 ...
#>  $ resp.h: int  25 29 12 3 10 11 7 8 19 1 ...
#>  $ fail.h: int  25 18 17 9 11 8 18 9 45 9 ...
#>  $ drop.h: int  2 22 1 0 0 0 1 0 2 0 ...
#>  $ resp.p: int  18 20 2 0 3 1 4 3 14 0 ...
#>  $ fail.p: int  33 14 28 12 19 14 21 10 50 10 ...
#>  $ drop.p: int  0 34 1 0 0 0 1 0 2 0 ...
#>  $ miss  : chr  "With missing data" "With missing data" "With missing data" "Without missing data" ...

Section ‘Fixed effect and random effects meta-analysis’

m.publ = metabin(resp.h, resp.h + fail.h, resp.p, resp.p + fail.p,
  data = joy, studlab = paste0(author, " (", year, ")"),
  label.e = "Haloperidol", label.c = "Placebo",
  label.left = "Favours placebo", label.right = "Favours haloperidol")

Print results of meta-analysis (Figure 1).

summary(m.publ)
#>                         RR         95%-CI %W(common) %W(random)
#> Arvanitis (1997)      1.42 [0.89;   2.25]       21.1       14.1
#> Beasley (1996)        1.05 [0.73;   1.50]       27.5       15.6
#> Bechelli (1983)       6.21 [1.52;  25.35]        2.3        4.7
#> Borison (1992)        7.00 [0.40; 121.94]        0.6        1.4
#> Chouinard (1993)      3.49 [1.11;  10.95]        3.5        6.3
#> Durost (1964)         8.68 [1.26;  59.95]        1.3        2.8
#> Garry (1962)          1.75 [0.58;   5.24]        4.7        6.7
#> Howard (1974)         2.04 [0.67;   6.21]        4.0        6.6
#> Marder (1994)         1.36 [0.75;   2.47]       16.6       12.2
#> Nishikawa (1982)      3.00 [0.14;  65.55]        0.6        1.2
#> Nishikawa (1984)      9.00 [0.57; 142.29]        0.8        1.5
#> Reschke (1974)        3.79 [1.06;  13.60]        3.4        5.4
#> Selman (1976)         1.48 [0.94;   2.35]       10.3       14.1
#> Serafetinides (1972)  8.38 [0.50; 141.44]        0.6        1.4
#> Simpson (1967)        2.27 [0.12;  41.77]        0.8        1.4
#> Spencer (1992)       11.00 [1.67;  72.40]        1.2        3.0
#> Vichaiya (1971)      19.00 [1.16; 311.71]        0.6        1.5
#> 
#> Number of studies: k = 17
#> Number of observations: o = 818 (o.e = 446, o.c = 372)
#> Number of events: e = 274
#> 
#>                        RR       95%-CI    z  p-value
#> Common effect model  2.09 [1.69; 2.59] 6.71 < 0.0001
#> Random effects model 2.15 [1.51; 3.06] 4.23 < 0.0001
#> 
#> Quantifying heterogeneity (with 95%-CIs):
#>  tau^2 = 0.1754 [0.0000; 1.0088]; tau = 0.4188 [0.0000; 1.0044]
#>  I^2 = 41.3% [0.0%; 67.0%]; H = 1.30 [1.00; 1.74]
#> 
#> Test of heterogeneity:
#>      Q d.f. p-value
#>  27.24   16  0.0388
#> 
#> Details of meta-analysis methods:
#> - Mantel-Haenszel method (common effect model)
#> - Inverse variance method (random effects model)
#> - Paule-Mandel estimator for tau^2
#> - Q-Profile method for confidence interval of tau^2 and tau
#> - Calculation of I^2 based on Q
#> - Continuity correction of 0.5 in studies with zero cell frequencies

Same printout (result not shown)

print(summary(m.publ))

Create Figure 2 (file ‘figure2.pdf’).

forest(m.publ, sortvar = year, prediction = TRUE,
  file = "figure2.pdf", width = 10)

Section ‘Assessing the impact of missing outcome data’

Subgroup analysis of studies with and without missing data

m.publ.sub = update(m.publ, subgroup = miss, print.subgroup.name = FALSE)
m.publ.sub
#> Number of studies: k = 17
#> Number of observations: o = 818 (o.e = 446, o.c = 372)
#> Number of events: e = 274
#> 
#>                        RR       95%-CI    z  p-value
#> Common effect model  2.09 [1.69; 2.59] 6.71 < 0.0001
#> Random effects model 2.15 [1.51; 3.06] 4.23 < 0.0001
#> 
#> Quantifying heterogeneity (with 95%-CIs):
#>  tau^2 = 0.1754 [0.0000; 1.0088]; tau = 0.4188 [0.0000; 1.0044]
#>  I^2 = 41.3% [0.0%; 67.0%]; H = 1.30 [1.00; 1.74]
#> 
#> Test of heterogeneity:
#>      Q d.f. p-value
#>  27.24   16  0.0388
#> 
#> Results for subgroups (common effect model):
#>                        k   RR       95%-CI     Q   I^2
#> With missing data     10 1.70 [1.35; 2.14] 13.73 34.4%
#> Without missing data   7 4.36 [2.45; 7.77]  3.35  0.0%
#> 
#> Test for subgroup differences (common effect model):
#>                   Q d.f. p-value
#> Between groups 8.82    1  0.0030
#> 
#> Results for subgroups (random effects model):
#>                        k   RR       95%-CI  tau^2    tau
#> With missing data     10 1.71 [1.13; 2.59] 0.1785 0.4224
#> Without missing data   7 3.80 [2.13; 6.80]      0      0
#> 
#> Test for subgroup differences (random effects model):
#>                   Q d.f. p-value
#> Between groups 4.80    1  0.0285
#> 
#> Details of meta-analysis methods:
#> - Mantel-Haenszel method (common effect model)
#> - Inverse variance method (random effects model)
#> - Paule-Mandel estimator for tau^2
#> - Q-Profile method for confidence interval of tau^2 and tau
#> - Calculation of I^2 based on Q
#> - Continuity correction of 0.5 in studies with zero cell frequencies

Create Figure 3 (file ‘figure3.pdf’).

forest(m.publ.sub, sortvar = year,
  xlim = c(0.1, 100), at = c(0.1, 0.3, 1, 3, 10, 30, 100),
  test.subgroup.common = FALSE,
  label.test.subgroup.random = "Test for subgroup differences:",
  file = "figure3.pdf", width = 10)

Use imputation methods

# Impute as no events (ICA-0) - default
mmiss.0 = metamiss(m.publ, drop.h, drop.p)
# Impute as events (ICA-1)
mmiss.1 = metamiss(m.publ, drop.h, drop.p, method = "1")
# Observed risk in control group (ICA-pc)
mmiss.pc = metamiss(m.publ, drop.h, drop.p, method = "pc")
# Observed risk in experimental group (ICA-pe)
mmiss.pe = metamiss(m.publ, drop.h, drop.p, method = "pe")
# Observed group-specific risks (ICA-p)
mmiss.p = metamiss(m.publ, drop.h, drop.p, method = "p")
# Best-case scenario (ICA-b)
mmiss.b = metamiss(m.publ, drop.h, drop.p, method = "b", small.values = "bad")
# Worst-case scenario (ICA-w)
mmiss.w = metamiss(m.publ, drop.h, drop.p, method = "w", small.values = "bad")
# Gamble-Hollis method
mmiss.gh = metamiss(m.publ, drop.h, drop.p, method = "GH")
# IMOR.e = 2 and IMOR.c = 2 (same as available case analysis)
mmiss.imor2 = metamiss(m.publ, drop.h, drop.p, method = "IMOR", IMOR.e = 2)
# IMOR.e = 0.5 and IMOR.c = 0.5
mmiss.imor0.5 = metamiss(m.publ, drop.h, drop.p, method = "IMOR", IMOR.e = 0.5)

Summarise results using R function metabind().

meths = c("Available case analysis (ACA)",
  "Impute no events (ICA-0)", "Impute events (ICA-1)",
  "Observed risk in control group (ICA-pc)",
  "Observed risk in experimental group (ICA-pe)",
  "Observed group-specific risks (ICA-p)",
  "Best-case scenario (ICA-b)", "Worst-case scenario (ICA-w)",
  "Gamble-Hollis analysis",
  "IMOR.e = 2, IMOR.c = 2", "IMOR.e = 0.5, IMOR.c = 0.5")
# Use inverse-variance method for pooling (which is used for
# imputation methods)
m.publ.iv = update(m.publ, method = "Inverse")
# Combine results (random effects)
mbr = metabind(m.publ.iv,
  mmiss.0, mmiss.1,
  mmiss.pc, mmiss.pe, mmiss.p,
  mmiss.b, mmiss.w, mmiss.gh,
  mmiss.imor2, mmiss.imor0.5,
  name = meths, pooled = "random")

Create Figure 4 (file ‘figure4.pdf’).

forest(mbr, xlim = c(0.5, 4),
  leftcols = c("studlab", "I2", "tau2", "Q", "pval.Q"),
  leftlab = c("Meta-Analysis Method", "I2", "Tau2", "Q", "P-value"),
  type = "diamond",
  digits.addcols = c(4, 2, 2, 2), just.addcols = "right",
  file = "figure4.pdf", width = 10)

Section ‘Assessing and accounting for small-study effects’

Funnel plot

funnel(m.publ)

Harbord’s score test for funnel plot asymmetry

metabias(m.publ, method.bias = "score")
#> Linear regression test of funnel plot asymmetry
#> 
#> Test result: t = 4.56, df = 15, p-value = 0.0004
#> Bias estimate: 2.21 (SE = 0.4853)
#> 
#> Details:
#> - multiplicative residual heterogeneity variance (tau^2 = 1.0948)
#> - predictor: standard error of score
#> - weight:    inverse variance of score
#> - reference: Harbord et al. (2006), Stat Med

Trim-and-fill method

tf.publ = trimfill(m.publ)
tf.publ
#> Number of studies: k = 26 (with 9 added studies)
#> Number of observations: o = 1174 (o.e = 645, o.c = 529)
#> Number of events: e = 374
#> 
#>                        RR       95%-CI    z p-value
#> Random effects model 1.40 [0.83; 2.38] 1.26  0.2063
#> 
#> Quantifying heterogeneity (with 95%-CIs):
#>  tau^2 = 1.0983 [0.2929; 3.1894]; tau = 1.0480 [0.5412; 1.7859]
#>  I^2 = 56.2% [32.1%; 71.8%]; H = 1.51 [1.21; 1.88]
#> 
#> Test of heterogeneity:
#>      Q d.f. p-value
#>  57.13   25  0.0003
#> 
#> Details of meta-analysis methods:
#> - Inverse variance method
#> - Paule-Mandel estimator for tau^2
#> - Q-Profile method for confidence interval of tau^2 and tau
#> - Calculation of I^2 based on Q
#> - Trim-and-fill method to adjust for funnel plot asymmetry (L-estimator)

summary(tf.publ)
#>                                 RR         95%-CI %W(random)
#> Arvanitis (1997)              1.42 [0.89;   2.25]        6.2
#> Beasley (1996)                1.05 [0.73;   1.50]        6.4
#> Bechelli (1983)               6.21 [1.52;  25.35]        4.5
#> Borison (1992)                7.00 [0.40; 121.94]        2.2
#> Chouinard (1993)              3.49 [1.11;  10.95]        5.0
#> Durost (1964)                 8.68 [1.26;  59.95]        3.5
#> Garry (1962)                  1.75 [0.58;   5.24]        5.1
#> Howard (1974)                 2.04 [0.67;   6.21]        5.1
#> Marder (1994)                 1.36 [0.75;   2.47]        6.0
#> Nishikawa (1982)              3.00 [0.14;  65.55]        2.0
#> Nishikawa (1984)              9.00 [0.57; 142.29]        2.3
#> Reschke (1974)                3.79 [1.06;  13.60]        4.7
#> Selman (1976)                 1.48 [0.94;   2.35]        6.2
#> Serafetinides (1972)          8.38 [0.50; 141.44]        2.3
#> Simpson (1967)                2.27 [0.12;  41.77]        2.2
#> Spencer (1992)               11.00 [1.67;  72.40]        3.6
#> Vichaiya (1971)              19.00 [1.16; 311.71]        2.3
#> Filled: Chouinard (1993)      0.50 [0.16;   1.55]        5.0
#> Filled: Reschke (1974)        0.46 [0.13;   1.64]        4.7
#> Filled: Bechelli (1983)       0.28 [0.07;   1.14]        4.5
#> Filled: Borison (1992)        0.25 [0.01;   4.31]        2.2
#> Filled: Serafetinides (1972)  0.21 [0.01;   3.49]        2.3
#> Filled: Durost (1964)         0.20 [0.03;   1.38]        3.5
#> Filled: Nishikawa (1984)      0.19 [0.01;   3.04]        2.3
#> Filled: Spencer (1992)        0.16 [0.02;   1.04]        3.6
#> Filled: Vichaiya (1971)       0.09 [0.01;   1.49]        2.3
#> 
#> Number of studies: k = 26 (with 9 added studies)
#> Number of observations: o = 1174 (o.e = 645, o.c = 529)
#> Number of events: e = 374
#> 
#>                        RR       95%-CI    z p-value
#> Random effects model 1.40 [0.83; 2.38] 1.26  0.2063
#> 
#> Quantifying heterogeneity (with 95%-CIs):
#>  tau^2 = 1.0983 [0.2929; 3.1894]; tau = 1.0480 [0.5412; 1.7859]
#>  I^2 = 56.2% [32.1%; 71.8%]; H = 1.51 [1.21; 1.88]
#> 
#> Test of heterogeneity:
#>      Q d.f. p-value
#>  57.13   25  0.0003
#> 
#> Details of meta-analysis methods:
#> - Inverse variance method
#> - Paule-Mandel estimator for tau^2
#> - Q-Profile method for confidence interval of tau^2 and tau
#> - Calculation of I^2 based on Q
#> - Trim-and-fill method to adjust for funnel plot asymmetry (L-estimator)

funnel(tf.publ)

Limit meta-analysis

l1.publ = limitmeta(m.publ)

Note, the printout for the limit meta-analysis is not shown in this vignette as the installation of R package metasens is optional.

l1.publ

Create Figure 5 (file ‘figure5.pdf’).

pdf("figure5.pdf", width = 10, height = 10)
#
par(mfrow = c(2, 2), pty = "s",
    oma = c(0, 0, 0, 0), mar = c(4.1, 3.1, 2.1, 1.1))
#
funnel(m.publ, xlim = c(0.05, 50), axes = FALSE)
axis(1, at = c(0.1, 0.2, 0.5, 1, 2, 5, 10, 50))
axis(2, at = c(0, 0.5, 1, 1.5))
box()
title(main = "Panel A: Funnel plot", adj = 0)
#
funnel(m.publ, xlim = c(0.05, 50), axes = FALSE,
       contour.levels = c(0.9, 0.95, 0.99),
       col.contour = c("darkgray", "gray", "lightgray"))
legend("topright",
       c("p < 1%", "1% < p < 5%", "5% < p < 10%", "p > 10%"),
       fill = c("lightgray", "gray", "darkgray", "white"),
       border = "white", bg = "white")
axis(1, at = c(0.1, 0.2, 0.5, 1, 2, 5, 10, 50))
axis(2, at = c(0, 0.5, 1, 1.5))
box()
title(main = "Panel B: Contour-enhanced funnel plot", adj = 0)
#
funnel(tf.publ, xlim = c(0.05, 50), axes = FALSE)
axis(1, at = c(0.1, 0.2, 0.5, 1, 2, 5, 10, 50))
axis(2, at = c(0, 0.5, 1, 1.5))
box()
title(main = "Panel C: Trim-and-fill method", adj = 0)
#
funnel(l1.publ, xlim = c(0.05, 50), axes = FALSE,
       col.line = 8, lwd.line = 3)
axis(1, at = c(0.1, 0.2, 0.5, 1, 2, 5, 10, 50))
axis(2, at = c(0, 0.5, 1, 1.5))
box()
title(main = "Panel D: Limit meta-analysis", adj = 0)
#
dev.off()

References

Balduzzi, Sara, Gerta Rücker, and Guido Schwarzer. 2019. “How to Perform a Meta-Analysis with R: A Practical Tutorial.” Evidence-Based Mental Health 22 (4): 153–60. doi:10.1136/ebmental-2019-300117.