Multivariate Meta-Analysis

Multivariate Meta-Analysis of Correlated Outcomes

Multivariate meta-analysis gives each study a vector of effect estimates and a covariance matrix, in place of the single estimate and variance used in an ordinary meta-analysis, and fits one model to all the outcomes at once. A trial that reports only progression-free survival can then still inform the pooled effect on overall survival, because the two effects move together. The method matters in health technology assessment because decision models need several clinical inputs at once, and because the NICE manual (PMG36) points to the bivariate methods in NICE Decision Support Unit Technical Support Document 20 for validating a surrogate endpoint. The gain depends on the strength of the correlation and on how many studies leave an outcome unreported.

Within-study and between-study correlation

Two kinds of correlation enter the model. Jackson, Riley and White explain that within-study correlation arises because the effects are calculated from the same patients, so effects on overall and disease-free survival in one trial are almost necessarily positively correlated. Between-study correlation describes how the true effects vary together across studies.

In diagnostic accuracy studies the within-study correlation between sensitivity and specificity is zero, as they come from different patients, but a less stringent positivity threshold raises one and lowers the other, giving a negative correlation across studies. Even so, each outcome keeps its own summary result, as Riley and colleagues note in the BMJ in 2017.

The multivariate random-effects model

The model has a within-study level and a between-study level, as in the univariate random effects model. In the notation of Jackson, Riley and White, the two levels and the resulting marginal model are as follows.

$$Y_i \mid \mu_i \sim N(\mu_i, S_i), \qquad \mu_i \sim N(\mu, \Sigma), \qquad Y_i \sim N(\mu, S_i + \Sigma)$$

Here $Y_i$ is the vector of effect estimates from study $i$, such as log hazard ratios for two outcomes, $\mu_i$ is that study's vector of true effects, $S_i$ is the within-study variance-covariance matrix, treated as known, $\mu$ is the vector of average effects and $\Sigma$ is the between-study covariance matrix. In the bivariate case $\Sigma$ holds two between-study variances and one between-study correlation, and if every entry is zero the model reduces to a common-effect model.

Given an estimate of $\Sigma$, the pooled effects are a generalised least squares average. The weights are inverse covariance matrices, which generalises inverse-variance weighting.

$$\hat\mu = \left( \sum_{i=1}^{n} (S_i + \hat\Sigma)^{-1} \right)^{-1} \sum_{i=1}^{n} (S_i + \hat\Sigma)^{-1} Y_i, \qquad \operatorname{Var}(\hat\mu) = \left( \sum_{i=1}^{n} (S_i + \hat\Sigma)^{-1} \right)^{-1}$$

Here $n$ is the number of studies and $\hat\Sigma$ is the estimated between-study covariance matrix. A study that reports only some outcomes enters through the part of the model for those outcomes, on the assumption that the unreported results are missing at random.

Estimating $\Sigma$ is the main statistical difficulty. Jackson, Riley and White describe restricted maximum likelihood (REML), the more usual choice, which helps to correct the downward bias of maximum likelihood estimation for variance components, and a multivariate method of moments that is the least computationally intensive and needs no between-study normality assumption. TSD 20 fits its models by Bayesian Markov chain Monte Carlo, and Jackson and colleagues advise sensitivity analysis to the prior distribution.

Handling unknown within-study correlations

The model treats $S_i$ as known, but published trials give standard errors for each outcome and almost never the correlation between the two estimates. Jackson, Riley and White call this perhaps the greatest practical difficulty, and TSD 20 adds that ignoring the within-study correlation may increase uncertainty around the pooled effects and potentially bias them. The remedies need different amounts of data.

• Individual participant data allow the correlation to be estimated, for example by bootstrapping, and imputed for other studies or used as an informative prior.
• A sensitivity analysis can impute values across the whole range between minus 1 and 1.
• An alternative model by Riley and colleagues merges the two covariance terms into one parameter, which TSD 20 notes cannot then be read as a within-study or a between-study correlation, limiting its use in surrogate evaluation.

Measuring the strength borrowed from correlated outcomes

Riley and colleagues summarise the extra information with two linked statistics. Both compare the variance of one pooled result with and without the correlated evidence.

$$E = \frac{\operatorname{Var}(\hat\mu_{\text{direct and related}})}{\operatorname{Var}(\hat\mu_{\text{direct only}})}, \qquad \text{BoS} = 100 \times (1 - E)$$

Here $E$ is the efficiency, the numerator uses direct and correlated evidence, the denominator comes from a univariate analysis of the direct evidence alone, and BoS (borrowing of strength) is the percentage reduction in variance due to the correlated evidence. The gain is equivalent to finding $n(1-E)/E$ extra studies of similar size, where $n$ is the number of studies with direct evidence.

In an endometrial cancer review they report, a univariate analysis of 6 studies gave a hazard ratio for cancer-specific survival of 0.61 (95% confidence interval 0.38 to 1.00), while the multivariate analysis of 17 studies, using a correlation of about 0.8 with progression-free survival, gave 0.48 (0.29 to 0.79) with BoS of 33%. TSD 20 adds that with complete data BoS will usually be small, for example under 10%, and cites anecdotal evidence that BoS is roughly bounded by the percentage of missing results.

Worked example with progression-free and overall survival

All numbers in this example are illustrative. Trials A and B of a new cancer drug report log hazard ratios for progression-free survival (PFS) and overall survival (OS), and trial C reports only PFS. Every PFS standard error is 0.10, every OS standard error is 0.15, the within-study correlation is assumed to be 0.6, and a common-effect bivariate model ($\Sigma$ set to zero) keeps the arithmetic short.

1. Combine the complete trials. A and B each have variances 0.01 and 0.0225 and covariance 0.6 × 0.10 × 0.15 = 0.009. Their average, minus 0.35 for PFS and minus 0.15 for OS, has half those values (0.005, 0.01125 and 0.0045).

2. Pool PFS. The A and B average gets weight 1/0.005 = 200 and trial C gets 1/0.01 = 100.

$$\hat\mu_{PFS} = \frac{200 \times (-0.35) + 100 \times (-0.50)}{200 + 100} = -0.40$$

Here $\hat\mu_{PFS}$ is the pooled PFS log hazard ratio, with variance 1/300. It matches the univariate result here because the OS data, all from A and B, carry no extra information about PFS.

3. Pool OS by borrowing from PFS. With $\Sigma$ at zero, the generalised least squares OS estimate equals the A and B average adjusted by the regression of OS on PFS within that average, with slope 0.0045/0.005 = 0.9.

$$\hat\mu_{OS} = \bar{Y}_{OS} + \beta(\hat\mu_{PFS} - \bar{Y}_{PFS}) = -0.15 + 0.9 \times (-0.40 + 0.35) = -0.195$$

Here $\bar{Y}_{OS}$ and $\bar{Y}_{PFS}$ are the A and B averages and $\beta$ is the slope. The univariate OS estimate, from A and B alone, is minus 0.15.

4. Variance and strength borrowed. The variance is the conditional variance of the OS average given the PFS average plus the slope squared times the PFS variance.

$$\operatorname{Var}(\hat\mu_{OS}) = 0.01125 - \frac{0.0045^2}{0.005} + 0.9^2 \times \frac{1}{300} = 0.0099$$

Here 0.01125 is the univariate variance, so the efficiency is 0.0099/0.01125 = 0.88 and BoS is 12%, equivalent to 2 × 0.12/0.88 = 0.27 of an extra trial with OS data. The pooled PFS and OS effects have covariance 0.9 × 1/300 = 0.003, a correlation of about 0.52.

5. Vary the assumed correlation. No trial published its correlation, so the analysis is repeated.

The OS hazard ratio moves to 0.82 from 0.86 because trial C's stronger PFS effect is carried across through the correlation, which is credible only if C's missing OS result behaves like those of A and B. In a partitioned survival model, the probabilistic analysis should draw the two log hazard ratios jointly with correlation 0.52, for example by Cholesky decomposition. A random-effects version would add $\Sigma$ and also borrow through the between-study correlation.

Uses in health technology assessment

The most visible role in appraisal is surrogate validation. PMG36 (section 4.6.7) wants evidence that the relative effect on a surrogate predicts the relative effect on the final outcome, preferably from a meta-analysis of randomised trials reporting both, using the bivariate methods in TSD 20. Section 4.6.11 asks that uncertainty in that relationship be quantified and carried through probabilistic sensitivity analysis.

TSD 20 gives two reasons why correlation between clinical inputs matters in an economic model. It changes the point estimate of the ICER when the model is a non-linear function of the endpoints, and it changes the estimated decision uncertainty, including the probability of cost-effectiveness and value of information. PMG36 (section 4.7.19) asks for evidence on correlation between parameters to be reflected in the probabilistic analysis.

TSD 20 describes a NICE appraisal in ankylosing spondylitis that synthesised disease activity and function scores jointly because both drove costs and QALYs. In a docetaxel analysis by Tan and colleagues that TSD 20 reports, bivariate prediction of an unreported PFS effect allowed a three-state Markov model in place of a two-state one, and the ICER changed to GBP 22,148 per QALY gained from GBP 30,026.

For diagnostic tests, Reitsma and colleagues proposed a bivariate random-effects model that analyses sensitivity and specificity jointly and lets covariates act on each separately. PMG36 (sections 3.4.23 and 4.7.18) notes that the two are usually correlated, and the summary ROC curve page covers that application.

Relation to network meta-analysis

Riley and colleagues contrast the two methods. A multivariate model borrows across outcomes through correlation, while a network meta-analysis borrows across treatments through the consistency assumption. In their experience more is gained across treatments than across outcomes, so a multivariate analysis is best kept for highly correlated outcomes chosen in advance. TSD 20 shows how the two combine in a multivariate network meta-analysis.

Assumptions and common misreadings

The borrowed strength rests on assumptions that are hard to check with few studies. Most of them concern missing results.

1. Assuming results are missing at random. Selective outcome reporting makes missingness not at random. Kirkham and colleagues found by simulation that borrowing can reduce the impact of outcome reporting bias, though Riley and colleagues note it cannot help with trials reporting none of the outcomes.
2. Trusting the normal and linear between-study model. Multivariate normality implies a linear relationship between true effects, and Jackson, Riley and White warn that one outlying result can then sway the estimates for every outcome.
3. Expecting large gains from complete data. Riley and colleagues' 2007 evaluation found similar standard errors to univariate analysis when all outcomes were reported. A formal comparison of effects on two outcomes still needs the joint model for correct confidence intervals.
4. Ignoring the pattern of missing data. In a p53 example from Jackson, Riley and White, the univariate disease-free survival hazard ratio of 0.45 (95% confidence interval 0.27 to 0.74) came from studies whose results differed from those reporting only overall survival, and the bivariate REML result was not statistically significant.
5. Adding many outcomes. Riley and colleagues report that applications beyond two or three outcomes are rare.

Sources

• Bujkiewicz S, Achana F, Papanikos T, Riley RD, Abrams KR. NICE DSU Technical Support Document 20: Multivariate meta-analysis of summary data for combining treatment effects on correlated outcomes and evaluating surrogate endpoints. Version 2. Decision Support Unit; 2019 (updated December 2022). https://www.sheffield.ac.uk/nice-dsu/tsds/multivariate-meta-analysis
• Jackson D, Riley R, White IR. Multivariate meta-analysis: potential and promise. Statistics in Medicine. 2011;30(20):2481-2498. https://doi.org/10.1002/sim.4172
• Kirkham JJ, Riley RD, Williamson PR. A multivariate meta-analysis approach for reducing the impact of outcome reporting bias in systematic reviews. Statistics in Medicine. 2012;31(20):2179-2195. https://doi.org/10.1002/sim.5356
• National Institute for Health and Care Excellence. NICE technology appraisal and highly specialised technologies guidance: the manual (PMG36). Published 31 January 2022, last updated 13 October 2026. Sections 3.4.23, 4.6.7, 4.6.11, 4.7.18 and 4.7.19. https://www.nice.org.uk/process/pmg36
• Reitsma JB, Glas AS, Rutjes AWS, Scholten RJPM, Bossuyt PM, Zwinderman AH. Bivariate analysis of sensitivity and specificity produces informative summary measures in diagnostic reviews. Journal of Clinical Epidemiology. 2005;58(10):982-990. https://doi.org/10.1016/j.jclinepi.2005.02.022
• Riley RD, Abrams KR, Lambert PC, Sutton AJ, Thompson JR. An evaluation of bivariate random-effects meta-analysis for the joint synthesis of two correlated outcomes. Statistics in Medicine. 2007;26(1):78-97. https://doi.org/10.1002/sim.2524
• Riley RD, Jackson D, Salanti G, Burke DL, Price M, Kirkham J, White IR. Multivariate and network meta-analysis of multiple outcomes and multiple treatments: rationale, concepts, and examples. BMJ. 2017;358:j3932. https://doi.org/10.1136/bmj.j3932