The metafor Package

A Meta-Analysis Package for R

User Tools

Site Tools

analyses:berkey1998

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
analyses:berkey1998 [2022/08/03 16:52] – Wolfgang Viechtbaueranalyses:berkey1998 [2023/06/22 11:42] (current) – Wolfgang Viechtbauer
Line 4: Line 4:
  
 Berkey et al. (1998) describe a meta-analytic multivariate model for the analysis of multiple correlated outcomes. The use of the model is illustrated with results from 5 trials comparing surgical and non-surgical treatments for medium-severity periodontal disease. Reported outcomes include the change in probing depth (PD) and attachment level (AL) one year after the treatment. The effect size measure used for this meta-analysis was the (raw) mean difference, calculated in such a way that positive values indicate that surgery was more effective than non-surgical treatment in decreasing the probing depth and increasing the attachment level. The data are provided in Table I in the article and are stored in the dataset ''dat.berkey1998'' that comes with the metafor package: Berkey et al. (1998) describe a meta-analytic multivariate model for the analysis of multiple correlated outcomes. The use of the model is illustrated with results from 5 trials comparing surgical and non-surgical treatments for medium-severity periodontal disease. Reported outcomes include the change in probing depth (PD) and attachment level (AL) one year after the treatment. The effect size measure used for this meta-analysis was the (raw) mean difference, calculated in such a way that positive values indicate that surgery was more effective than non-surgical treatment in decreasing the probing depth and increasing the attachment level. The data are provided in Table I in the article and are stored in the dataset ''dat.berkey1998'' that comes with the metafor package:
 +
 <code rsplus> <code rsplus>
 library(metafor) library(metafor)
Line 9: Line 10:
 dat dat
 </code> </code>
 +
 (I copy the dataset into ''dat'', which is a bit shorter and therefore easier to type further below). The contents of the dataset are: (I copy the dataset into ''dat'', which is a bit shorter and therefore easier to type further below). The contents of the dataset are:
 +
 <code output> <code output>
    trial           author year ni outcome    yi    v1i    v2i    trial           author year ni outcome    yi    v1i    v2i
Line 23: Line 26:
 10     5    Becker et al. 1988 16      AL -0.39 0.0072 0.0304 10     5    Becker et al. 1988 16      AL -0.39 0.0072 0.0304
 </code> </code>
 +
 So, the results from the various trials indicate that surgery is preferable for reducing the probing depth, while non-surgical treatment is preferable for increasing the attachment level. So, the results from the various trials indicate that surgery is preferable for reducing the probing depth, while non-surgical treatment is preferable for increasing the attachment level.
  
Line 30: Line 34:
  
 Before we can proceed with the model fitting, we need to construct the full (block-diagonal) variance-covariance for all studies from these two variables. We can do this using the ''vcalc()'' function in one line of code: Before we can proceed with the model fitting, we need to construct the full (block-diagonal) variance-covariance for all studies from these two variables. We can do this using the ''vcalc()'' function in one line of code:
 +
 <code rsplus> <code rsplus>
 V <- vcalc(vi=1, cluster=author, rvars=c(v1i, v2i), data=dat) V <- vcalc(vi=1, cluster=author, rvars=c(v1i, v2i), data=dat)
 </code> </code>
 +
 The ''V'' matrix is then equal to: The ''V'' matrix is then equal to:
 +
 <code output> <code output>
         [,1]   [,2]   [,3]   [,4]   [,5]   [,6]   [,7]   [,8]   [,9]  [,10]         [,1]   [,2]   [,3]   [,4]   [,5]   [,6]   [,7]   [,8]   [,9]  [,10]
Line 50: Line 57:
 ==== Multivariate Random-Effects Model ==== ==== Multivariate Random-Effects Model ====
  
-A multivariate random-effects model can now be used to meta-analyze the two outcomes simultaneously.+A multivariate random-effects model can now be used to meta-analyze the two outcomes simultaneously. The model is given by $$\left[ \begin{array}{c} 
 +   y_{PD,i}\\ 
 +   y_{AL,i} 
 +\end{array} \right] 
 +   = 
 +\left[ \begin{array}{c} 
 +   \mu_{PD} \\ 
 +   \mu_{AL} 
 +\end{array} \right] 
 +   +  
 +\left[ \begin{array}{c} 
 +   u_{PD,i} \\ 
 +   u_{AL,i} 
 +\end{array} \right] 
 +   +  
 +\left[ \begin{array}{c} 
 +   \epsilon_{PD,i} \\ 
 +   \epsilon_{AL,i} 
 +\end{array} \right],$$ where $$G = \mbox{Var} 
 +\left[ \begin{array}{c} 
 +   u_{PD,i} \\ 
 +   u_{AL,i} 
 +\end{array} \right] 
 +   = 
 +\left[ \begin{array}{cc} 
 +   \tau^2_{PD}              & \rho \tau_{PD} \tau_{AL} \\ 
 +   \rho \tau_{PD} \tau_{AL} & \tau^2_{AL} 
 +\end{array} \right]$$ and $$V_i = \mbox{Var} 
 +\left[ \begin{array}{c} 
 +   \epsilon_{PD,i} \\ 
 +   \epsilon_{AL,i} 
 +\end{array} \right] 
 +   = 
 +\left[ \begin{array}{cc} 
 +   \mbox{var}_{PD,i}      & \mbox{cov}_{PD,i;AL,i} \\ 
 +   \mbox{cov}_{PD,i;AL,i} & \mbox{var}_{AL,i} 
 +\end{array} \right],$$ where the elements in this second variance-covariance matrix are given by the $2 \times 2$ blocks along the diagonal in the ''V'' matrix. Therefore, $$V = \left[ \begin{array}{ccc} 
 +   V_1 &        & \\ 
 +       & \ddots & \\ 
 +       &        & V_5 
 +\end{array} \right].$$  
 + 
 +We can fit this model with: 
 <code rsplus> <code rsplus>
 res <- rma.mv(yi, V, mods = ~ outcome - 1, res <- rma.mv(yi, V, mods = ~ outcome - 1,
Line 88: Line 138:
 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>
 +
 This is what Berkey et al. (1998) call a multivariate maximum likelihood (MML) random-effects model.((Berkey et al. (1998) claim that this approach can only be used when all the trials in the meta-analysis provide results for all the outcomes considered. This does not apply to the way ''rma.mv()'' estimates this model. Therefore, if a particular study had only measured one of the two outcomes, it could still be included in the analysis. The variance-covariance matrix of that study would then simply be a $1 \times 1$ matrix, with the value of that matrix equal to the sampling variance of the observed outcome.)) Note that ''rma.mv()'' uses REML estimation by default, so ''method="ML"'' must be explicitly requested. The ''random = ~ outcome | trial'' part adds random effects for each outcome within each trial to the model. With ''struct="UN"'', the random effects are allowed to have different variances for each outcome and are allowed to be correlated. This is what Berkey et al. (1998) call a multivariate maximum likelihood (MML) random-effects model.((Berkey et al. (1998) claim that this approach can only be used when all the trials in the meta-analysis provide results for all the outcomes considered. This does not apply to the way ''rma.mv()'' estimates this model. Therefore, if a particular study had only measured one of the two outcomes, it could still be included in the analysis. The variance-covariance matrix of that study would then simply be a $1 \times 1$ matrix, with the value of that matrix equal to the sampling variance of the observed outcome.)) Note that ''rma.mv()'' uses REML estimation by default, so ''method="ML"'' must be explicitly requested. The ''random = ~ outcome | trial'' part adds random effects for each outcome within each trial to the model. With ''struct="UN"'', the random effects are allowed to have different variances for each outcome and are allowed to be correlated.
  
Line 95: Line 146:
  
 The results given in Table II in the paper actually are based on a meta-regression model, using year of publication as a potential moderator. To replicate those analyses, we use: The results given in Table II in the paper actually are based on a meta-regression model, using year of publication as a potential moderator. To replicate those analyses, we use:
 +
 <code rsplus> <code rsplus>
 res <- rma.mv(yi, V, mods = ~ outcome + outcome:I(year - 1983) - 1, res <- rma.mv(yi, V, mods = ~ outcome + outcome:I(year - 1983) - 1,
Line 134: Line 186:
 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 publication year was centered at 1983, as was done by the authors. These results correspond to those given in the rightmost column in Table II on page 2545 (column "multiple outcomes MML"). The output above directly provides the correlation among the true effects. We can obtain the covariance with: Note that publication year was centered at 1983, as was done by the authors. These results correspond to those given in the rightmost column in Table II on page 2545 (column "multiple outcomes MML"). The output above directly provides the correlation among the true effects. We can obtain the covariance with:
 +
 <code rsplus> <code rsplus>
 round(res$G[1,2], digits=3) round(res$G[1,2], digits=3)
Line 143: Line 197:
  
 To test whether the slope of publication year actually differs for the two outcomes, we can fit the same model with: To test whether the slope of publication year actually differs for the two outcomes, we can fit the same model with:
 +
 <code rsplus> <code rsplus>
 res <- rma.mv(yi, V, mods = ~ outcome*I(year - 1983) - 1, res <- rma.mv(yi, V, mods = ~ outcome*I(year - 1983) - 1,
Line 149: Line 204:
 print(res, digits=3) print(res, digits=3)
 </code> </code>
 +
 The output is identical, except for the last part, which is now equal to: The output is identical, except for the last part, which is now equal to:
 +
 <code output> <code output>
 Model Results: Model Results:
Line 162: Line 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>
 +
 Therefore, the slope is actually not significantly different for the two outcomes ($p = .553$). In fact, it does not appear as if publication year is at all related to the two outcomes. Therefore, the slope is actually not significantly different for the two outcomes ($p = .553$). In fact, it does not appear as if publication year is at all related to the two outcomes.
  
Line 167: Line 225:
  
 One could actually consider a simpler model for these data, which assumes a compound symmetry structure for the random effects (this would imply that the amount of heterogeneity is the same for the two outcomes). A formal comparison of the two models can be conducted using a likelihood ratio test: One could actually consider a simpler model for these data, which assumes a compound symmetry structure for the random effects (this would imply that the amount of heterogeneity is the same for the two outcomes). A formal comparison of the two models can be conducted using a likelihood ratio test:
 +
 <code rsplus> <code rsplus>
 res1 <- rma.mv(yi, V, mods = ~ outcome - 1, res1 <- rma.mv(yi, V, mods = ~ outcome - 1,
Line 181: Line 240:
 Reduced  4 -2.4111 -1.2008  5.5889 5.2056 1.2702 0.2597 128.2267 Reduced  4 -2.4111 -1.2008  5.5889 5.2056 1.2702 0.2597 128.2267
 </code> </code>
 +
 Since model ''res1'' has one more parameter than model ''res0'', the test statistic (''LRT'') follows (approximately) a chi-square distribution with 1 degree of freedom. The p-value (0.2597) suggests that the more complex model (with ''struct="UN"'') does not actually yield a significantly better fit than the simpler model (with ''struct="CS"''). However, with only 5 studies, this test is likely to have very low power. Since model ''res1'' has one more parameter than model ''res0'', the test statistic (''LRT'') follows (approximately) a chi-square distribution with 1 degree of freedom. The p-value (0.2597) suggests that the more complex model (with ''struct="UN"'') does not actually yield a significantly better fit than the simpler model (with ''struct="CS"''). However, with only 5 studies, this test is likely to have very low power.
  
analyses/berkey1998.txt · Last modified: 2023/06/22 11:42 by Wolfgang Viechtbauer