Meta-analysis for multiple outcomes with random effects model

Question

library(metafor)
dat <- get(data(dat.berkey1998))
V <- bldiag(lapply(split(dat[,c("v1i", "v2i")], dat$trial), as.matrix))
res <- rma.mv(yi, V, mods = ~ outcome - 1, random = ~ outcome | trial, struct="UN", data=dat, method="ML")
print(res, digits=3)

> dat
   trial           author year ni outcome      yi     vi    v1i    v2i
1      1 Pihlstrom et al. 1983 14      PD  0.4700 0.0075 0.0075 0.0030
2      1 Pihlstrom et al. 1983 14      AL -0.3200 0.0077 0.0030 0.0077
3      2    Lindhe et al. 1982 15      PD  0.2000 0.0057 0.0057 0.0009
4      2    Lindhe et al. 1982 15      AL -0.6000 0.0008 0.0009 0.0008
5      3   Knowles et al. 1979 78      PD  0.4000 0.0021 0.0021 0.0007
6      3   Knowles et al. 1979 78      AL -0.1200 0.0014 0.0007 0.0014
7      4  Ramfjord et al. 1987 89      PD  0.2600 0.0029 0.0029 0.0009
8      4  Ramfjord et al. 1987 89      AL -0.3100 0.0015 0.0009 0.0015
9      5    Becker et al. 1988 16      PD  0.5600 0.0148 0.0148 0.0072
10     5    Becker et al. 1988 16      AL -0.3900 0.0304 0.0072 0.0304

I am trying to understand this example documented here. Each study has 2 correlated outcomes (PD and AL). And the goal of this analysis is to estimate a combined effect size fo PD, and to estimate a combined effect size for AL, correct?

In the following forest plot:

What exactly are the grey diamonds overlapping the effect size & CI for each study? And why is there no combined effect at the bottom of the forest plot? In my mind I was expecting 2 forest plots: 1 for the combined effect size of PD, and one of the combined effect size of AL, am I understanding the analysis correctly?

res2

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.0261  0.1617      5     no     AL
tau^2.2    0.0070  0.0837      5     no     PD

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

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

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

Model Results:

           estimate      se     zval    pval    ci.lb    ci.ub     
outcomeAL   -0.3379  0.0798  -4.2368  <.0001  -0.4943  -0.1816  ***
outcomePD    0.3448  0.0495   6.9721  <.0001   0.2479   0.4418  ***

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

A quick look at the summary of res shows me an estimate of AL and PD. Are these the combined effect sizes for AL, and PD, respectively? And what exactly does mods = ~ outcome - 1 do?

Answer 1

Yes, the values given for outcomeAL and outcomePD are the estimated average effects (in this case, mean differences) for the two outcomes.

By using mods = ~ outcome - 1, you get the two estimated average effects directly. If we had used mods = ~ outcome, we would get the estimated average effect for outcome AL and then a coefficient that is the difference (i.e., contrast) between the estimated average effect for outcome AL and outcome PD.

The gray diamonds/polygons in the forest plot show the predicted effects (i.e., fitted values) for each study. That is what is automatically shown when drawing a forest plot based on a model with moderator variables. In this example, the fitted values are just the values given for outcomeAL and outcomePD (so, for "studies" 1, 3, 5, 7, and 9, the value for outcomePD; for "studies" 2, 4, 6, 8, and 19, the value for outcomeAL). In this case, it probably would make things clearer by not drawing the gray polygons (use addfit=FALSE) and instead drawing two summary polygons at the bottom of the forest plot (this can be done using the addpoly() function).