The metafor Package

A Meta-Analysis Package for R

User Tools

Site Tools

tips:comp_two_independent_estimates

Differences

This shows you the differences between two versions of the page.

Link to this comparison view

Both sides previous revisionPrevious revision
Next revision		Previous revision
Next revisionBoth sides next revision
tips:comp_two_independent_estimates [2021/11/08 15:49] – Wolfgang Viechtbauertips:comp_two_independent_estimates [2023/02/18 17:35] – Wolfgang Viechtbauer
Line 6: Line 6:
  
 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. 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.
 +
 <code rsplus> <code rsplus>
 library(metafor) library(metafor)
Line 12: Line 13:
 dat dat
 </code> </code>
 +
 <code output> <code output>
    trial               author year tpos  tneg cpos  cneg ablat  alloc      yi     vi    trial               author year tpos  tneg cpos  cneg ablat  alloc      yi     vi
Line 32: Line 34:
  
 First, we fit two separate random-effects models within each subset defined by the ''alloc'' variable: First, we fit two separate random-effects models within each subset defined by the ''alloc'' variable:
 +
 <code rsplus> <code rsplus>
 res1 <- rma(yi, vi, data=dat, subset=alloc=="random") res1 <- rma(yi, vi, data=dat, subset=alloc=="random")
Line 38: Line 41:
  
 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. 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.
 +
 <code rsplus> <code rsplus>
 dat.comp <- data.frame(estimate = c(coef(res1), coef(res2)), stderror = c(res1$se, res2$se), dat.comp <- data.frame(estimate = c(coef(res1), coef(res2)), stderror = c(res1$se, res2$se),
Line 43: Line 47:
 dat.comp dat.comp
 </code> </code>
 +
 <code output> <code output>
     estimate  stderror   meta  tau2     estimate  stderror   meta  tau2
Line 50: Line 55:
  
 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 meta-regression model for this purpose, because the (residual) heterogeneity within each subset has already been accounted for by fitting random-effects models above. 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 meta-regression model for this purpose, because the (residual) heterogeneity within each subset has already been accounted for by fitting random-effects models above.
 +
 <code rsplus> <code rsplus>
 rma(estimate, sei=stderror, mods = ~ meta, method="FE", data=dat.comp, digits=3) rma(estimate, sei=stderror, mods = ~ meta, method="FE", data=dat.comp, digits=3)
 </code> </code>
 +
 <code output> <code output>
 Fixed-Effects with Moderators Model (k = 2) Fixed-Effects with Moderators Model (k = 2)
  
-Test for Residual Heterogeneity: +Test for Residual Heterogeneity:
 QE(df = 0) = 0.000, p-val = 1.000 QE(df = 0) = 0.000, p-val = 1.000
  
-Test of Moderators (coefficient(s) 2): +Test of Moderators (coefficient(s) 2):
 QM(df = 1) = 1.946, p-val = 0.163 QM(df = 1) = 1.946, p-val = 0.163
  
 Model Results: Model Results:
  
-            estimate     se    zval   pval   ci.lb   ci.ub   +            estimate     se    zval   pval   ci.lb   ci.ub
 intrcpt       -0.481  0.217  -2.218  0.027  -0.907  -0.056  * 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   +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 Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
 </code> </code>
 +
 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$). 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 [[https://en.wikipedia.org/wiki/Wald_test|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: The test of the difference between the two estimates is really just a [[https://en.wikipedia.org/wiki/Wald_test|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:
 +
 <code rsplus> <code rsplus>
 with(dat.comp, round(c(zval = (estimate[1] - estimate[2])/sqrt(stderror[1]^2 + stderror[2]^2)), 3)) with(dat.comp, round(c(zval = (estimate[1] - estimate[2])/sqrt(stderror[1]^2 + stderror[2]^2)), 3))
 </code> </code>
 +
 <code output> <code output>
-  zval +  zval
 -1.395 -1.395
 </code> </code>
 +
 This is the same value that we obtained above. This is the same value that we obtained above.
  
Line 86: Line 97:
  
 Now let's take a different approach, fitting a meta-regression model with ''alloc'' as a categorical moderator based on all studies: Now let's take a different approach, fitting a meta-regression model with ''alloc'' as a categorical moderator based on all studies:
 +
 <code rsplus> <code rsplus>
 rma(yi, vi, mods = ~ alloc, data=dat, digits=3) rma(yi, vi, mods = ~ alloc, data=dat, digits=3)
 </code> </code>
 +
 <code output> <code output>
 Mixed-Effects Model (k = 13; tau^2 estimator: REML) Mixed-Effects Model (k = 13; tau^2 estimator: REML)
Line 98: Line 111:
 R^2 (amount of heterogeneity accounted for):            0.00% R^2 (amount of heterogeneity accounted for):            0.00%
  
-Test for Residual Heterogeneity: +Test for Residual Heterogeneity:
 QE(df = 11) = 138.511, p-val < .001 QE(df = 11) = 138.511, p-val < .001
  
-Test of Moderators (coefficient(s) 2): +Test of Moderators (coefficient(s) 2):
 QM(df = 1) = 1.833, p-val = 0.176 QM(df = 1) = 1.833, p-val = 0.176
  
 Model Results: Model Results:
  
-             estimate     se    zval   pval   ci.lb  ci.ub   +             estimate     se    zval   pval   ci.lb  ci.ub
 intrcpt        -0.467  0.257  -1.816  0.069  -0.972  0.037  . 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   +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 Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
 </code> </code>
 +
 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$). 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$).
  
Line 120: Line 134:
  
 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: 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:
 +
 <code rsplus> <code rsplus>
 rma.mv(yi, vi, mods = ~ alloc, random = ~ alloc | trial, struct="DIAG", data=dat, digits=3) rma.mv(yi, vi, mods = ~ alloc, random = ~ alloc | trial, struct="DIAG", data=dat, digits=3)
 </code> </code>
 +
 <code output> <code output>
 Multivariate Meta-Analysis Model (k = 13; method: REML) Multivariate Meta-Analysis Model (k = 13; method: REML)
  
-Variance Components: +Variance Components:
  
 outer factor: trial (nlvls = 13) outer factor: trial (nlvls = 13)
Line 135: Line 151:
 tau^2.2    0.393  0.627      7     no  random tau^2.2    0.393  0.627      7     no  random
  
-Test for Residual Heterogeneity: +Test for Residual Heterogeneity:
 QE(df = 11) = 138.511, p-val < .001 QE(df = 11) = 138.511, p-val < .001
  
-Test of Moderators (coefficient(s) 2): +Test of Moderators (coefficient(s) 2):
 QM(df = 1) = 1.946, p-val = 0.163 QM(df = 1) = 1.946, p-val = 0.163
  
 Model Results: Model Results:
  
-             estimate     se    zval   pval   ci.lb   ci.ub   +             estimate     se    zval   pval   ci.lb   ci.ub
 intrcpt        -0.481  0.217  -2.218  0.027  -0.907  -0.056  * 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   +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 Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
 </code> </code>
-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.+ 
 +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.((See also [[tips:different_tau2_across_subgroups|here]] for another example to illustrate various approaches for allowing $\tau^2$ to differ across subgroups.))
  
 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: 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:
Line 163: Line 180:
 anova(res1, res0) anova(res1, res0)
 </code> </code>
 +
 <code output> <code output>
-        df     AIC     BIC    AICc   logLik    LRT   pval       QE  +        df     AIC     BIC    AICc   logLik    LRT   pval       QE 
-Full     4 29.2959 30.8875 35.9626 -10.6480               138.5113 +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 Reduced  3 27.5948 28.7885 31.0234 -10.7974 0.2989 0.5845 138.5113
 </code> </code>
  
 So in this example, we would not reject the null hypothesis $H_0: \tau^2_1 = \tau^2_2$ ($p = .58$). 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: 2024/04/18 11:36 by Wolfgang Viechtbauer