The metafor Package

Berkey et al. (1995) describe the meta-analytic random- and mixed-effects models and provide the equation for the empirical Bayes estimator for the amount of (residual) heterogeneity (p. 398). The models and methods are illustrated with the BCG vaccine dataset (Colditz et al., 1994). The data are provided in Table 1 in the article and can be loaded with:

library(metafor)
dat.bcg

The contents of the dataset are:

   trial               author year tpos  tneg cpos  cneg ablat      alloc
1      1              Aronson 1948    4   119   11   128    44     random
2      2     Ferguson & Simes 1949    6   300   29   274    55     random
3      3      Rosenthal et al 1960    3   228   11   209    42     random
4      4    Hart & Sutherland 1977   62 13536  248 12619    52     random
5      5 Frimodt-Moller et al 1973   33  5036   47  5761    13  alternate
6      6      Stein & Aronson 1953  180  1361  372  1079    44  alternate
7      7     Vandiviere et al 1973    8  2537   10   619    19     random
8      8           TPT Madras 1980  505 87886  499 87892    13     random
9      9     Coetzee & Berjak 1968   29  7470   45  7232    27     random
10    10      Rosenthal et al 1961   17  1699   65  1600    42 systematic
11    11       Comstock et al 1974  186 50448  141 27197    18 systematic
12    12   Comstock & Webster 1969    5  2493    3  2338    33 systematic
13    13       Comstock et al 1976   27 16886   29 17825    33 systematic

Next, we can calculate the log risk ratios and corresponding sampling variances with:

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

The metafor package calculates the sampling variances of the log risk ratios as described by Berkey et al. (1995) on page 399 (right under the 2×2 table). However, note that Berkey et al. (1995) actually use 'smoothed' estimates of the sampling variances for their analyses (to reduce the correlation between the log risk ratios and corresponding sampling variances). We can compute these estimates with:

dat$vi <- with(dat, mean(tneg/tpos)/(tneg+tpos) +
                    mean(cneg/cpos)/(cneg+cpos))

Now we can fit a random-effects model to these data (using the empirical Bayes estimator for the amount of heterogeneity) with:

res.RE <- rma(yi, vi, data=dat, method="EB")
res.RE

Random-Effects Model (k = 13; tau^2 estimator: EB)
 
tau^2 (estimated amount of total heterogeneity): 0.2682 (SE = 0.1801)
tau (square root of estimated tau^2 value):      0.5178
I^2 (total heterogeneity / total variability):   87.49%
H^2 (total variability / sampling variability):  7.99
 
Test for Heterogeneity:
Q(df = 12) = 85.8625, p-val < .0001
 
Model Results:
 
estimate      se     zval    pval    ci.lb    ci.ub
 -0.5429  0.1842  -2.9474  0.0032  -0.9040  -0.1819  **
 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

These results match exactly what Berkey et al. (1995) report on page 408: The amount of heterogeneity ('between-trial variance') is estimated to be $\hat{\tau}^2 = 0.268$ and the pooled estimate is $\hat{\mu} = -0.5429$ with a standard error of $SE[\hat{\mu}] = 0.1842$.

A mixed-effects meta-regression model with absolute latitude as moderators (centered at 33.46 degrees) can be fitted with:

res.ME <- rma(yi, vi, mods=~I(ablat-33.46), data=dat, method="EB")
res.ME

Mixed-Effects Model (k = 13; tau^2 estimator: EB)
 
tau^2 (estimated amount of residual heterogeneity):     0.1572 (SE = 0.1262)
tau (square root of estimated tau^2 value):             0.3965
I^2 (residual heterogeneity / unaccounted variability): 77.42%
H^2 (unaccounted variability / sampling variability):   4.43
R^2 (amount of heterogeneity accounted for):            41.38%
 
Test for Residual Heterogeneity:
QE(df = 11) = 35.8827, p-val = 0.0002
 
Test of Moderators (coefficient 2):
QM(df = 1) = 5.9018, p-val = 0.0151
 
Model Results:
 
                  estimate      se     zval    pval    ci.lb    ci.ub
intrcpt            -0.6303  0.1591  -3.9630  <.0001  -0.9421  -0.3186  ***
I(ablat - 33.46)   -0.0268  0.0110  -2.4294  0.0151  -0.0484  -0.0052    *
 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Again, the results match the findings from Berkey et al. (1995): The residual amount of heterogeneity is now $\hat{\tau}^2 = 0.157$ and the estimated model is $log(RR) = -0.6303 - 0.0268 (x - 33.46)$, where $x$ is the distance from the equator (in degrees latitude). The standard errors of the model coefficients are $SE[b_0] = 0.1591$ and $SE[b_1] = 0.0110$.

The amount of variance (heterogeneity) accounted for by the absolute latitude moderator is provided in the output above. It can also be obtained with:¹⁾

anova(res.RE, res.ME)

        df     AIC     BIC    AICc   logLik    LRT   pval      QE  tau^2      R^2
Full     3 29.3011 30.9960 31.9678 -11.6506               35.8827 0.1572
Reduced  2 32.2578 33.3877 33.4578 -14.1289 4.9566 0.0260 85.8625 0.2682 41.3844%

The value below R^2 indicates that approximately 41% of the heterogeneity has been accounted for.

The predicted average risk ratios at 33.46 and 42 degrees reported by Berkey et al. (1995) can be computed with:

predict(res.ME, newmods=c(33.46,42)-33.46, transf=exp, digits=2)

  pred ci.lb ci.ub pi.lb pi.ub
1 0.53  0.39  0.73  0.23  1.23
2 0.42  0.28  0.64  0.18  1.02

The results from a fixed-effects meta-regression model with absolute latitude as the predictor can be obtained with:

res.FE <- rma(yi, vi, mods=~I(ablat-33.46), data=dat, method="FE")
res.FE

Fixed-Effects with Moderators Model (k = 13)
 
I^2 (residual heterogeneity / unaccounted variability): 69.34%
H^2 (unaccounted variability / sampling variability):   3.26
 
Test for Residual Heterogeneity:
QE(df = 11) = 35.8827, p-val = 0.0002
 
Test of Moderators (coefficient 2):
QM(df = 1) = 49.9798, p-val < .0001
 
Model Results:
 
                  estimate      se     zval    pval    ci.lb    ci.ub
intrcpt            -0.5949  0.0696  -8.5462  <.0001  -0.7314  -0.4585  ***
I(ablat - 33.46)   -0.0282  0.0040  -7.0696  <.0001  -0.0360  -0.0204  ***
 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

The predicted risk ratios at 33.46 and 42 degrees are now:

predict(res.FE, newmods=c(33.46,42)-33.46, transf=exp, digits=2)

  pred ci.lb ci.ub
1 0.55  0.48  0.63
2 0.43  0.36  0.52

Berkey, C. S., Hoaglin, D. C., Mosteller, F., & Colditz, G. A. (1995). A random-effects regression model for meta-analysis. Statistics in Medicine, 14(4), 395–411.

Colditz, G. A., Brewer, T. F., Berkey, C. S., Wilson, M. E., Burdick, E., Fineberg, H. V., et al. (1994). Efficacy of BCG vaccine in the prevention of tuberculosis: Meta-analysis of the published literature. Journal of the American Medical Association, 271(9), 698–702.