The metafor Package

The main workhorse of the metafor package is the rma.uni() function (same as rma()), which fits meta-analytic equal-, fixed-, random-, and mixed-effects models using what is often referred to as the "inverse-variance method" (for random/mixed-effects models, this is also called the "normal-normal model", as it assumes normally distributed sampling distributions with known sampling variances and normally distributed random effects for the total/residual amount of heterogeneity).

For multilevel, multivariate, and network meta-analyses, the package also provides the rma.mv() function, which allows for correlated sampling errors and/or true effects. Of course, one can also use the rma.mv() function to fit the same models as the rma.uni() function, since those models are really just special cases of the more general models handled by the rma.mv() function. This is illustrated below.

Note: Alternative modeling approaches for particular types of data (e.g., 2×2 table data, two-group person-time data, proportions, incidence rates) are also available, including the Mantel-Haenszel method (rma.mh() function), Peto's (one-step) method (rma.peto() function), and a variety of suitable mixed-effects (conditional) logistic and Poisson regression models (rma.glmm() function).

To illustrate the use of the rma.uni() and rma.mv() functions for fitting equal-effects models, let's use the BCG vaccine data. We can compute the log risk ratios and corresponding sampling variances with:

library(metafor)
dat <- escalc(measure="RR", ai=tpos, bi=tneg, ci=cpos, di=cneg, data=dat.bcg)

   trial               author year tpos  tneg cpos  cneg ablat      alloc      yi     vi
1      1              Aronson 1948    4   119   11   128    44     random -0.8893 0.3256
2      2     Ferguson & Simes 1949    6   300   29   274    55     random -1.5854 0.1946
3      3      Rosenthal et al 1960    3   228   11   209    42     random -1.3481 0.4154
4      4    Hart & Sutherland 1977   62 13536  248 12619    52     random -1.4416 0.0200
5      5 Frimodt-Moller et al 1973   33  5036   47  5761    13  alternate -0.2175 0.0512
6      6      Stein & Aronson 1953  180  1361  372  1079    44  alternate -0.7861 0.0069
7      7     Vandiviere et al 1973    8  2537   10   619    19     random -1.6209 0.2230
8      8           TPT Madras 1980  505 87886  499 87892    13     random  0.0120 0.0040
9      9     Coetzee & Berjak 1968   29  7470   45  7232    27     random -0.4694 0.0564
10    10      Rosenthal et al 1961   17  1699   65  1600    42 systematic -1.3713 0.0730
11    11       Comstock et al 1974  186 50448  141 27197    18 systematic -0.3394 0.0124
12    12   Comstock & Webster 1969    5  2493    3  2338    33 systematic  0.4459 0.5325
13    13       Comstock et al 1976   27 16886   29 17825    33 systematic -0.0173 0.0714

An equal-effects model can be fitted to these data using the rma.uni() function with:

rma(yi, vi, method="EE", data=dat)

Equal-Effects Model (k = 13)
 
I^2 (total heterogeneity / total variability):   92.12%
H^2 (total variability / sampling variability):  12.69
 
Test for Heterogeneity:
Q(df = 12) = 152.2330, p-val < .0001
 
Model Results:
 
estimate      se      zval    pval    ci.lb    ci.ub
 -0.4303  0.0405  -10.6247  <.0001  -0.5097  -0.3509  ***

Again, using rma() is the same as invoking the rma.uni() function, but the former is quicker to type.

The same model can be fitted to these data using the rma.mv() function with:

rma.mv(yi, vi, data=dat)

Multivariate Meta-Analysis Model (k = 13; method: REML)
 
Variance Components: none
 
Test for Heterogeneity:
Q(df = 12) = 152.2330, p-val < .0001
 
Model Results:
 
estimate       se     zval     pval    ci.lb    ci.ub
 -0.4303   0.0405 -10.6247   <.0001  -0.5097  -0.3509      ***
 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

The results are exactly the same.

In both functions, predictors/covariates/moderators can be included in the model via the mods argument. For example:

rma(yi, vi, mods = ~ ablat + year, method="FE", data=dat)
rma.mv(yi, vi, mods = ~ ablat + year, data=dat)

Again, the results are the same (output not shown).

A random-effects model can be fitted to these data using the rma.uni() function with:

rma(yi, vi, data=dat)

Random-Effects Model (k = 13; tau^2 estimator: REML)
 
tau^2 (estimated amount of total heterogeneity): 0.3132 (SE = 0.1664)
tau (square root of estimated tau^2 value):      0.5597
I^2 (total heterogeneity / total variability):   92.22%
H^2 (total variability / sampling variability):  12.86
 
Test for Heterogeneity:
Q(df = 12) = 152.2330, p-val < .0001
 
Model Results:
 
estimate       se     zval     pval    ci.lb    ci.ub
 -0.7145   0.1798  -3.9744   <.0001  -1.0669  -0.3622      ***
 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

The default for the method argument is to use REML estimation for the amount of heterogeneity and a random-effects model is then automatically fitted.

When using the rma.mv() function, random effects must be explicitly added to the model via the random argument. For a standard random-effects model, we need to add random effects for the trials, which can be done with:

rma.mv(yi, vi, random = ~ 1 | trial, data=dat)

Multivariate Meta-Analysis Model (k = 13; method: REML)
 
Variance Components:
 
            estim    sqrt  nlvls  fixed  factor
sigma^2    0.3132  0.5597     13     no   trial
 
Test for Heterogeneity:
Q(df = 12) = 152.2330, p-val < .0001
 
Model Results:
 
estimate       se     zval     pval    ci.lb    ci.ub
 -0.7145   0.1798  -3.9745   <.0001  -1.0669  -0.3622      ***
 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

These are the same results as obtained earlier.¹⁾

In essence, the meta-analytic random-effects model can be conceptualized as a multilevel model with the true effects at level 2 and the observed effects at level 1. Using typical multilevel model terminology, the random = ~ 1 | trial argument adds random intercepts at level 2 to the model. Note that the variance of the true effects (commonly denoted as $\tau^2$ in the meta-analytic literature) is denoted as $\sigma^2$ in the output above (i.e., the value of $\sigma^2$ given in the output is the estimated variance of the random intercepts).

Again, predictors/covariates/moderators can be included in the model via the mods argument. For example:

rma(yi, vi, mods = ~ ablat + year, data=dat)
rma.mv(yi, vi, mods = ~ ablat + year, random = ~ 1 | trial, data=dat)

The results are again the same (output not shown).