tips:comp_two_independent_estimates

Suppose we have summary estimates (e.g., estimated average effects) obtained from two independent meta-analyses or two subgroups of studies within the same meta-analysis 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.

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

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 models 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 on average larger (i.e., 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 that we obtained above.

Now let's take a different approach, fitting a meta-regression model with `alloc`

as a categorical moderator based on 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 corresponding to 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).

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.

A discussion/comparison of these two approaches (i.e., assuming a single $\tau^2$ value or allowing $\tau^2$ to differ across subsets) can be found in the following article:

Rubio-Aparicio, M., López-López, J. A., Viechtbauer, W., Marín-Martínez, F., Botella, J., & Sánchez-Meca, J. (2020). Testing categorical moderators in mixed-effects meta-analysis in the presence of heteroscedasticity. *Journal of Experimental Education, 88*(2), 288-310. https://doi.org/10.1080/00220973.2018.1561404

We can also do a likelihood ratio test (LRT) to examine whether there are significant differences in the $\tau^2$ values across subsets. This can be done with:

res1 <- rma.mv(yi, vi, mods = ~ alloc, random = ~ alloc | trial, struct="DIAG", data=dat) res0 <- rma.mv(yi, vi, mods = ~ alloc, random = ~ alloc | trial, struct="ID", data=dat) anova(res1, res0)

df AIC BIC AICc logLik LRT pval QE Full 4 29.2959 30.8875 35.9626 -10.6480 138.5113 Reduced 3 27.5948 28.7885 31.0234 -10.7974 0.2989 0.5845 138.5113

So in this example, we would not reject the null hypothesis $H_0: \tau^2_1 = \tau^2_2$ ($p = .58$).

tips/comp_two_independent_estimates.txt · Last modified: 2021/05/04 11:06 by Wolfgang Viechtbauer

Except where otherwise noted, content on this wiki is licensed under the following license: CC Attribution-Noncommercial-Share Alike 4.0 International