# The metafor Package

A Meta-Analysis Package for R

### Sidebar

tips:comp_two_independent_estimates

## Comparing Estimates of Independent Meta-Analyses or Subgroups

Suppose we have summary estimates (e.g., estimated average effects) obtained from two independent meta-analyses or subgroups of studies and we want to test whether the estimates are different from each other. A Wald-type test can be used for this purpose. Alternatively, one could run a single meta-regression model including all studies and using a dichotomous moderator to distinguish the two sets. Both approaches are conceptually very similar with a subtle difference that will be illustrated below with an example.

### Data Preparation

We will use the 'famous' BCG vaccine meta-analysis for this illustration. First, we compute the log risk ratios (and corresponding sampling variances) for each study and then dichotomize the alloc variable.

library(metafor)
dat <- escalc(measure="RR", ai=tpos, bi=tneg, ci=cpos, di=cneg, data=dat.bcg)
dat$alloc <- ifelse(dat$alloc == "random", "random", "other")
dat
   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  other -0.2175 0.0512
6      6      Stein & Aronson 1953  180  1361  372  1079    44  other -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  other -1.3713 0.0730
11    11       Comstock et al 1974  186 50448  141 27197    18  other -0.3394 0.0124
12    12   Comstock & Webster 1969    5  2493    3  2338    33  other  0.4459 0.5325
13    13       Comstock et al 1976   27 16886   29 17825    33  other -0.0173 0.0714

### Separate Meta-Analyses

First, we fit two separate random-effects models within each subset defined by the alloc variable:

res1 <- rma(yi, vi, data=dat, subset=alloc=="random")
res2 <- rma(yi, vi, data=dat, subset=alloc=="other")

We then combine the estimates and standard errors from each model into a data frame. We also add a variable to distinguish the two models and, for reasons to be explained in more detail below, we add the estimated amounts of heterogeneity within each subset to the data frame.

dat.comp <- data.frame(estimate = c(coef(res1), coef(res2)), stderror = c(res1$se, res2$se),
meta = c("random","other"), tau2 = round(c(res1$tau2, res2$tau2),3))
dat.comp
    estimate  stderror   meta  tau2
1 -0.9709645 0.2759557 random 0.393
2 -0.4812706 0.2169886  other 0.212

We can now compare the two estimates (i.e., the estimated average log risk ratios) by feeding them back to the rma() function and using the variable to distinguish the two estimates as a moderator. We use a fixed-effects model, because the (residual) heterogeneity within each subset has already been accounted for by fitting random-effects model above.

rma(estimate, sei=stderror, mods = ~ meta, method="FE", data=dat.comp, digits=3)
Fixed-Effects with Moderators Model (k = 2)

Test for Residual Heterogeneity:
QE(df = 0) = 0.000, p-val = 1.000

Test of Moderators (coefficient(s) 2):
QM(df = 1) = 1.946, p-val = 0.163

Model Results:

estimate     se    zval   pval   ci.lb   ci.ub
intrcpt       -0.481  0.217  -2.218  0.027  -0.907  -0.056  *
metarandom    -0.490  0.351  -1.395  0.163  -1.178   0.198

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

While we find that studies using random assignment obtain larger (more negative) effects than studies not using random assignment ($b_1 = -0.490$, $SE = 0.351$), the difference between the two estimates is not significant ($z = -1.395$, $p = .163$).

The test of the difference between the two estimates is really just a Wald-type test, given by the equation $$z = \frac{\hat{\mu}_1 - \hat{\mu}_2}{\sqrt{SE[\hat{\mu}_1]^2 + SE[\hat{\mu}_2]^2}},$$ where $\hat{\mu}_1$ and $\hat{\mu}_2$ are the two estimates and $SE[\hat{\mu}_1]$ and $SE[\hat{\mu}_2]$ the corresponding standard errors. The test statistics can therefore also be computed with:

with(dat.comp, round(c(zval = (estimate[1] - estimate[2])/sqrt(stderror[1]^2 + stderror[2]^2)), 3))
   zval
-1.395

This is the same value as obtained above.

### Meta-Regression with All Studies

Now let's take a different approach, fitting a meta-regression model using all studies:

rma(yi, vi, mods = ~ alloc, data=dat, digits=3)
Mixed-Effects Model (k = 13; tau^2 estimator: REML)

tau^2 (estimated amount of residual heterogeneity):     0.318 (SE = 0.178)
tau (square root of estimated tau^2 value):             0.564
I^2 (residual heterogeneity / unaccounted variability): 89.92%
H^2 (unaccounted variability / sampling variability):   9.92
R^2 (amount of heterogeneity accounted for):            0.00%

Test for Residual Heterogeneity:
QE(df = 11) = 138.511, p-val < .001

Test of Moderators (coefficient(s) 2):
QM(df = 1) = 1.833, p-val = 0.176

Model Results:

estimate     se    zval   pval   ci.lb  ci.ub
intrcpt        -0.467  0.257  -1.816  0.069  -0.972  0.037  .
allocrandom    -0.490  0.362  -1.354  0.176  -1.199  0.219

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

The result is very similar to what we saw earlier: The coefficient for the alloc dummy is equal to $b_1 = -0.490$ ($SE = 0.362$) and not significant ($p = .176$).

However, the results are not exactly identical. The reason for this is as follows. When we fit separate random-effects models in the two subsets, we are allowing the amount of heterogeneity within each set to be different (as shown earlier, the estimates were $\hat{\tau}^2 = 0.393$ and $\hat{\tau}^2 = 0.212$ for studies using and not using random assignment, respectively). On the other hand, the mixed-effects meta-regression model fitted above has a single variance component for the amount of residual heterogeneity, which implies that the amount of heterogeneity within each subset is assumed to be the same ($\hat{\tau}^2 = 0.318$ in this example).

### Meta-Regression with All Studies but Different Amounts of (Residual) Heterogeneity

Using the rma.mv() function, we can easily fit a meta-regression model using all studies where we allow the amount of residual heterogeneity to be different in each subset:

rma.mv(yi, vi, mods = ~ alloc, random = ~ alloc | trial, struct="DIAG", data=dat, digits=3)
Multivariate Meta-Analysis Model (k = 13; method: REML)

Variance Components:

outer factor: trial (nlvls = 13)
inner factor: alloc (nlvls = 2)

estim   sqrt  k.lvl  fixed   level
tau^2.1    0.212  0.460      6     no   other
tau^2.2    0.393  0.627      7     no  random

Test for Residual Heterogeneity:
QE(df = 11) = 138.511, p-val < .001

Test of Moderators (coefficient(s) 2):
QM(df = 1) = 1.946, p-val = 0.163

Model Results:

estimate     se    zval   pval   ci.lb   ci.ub
intrcpt        -0.481  0.217  -2.218  0.027  -0.907  -0.056  *
allocrandom    -0.490  0.351  -1.395  0.163  -1.178   0.198

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

Note that the two estimates of $\tau^2$ are now identical to the ones we obtained earlier from the separate random-effects models. Also, the coefficient, standard error, and p-value for the moderator now matches the results obtained earlier.