# The metafor Package

A Meta-Analysis Package for R

### Sidebar

analyses:stijnen2010

## Stijnen et al. (2010)

### The Methods and Data

Dichotomous and event count data (based on which one can calculate effect size or outcome measures, such as odds ratios, incidence rate ratios, proportions, and incidence rates) are often assumed to arise from binomial and Poisson distributed data. Stijnen et al. (2010) describe meta-analytic models that are directly based on these distibutions. The collection of models are essentially special cases of generalized linear (mixed-effects) models (i.e., mixed-effects logistic and Poisson regression models). For 2×2 table data, a mixed-effects conditional logistic model (based on the non-central hypergeometric distribution) can also be used. Mixed-effects logistic regression models for dichotomous data are often referred to as "binomial-normal models" in the meta-analytic literature. Analogously, mixed-effects Poisson regression models for event count data could be referred to as "Poisson-normal models". The results obtained from such models can be contrasted with those obtained from the standard (inverse-variance) approach (sometimes called the "normal-normal model"). The various models described in the article are illustrated with datasets from two meta-analyses comparing the risk of catheter-related bloodstream infection (CRBSI) when using anti-infective-treated versus standard catheters (Niel-Weise et al., 2007, 2008).

### Binomial-Normal Model for the Meta-Analysis of Proportions

Niel-Weise et al. (2007) conducted a meta-analysis comparing the CRBSI risk when using anti-infective-treated versus standard catheters in the acute care setting. These data can be used to illustrate the binomial-normal model for the meta-analysis of proportions (i.e., infection risks) within the standard catheter and anti-infective catheter groups. The data can be extracted from Figure 2 (p. 2063) in Niel-Weise et al. (2007) and can be loaded with:

library(metafor)
dat <- dat.nielweise2007
dat

(I copy the dataset into 'dat', which is a bit shorter and therefore easier to type further below). The contents of the dataset are:

   study         author year ai n1i ci n2i
1      1           Bach 1996  0 116  3 117
2      2         George 1997  1  44  3  35
3      3           Maki 1997  2 208  9 195
4      4           Raad 1997  0 130  7 136
5      5          Heard 1998  5 151  6 157
6      6         Collin 1999  1  98  4 139
7      7         Hannan 1999  1 174  3 177
8      8          Marik 1999  1  74  2  39
9      9         Pierce 2000  1  97 19 103
10    10          Sheng 2000  1 113  2 122
11    11 Chatzinikolaou 2003  0  66  7  64
12    12         Corral 2003  0  70  1  58
13    13   Brun-Buisson 2004  3 188  5 175
14    14           Leon 2004  6 187 11 180
15    15          Yucel 2004  0 118  0 105
16    16        Moretti 2005  0 252  1 262
17    17           Rupp 2005  1 345  3 362
18    18           Osma 2006  4  64  1  69

Variables n1i and n2i denote the number of patients receiving an anti-infective or standard catheter, respectively, and ai and ci the number of CRBSIs in the respective groups. We can now fit a "normal-normal model" to the data for the two groups separately. Here, we use the logit transformed proportions (i.e., log odds) for the meta-analysis. First, for the standard catheter groups, this can be done with:

res <- rma(measure="PLO", xi=ci, ni=n2i, data=dat)
print(res, digits=3)
Random-Effects Model (k = 18; tau^2 estimator: REML)

tau^2 (estimated amount of total heterogeneity): 0.663 (SE = 0.338)
tau (square root of estimated tau^2 value):      0.814
I^2 (total heterogeneity / total variability):   74.01%
H^2 (total variability / sampling variability):  3.85

Test for Heterogeneity:
Q(df = 17) = 72.166, p-val < .001

Model Results:

estimate       se     zval     pval    ci.lb    ci.ub
-3.302    0.238  -13.883    <.001   -3.768   -2.836      ***

---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

The estimated average log odds (i.e., $\hat{\mu} = -3.302$) can be back-transformed using the inverse logit transformation with:

predict(res, transf=transf.ilogit, digits=3)
  pred ci.lb ci.ub cr.lb cr.ub
0.036 0.023 0.055 0.007 0.163

As pointed out in the article, the proper interpretation of the back-transformed value (i.e., $0.036$) is that it reflects the median infection risk (although it would often be interpreted as the average risk).

Next, for the anti-infective catheter groups, we obtain:

res <- rma(measure="PLO", xi=ai, ni=n1i, data=dat)
print(res, digits=3)
Random-Effects Model (k = 18; tau^2 estimator: REML)

tau^2 (estimated amount of total heterogeneity): 0.393 (SE = 0.374)
tau (square root of estimated tau^2 value):      0.627
I^2 (total heterogeneity / total variability):   37.60%
H^2 (total variability / sampling variability):  1.60

Test for Heterogeneity:
Q(df = 17) = 23.671, p-val = 0.129

Model Results:

estimate       se     zval     pval    ci.lb    ci.ub
-4.260    0.259  -16.457    <.001   -4.768   -3.753      ***

---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

And again, using the back-transformation:

predict(res, transf=transf.ilogit, digits=3)
  pred ci.lb ci.ub cr.lb cr.ub
0.014 0.008 0.023 0.004 0.051

Note that REML estimation is the default for the rma() function. Also, calculation of the log odds (and the corresponding sampling variance) is problematic when a group has zero events, as is the case in some of the studies. The standard continuity correction (adding 1/2 to the number of patients with and without an event) is automatically applied inside the function. These results match what is reported by Stijnen et al. (2010) in Table II (p. 3050).

Instead of the normal-normal model, we can use a "binomial-normal model", which in this case is just a logistic regression model with a random intercept. This model can be fitted to the data for the standard catheter groups with:

res <- rma.glmm(measure="PLO", xi=ci, ni=n2i, data=dat)
print(res, digits=3)
Random-Effects Model (k = 18; tau^2 estimator: ML)

tau^2 (estimated amount of total heterogeneity): 0.812
tau (square root of estimated tau^2 value):      0.901
I^2 (total heterogeneity / total variability):   77.73%
H^2 (total variability / sampling variability):  4.49

Tests for Heterogeneity:
Wld(df = 17) = 69.103, p-val < .001
LRT(df = 17) = 85.284, p-val < .001

Model Results:

estimate       se     zval     pval    ci.lb    ci.ub
-3.496    0.257  -13.605    <.001   -4.000   -2.993      ***

---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
predict(res, transf=transf.ilogit, digits=3)
  pred ci.lb ci.ub cr.lb cr.ub
0.029 0.018 0.048 0.005 0.160

And again, for the anti-infective catheter groups:

res <- rma.glmm(measure="PLO", xi=ai, ni=n1i, data=dat)
print(res, digits=3)
Random-Effects Model (k = 18; tau^2 estimator: ML)

tau^2 (estimated amount of total heterogeneity): 0.827
tau (square root of estimated tau^2 value):      0.909
I^2 (total heterogeneity / total variability):   55.91%
H^2 (total variability / sampling variability):  2.27

Tests for Heterogeneity:
Wld(df = 17) = 16.145, p-val = 0.514
LRT(df = 17) = 37.990, p-val = 0.002

Model Results:

estimate       se     zval     pval    ci.lb    ci.ub
-4.812    0.355  -13.537    <.001   -5.509   -4.115      ***

---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
predict(res, transf=transf.ilogit, digits=3)
  pred ci.lb ci.ub cr.lb cr.ub
0.008 0.004 0.016 0.001 0.052

By default, a random/mixed-effects model using ML estimation is fitted by rma.glmm() as was done in the article. Again, these results match what is reported in the article.

### Hypergeometric-Normal Model for the Meta-Analysis of Odds Ratios

Next, the authors describe a "hypergeometric-normal model" for the meta-analysis of odds ratios. This is essentially a mixed-effects conditional logistic model. The model results from conditioning on the total number of cases within each study, leading to the non-central hypergeometric distribution for the 2×2 tables. Since fitting this model can be difficult and computationally expensive, one can approximate the exact likelihood by a binomial distribution, which works well when the number of cases in each study is small relative to the group sizes (as is the case here). Again, the results obtained from this model can be compared with the normal-normal model approach.

In fact, again, we start with the normal-normal model, using the log odds ratio as the outcome measure. This model can be fitted with:

res <- rma(measure="OR", ai=ai, n1i=n1i, ci=ci, n2i=n2i, data=dat, drop00=TRUE)
print(res, digits=3)
Random-Effects Model (k = 17; tau^2 estimator: REML)

tau^2 (estimated amount of total heterogeneity): 0.036 (SE = 0.299)
tau (square root of estimated tau^2 value):      0.189
I^2 (total heterogeneity / total variability):   3.46%
H^2 (total variability / sampling variability):  1.04

Test for Heterogeneity:
Q(df = 16) = 15.812, p-val = 0.466

Model Results:

estimate       se     zval     pval    ci.lb    ci.ub
-0.980    0.243   -4.027    <.001   -1.458   -0.503      ***

---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

The estimated average log odds ratio (i.e., $\hat{\mu} = -0.980$) can be back-transformed through exponentiation with:

predict(res, transf=exp, digits=2)
 pred ci.lb ci.ub cr.lb cr.ub
0.38  0.23  0.60  0.21  0.69

Note that there is one study with zero events in both arms. Such studies can be explicitly dropped from the analysis by setting drop00=TRUE (careful: this is not the default in the rma() function).

We can fit the hypergeometric-normal model with:

res <- rma.glmm(measure="OR", ai=ai, n1i=n1i, ci=ci, n2i=n2i, data=dat, model="CM.EL")
print(res, digits=3)

Model fitting may take a few moments to complete. The results are:

Random-Effects Model (k = 17; tau^2 estimator: ML)
Model Type: Conditional Model with Exact Likelihood

tau^2 (estimated amount of total heterogeneity): 0.693 (SE = 0.640)
tau (square root of estimated tau^2 value):      0.833
I^2 (total heterogeneity / total variability):   41.13%
H^2 (total variability / sampling variability):  1.70

Tests for Heterogeneity:
Wld(df = 16) = 11.866, p-val = 0.753
LRT(df = 16) = 28.609, p-val = 0.027

Model Results:

estimate       se     zval     pval    ci.lb    ci.ub
-1.353    0.351   -3.855    <.001   -2.041   -0.665      ***

---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

And we can back-transform the estimated log odds ratio with:

predict(res, transf=exp, digits=2)
 pred ci.lb ci.ub cr.lb cr.ub
0.26  0.13  0.51  0.04  1.52

Since the likelihood for studies with zero events is flat, such studies are automatically dropped by the rma.glmm() function (i.e., drop00=TRUE by default for the rma.glmm() function).

Finally, using the approximation to the exact likelihood, we can fit the same model with:

res <- rma.glmm(measure="OR", ai=ai, n1i=n1i, ci=ci, n2i=n2i, data=dat, model="CM.AL")
print(res, digits=3)
Random-Effects Model (k = 17; tau^2 estimator: ML)
Model Type: Conditional Model with Approximate Likelihood

tau^2 (estimated amount of total heterogeneity): 0.601
tau (square root of estimated tau^2 value):      0.775
I^2 (total heterogeneity / total variability):   37.70%
H^2 (total variability / sampling variability):  1.61

Tests for Heterogeneity:
Wld(df = 16) = 11.115, p-val = 0.802
LRT(df = 16) = 27.392, p-val = 0.037

Model Results:

estimate       se     zval     pval    ci.lb    ci.ub
-1.303    0.339   -3.847    <.001   -1.966   -0.639      ***

---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
predict(res, transf=exp, digits=2)
 pred ci.lb ci.ub cr.lb cr.ub
0.27  0.14  0.53  0.05  1.43

All of these results match what is reported in Table III (p. 3052).

### Poisson-Normal Model for the Meta-Analysis of Incidence Rates

Niel-Weise et al. (2008) conducted a meta-analysis comparing the CRBSI risk when using anti-infective-treated versus standard catheters for total parenteral nutrition or chemotherapy. These data can be used to illustrate Poisson-normal models for the meta-analysis of incidence rates and incidence rate ratios. The data can be extracted from Table III (p. 119) in Niel-Weise et al. (2008) and can be loaded with:

dat <- dat.nielweise2008
dat

(I copy the dataset into 'dat', which is a bit shorter and therefore easier to type further below). The contents of the dataset are:

  study          authors year x1i   t1i x2i   t2i
1     1      Bong et al. 2003   7  1344  11  1988
2     2    Ciresi et al. 1996   8  1600   8  1461
3     3     Hanna et al. 2004   3 12012  14 10962
4     4    Harter et al. 2002   6  1536  10  1503
5     5    Jaeger et al. 2001   1   370   1   483
6     6    Jaeger et al. 2005   1   729   8   913
7     7    Logghe et al. 1997  17  6760  15  6840
8     8 Ostendorf et al. 2005   3  1107   7  1015
9     9 Pemberton et al. 1996   2   320   3   440

Here, the number of infections (x1i and x2i) and the total number of catheter days (t1i and t2i) are given for patients in the anti-infective-treated and standard catheter groups, respectively. Note that these data differ very slighty (for a few studies) from the data reported by Stijnen et al. (2010) in the appendix of their article (p. 3062)1) but these differences are inconsequential for the results.

Following the authors, we first divide the total number of catheter days by 1000, so that the estimated (average) incidence rates reflect the expected number of events per 1000 days:

dat$t1i <- dat$t1i/1000
dat$t2i <- dat$t2i/1000

Next, the normal-normal model can be fitted using the log incidence rate as the outcome measure. For patients in the standard catheter groups, this can be done with:

res <- rma(measure="IRLN", xi=x2i, ti=t2i, data=dat)
print(res, digits=3)
Random-Effects Model (k = 9; tau^2 estimator: REML)

tau^2 (estimated amount of total heterogeneity): 0.370 (SE = 0.259)
tau (square root of estimated tau^2 value):      0.608
I^2 (total heterogeneity / total variability):   75.38%
H^2 (total variability / sampling variability):  4.06

Test for Heterogeneity:
Q(df = 8) = 36.384, p-val < .001

Model Results:

estimate       se     zval     pval    ci.lb    ci.ub
1.468    0.243    6.051    <.001    0.992    1.943      ***

---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

and the estimated average log incidence rate (i.e., $\hat{\mu} = 1.468$) can be back-transformed through exponentiation:

predict(res, transf=exp, digits=2)
 pred ci.lb ci.ub cr.lb cr.ub
4.34  2.70  6.98  1.20 15.66

The same model can be fitted to patients in the anti-infective-treated catheter groups with:

res <- rma(measure="IRLN", xi=x1i, ti=t1i, data=dat)
print(res, digits=3)
Random-Effects Model (k = 9; tau^2 estimator: REML)

tau^2 (estimated amount of total heterogeneity): 0.639 (SE = 0.467)
tau (square root of estimated tau^2 value):      0.800
I^2 (total heterogeneity / total variability):   75.41%
H^2 (total variability / sampling variability):  4.07

Test for Heterogeneity:
Q(df = 8) = 25.427, p-val = 0.001

Model Results:

estimate       se     zval     pval    ci.lb    ci.ub
0.981    0.326    3.009    0.003    0.342    1.620       **

---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
predict(res, transf=exp, digits=2)
 pred ci.lb ci.ub cr.lb cr.ub
2.67  1.41  5.05  0.49 14.49

Now we fit the Poisson-normal model to the same data. For patients in the standard catheter groups, this can be done with:

res <- rma.glmm(measure="IRLN", xi=x2i, ti=t2i, data=dat)
print(res, digits=3)
Random-Effects Model (k = 9; tau^2 estimator: ML)

tau^2 (estimated amount of total heterogeneity): 0.317
tau (square root of estimated tau^2 value):      0.563
I^2 (total heterogeneity / total variability):   72.38%
H^2 (total variability / sampling variability):  3.62

Tests for Heterogeneity:
Wld(df = 8) = 36.384, p-val < .001
LRT(df = 8) = 38.330, p-val < .001

Model Results:

estimate       se     zval     pval    ci.lb    ci.ub
1.401    0.231    6.063    <.001    0.948    1.853      ***

---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
predict(res, transf=exp, digits=2)
 pred ci.lb ci.ub cr.lb cr.ub
4.06  2.58  6.38  1.23 13.37

For patients in the anti-infective-treated catheter groups, we obtain:

res <- rma.glmm(measure="IRLN", xi=x1i, ti=t1i, data=dat)
print(res, digits=3)
Random-Effects Model (k = 9; tau^2 estimator: ML)

tau^2 (estimated amount of total heterogeneity): 0.654
tau (square root of estimated tau^2 value):      0.809
I^2 (total heterogeneity / total variability):   75.84%
H^2 (total variability / sampling variability):  4.14

Tests for Heterogeneity:
Wld(df = 8) = 25.427, p-val = 0.001
LRT(df = 8) = 44.488, p-val < .001

Model Results:

estimate       se     zval     pval    ci.lb    ci.ub
0.849    0.330    2.572    0.010    0.202    1.497        *

---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
predict(res, transf=exp, digits=2)
 pred ci.lb ci.ub cr.lb cr.ub
2.34  1.22  4.47  0.42 12.96

These findings match what is reported in Table VII (p. 3054).

### Binomial-Normal Model for the Meta-Analysis of Incidence Rate Ratios

Next, we consider the mixed-effects conditional Poisson regression model described by the authors. Through the conditioning on the total number of events within each study, we actually obtain a binomial-normal model for the data within each study. But first, we start again with the normal-normal model, using the log incidence rate ratio as the outcome measure:

res <- rma(measure="IRR", x1i=x1i, t1i=t1i, x2i=x2i, t2i=t2i, data=dat)
print(res, digits=3)
Random-Effects Model (k = 9; tau^2 estimator: REML)

tau^2 (estimated amount of total heterogeneity): 0.094 (SE = 0.212)
tau (square root of estimated tau^2 value):      0.306
I^2 (total heterogeneity / total variability):   20.88%
H^2 (total variability / sampling variability):  1.26

Test for Heterogeneity:
Q(df = 8) = 9.698, p-val = 0.287

Model Results:

estimate       se     zval     pval    ci.lb    ci.ub
-0.396    0.227   -1.747    0.081   -0.841    0.048        .

---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

And the estimated average log incidence rate ratio (i.e., $\hat{\mu} = -0.396$) can again be back-transformed through exponentiation:

predict(res, transf=exp, digits=2)
 pred ci.lb ci.ub cr.lb cr.ub
0.67  0.43  1.05  0.32  1.42

The mixed-effects conditional Poisson regression model (i.e., the binomial-normal model) can be fitted with:

res <- rma.glmm(measure="IRR", x1i=x1i, t1i=t1i, x2i=x2i, t2i=t2i, data=dat, model="CM.EL")
print(res, digits=3)
Random-Effects Model (k = 9; tau^2 estimator: ML)
Model Type: Conditional Model with Exact Likelihood

tau^2 (estimated amount of total heterogeneity): 0.123
tau (square root of estimated tau^2 value):      0.350
I^2 (total heterogeneity / total variability):   25.68%
H^2 (total variability / sampling variability):  1.35

Tests for Heterogeneity:
Wld(df = 8) =  9.698, p-val = 0.287
LRT(df = 8) = 11.602, p-val = 0.170

Model Results:

estimate       se     zval     pval    ci.lb    ci.ub
-0.476    0.238   -2.004    0.045   -0.942   -0.010        *

---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
predict(res, transf=exp, digits=2)
 pred ci.lb ci.ub cr.lb cr.ub
0.62  0.39  0.99  0.27  1.42

These results match up with Table VIII (p. 3055) in the article.

### References

Niel-Weise, B. S., Stijnen, T., & van den Broek, P. J. (2007). Anti-infective-treated central venous catheters: A systematic review of randomized controlled trials. Intensive Care Medicine, 33(12), 2058–2068.

Niel-Weise, B. S., Stijnen, T., & van den Broek, P. J. (2008). Anti-infective-treated central venous catheters for total parenteral nutrition or chemotherapy: A systematic review. Journal of Hospital Infection, 69(2), 114–123.

Stijnen, T., Hamza, T. H., & Ozdemir, P. (2010). Random effects meta-analysis of event outcome in the framework of the generalized linear mixed model with applications in sparse data. Statistics in Medicine, 29(29), 3046–3067.

1)
And in fact, these data differ from what is reported in Table VI (p. 3054) in the same article. However, Stijnen et al. (2010) clearly used the data from the appendix for the analysis, so Table VI can be ignored.