# The metafor Package

A Meta-Analysis Package for R

### Site Tools

analyses:berkey1998

## Berkey et al. (1998)

### The Methods and Data

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:

library(metafor)
dat <- dat.berkey1998
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:

   trial           author year ni outcome    yi    v1i    v2i
1      1 Pihlstrom et al. 1983 14      PD  0.47 0.0075 0.0030
2      1 Pihlstrom et al. 1983 14      AL -0.32 0.0030 0.0077
3      2    Lindhe et al. 1982 15      PD  0.20 0.0057 0.0009
4      2    Lindhe et al. 1982 15      AL -0.60 0.0009 0.0008
5      3   Knowles et al. 1979 78      PD  0.40 0.0021 0.0007
6      3   Knowles et al. 1979 78      AL -0.12 0.0007 0.0014
7      4  Ramfjord et al. 1987 89      PD  0.26 0.0029 0.0009
8      4  Ramfjord et al. 1987 89      AL -0.31 0.0009 0.0015
9      5    Becker et al. 1988 16      PD  0.56 0.0148 0.0072
10     5    Becker et al. 1988 16      AL -0.39 0.0072 0.0304

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.

Since each trial provides effect size estimates for both outcomes, the estimates are correlated. The v1i and v2i values are the variances and covariances of the observed effects. In particular, for each study, variables v1i and v2i form a 2×2 variance-covariance matrix of the observed effects, with the diagonal elements corresponding to the sampling variances of the mean differences (the first for probing depth, the second for attachment level) and the off-diagonal value corresponding to the covariance of the two mean differences.

### Variance-Covariance Matrix

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:

V <- vcalc(vi=1, cluster=author, rvars=c(v1i, v2i), data=dat)

The V matrix is then equal to:

        [,1]   [,2]   [,3]   [,4]   [,5]   [,6]   [,7]   [,8]   [,9]  [,10]
[1,] 0.0075 0.0030 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000
[2,] 0.0030 0.0077 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000
[3,] 0.0000 0.0000 0.0057 0.0009 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000
[4,] 0.0000 0.0000 0.0009 0.0008 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000
[5,] 0.0000 0.0000 0.0000 0.0000 0.0021 0.0007 0.0000 0.0000 0.0000 0.0000
[6,] 0.0000 0.0000 0.0000 0.0000 0.0007 0.0014 0.0000 0.0000 0.0000 0.0000
[7,] 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0029 0.0009 0.0000 0.0000
[8,] 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0009 0.0015 0.0000 0.0000
[9,] 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0148 0.0072
[10,] 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0072 0.0304

### Multivariate Random-Effects Model

A multivariate random-effects model can now be used to meta-analyze the two outcomes simultaneously.

res <- rma.mv(yi, V, mods = ~ outcome - 1,
random = ~ outcome | trial, struct="UN",
data=dat, method="ML")
print(res, digits=3)
Multivariate Meta-Analysis Model (k = 10; method: ML)

Variance Components:

outer factor: trial   (nlvls = 5)
inner factor: outcome (nlvls = 2)

estim   sqrt  k.lvl  fixed  level
tau^2.1    0.026  0.162      5     no     AL
tau^2.2    0.007  0.084      5     no     PD

rho.AL  rho.PD    AL  PD
AL       1   0.699     -  no
PD   0.699       1     5   -

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

Test of Moderators (coefficient(s) 1,2):
QM(df = 2) = 155.772, p-val < .001

Model Results:

estimate     se    zval   pval   ci.lb   ci.ub
outcomeAL    -0.338  0.080  -4.237  <.001  -0.494  -0.182  ***
outcomePD     0.345  0.049   6.972  <.001   0.248   0.442  ***

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

This is what Berkey et al. (1998) call a multivariate maximum likelihood (MML) random-effects model.1) 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.

The results show that the amount of heterogeneity in the attachment level (AL) outcome (i.e., $\hat{\tau}_{AL}^2 = .026$) is larger than the amount of heterogeneity in the probing depth (PD) outcome (i.e., $\hat{\tau}_{PD}^2 = .007$). Furthermore, the true outcomes appear to correlate quite strongly (i.e., $\hat{\rho} = .70$). On average, surgery is estimated to lead to significantly greater decreases in probing depth (i.e., $\hat{\mu}_{PB} = .35$), but non-surgery is more effective in increasing the attachment level (i.e., $\hat{\mu}_{AL} = -.34$).

### Multivariate Mixed-Effects Meta-Regression Model

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:

res <- rma.mv(yi, V, mods = ~ outcome + outcome:I(year - 1983) - 1,
random = ~ outcome | trial, struct="UN",
data=dat, method="ML")
print(res, digits=3)
Multivariate Meta-Analysis Model (k = 10; method: ML)

Variance Components:

outer factor: trial   (nlvls = 5)
inner factor: outcome (nlvls = 2)

estim   sqrt  k.lvl  fixed  level
tau^2.1    0.025  0.158      5     no     AL
tau^2.2    0.008  0.090      5     no     PD

rho.AL  rho.PD    AL  PD
AL       1   0.659     -  no
PD   0.659       1     5   -

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

Test of Moderators (coefficient(s) 1,2,3,4):
QM(df = 4) = 143.439, p-val < .001

Model Results:

estimate     se    zval   pval   ci.lb   ci.ub
outcomeAL                   -0.335  0.079  -4.261  <.001  -0.489  -0.181  ***
outcomePD                    0.348  0.052   6.694  <.001   0.246   0.450  ***
outcomeAL:I(year - 1983)    -0.011  0.024  -0.445  0.656  -0.059   0.037
outcomePD:I(year - 1983)     0.001  0.015   0.063  0.950  -0.029   0.031

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

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: