»  Home »  Products »  Stata 17 »  Multivariate meta-analysis

Multivariate meta-analysis


  • Multivariate meta-analysis

    • Fixed effects
    • Random effects
    • Four estimation methods
    • Adjust for moderators
    • Sensitivity analysis
    • Between-study covariance structures
    • Jackson–Riley standard-error adjustments
    • Multivariate \( Q \) statistic and test
  • Heterogeneity

    • Cochran's multivariate statistics
    • Jackson–White–Riley multivariate statistics
    • White's multivariate statistics
  • Postestimation

    • Predict random effects
    • Variance–covariance matrix of random effects
    • Residuals
    • Standardized residuals

Univariate meta-analysis deals with a single effect reported by each study. However, there are many cases in practice where a study may report multiple effect sizes. For example, does the keto diet, the high-protein diet, the vegan diet, or intermittent fasting achieve the highest amount of weight loss? Each of these comparisons will generate an effect size. Modeling each effect separately ignores the fact that they may be correlated. Multivariate meta-analysis models the effects jointly and accounts for their dependence.

Let's see it work

Multivariate meta-analysis

Suppose we are interested in investigating the effect of multiple dietary regimens on weight loss. Multiple effect sizes that compare each of these diets with a control group (not following any diet) can be computed. These effect sizes are usually correlated because they share a common control group. Or perhaps we want to explore the impact of a new teaching technique on math (outcome 1), physics (outcome 2), and chemistry (outcome 3) testing scores. Three effect sizes that compare the three testing scores across two groups of students (those who were taught using the new technique and those who were not) are computed. These three effect sizes are dependent because they were calculated using the same set of students.

Example dataset: Treatment of moderate periodontal disease

Consider a dataset from Antczak-Bouckoms et al. (1993) of five randomized controlled trials that explored the impact of two (surgical and nonsurgical) procedures on treating periodontal disease. Two outcomes of interest are improvements from baseline (pretreatment) in probing depth (y1) and attachment level (y2) around the teeth. The main objectives of the periodontal treatment were to reduce probing depths and increase attachment levels (Berkey et al. 1998). Because the two outcomes y1 and y2 are measured on the same subject, they should not be treated as independent.

. webuse periodontal
(Treatment of moderate periodontal disease)

. describe

Contains data from https://www.stata-press.com/data/r17/periodontal.dta
 Observations:             5                  Treatment of moderate periodontal disease
    Variables:             9                  13 Jan 2021 18:11
                                              (_dta has notes)
Variable Storage Display Value
name type format label Variable label
trial str23 %23s Trial label
pubyear byte %9.0g Publication year centered at 1983
y1 float %6.2f Mean improvement in probing depth (mm)
y2 float %6.2f Mean improvement in attachment level (mm)
v11 float %6.4f Variance of y1
v12 float %6.4f Covariance of y1 and y2
v22 float %6.4f Variance of y2
s1 double %10.0g Standard error of y1
s2 double %10.0g Standard error of y2
Sorted by:

Variables v11, v12, and v22 define the within-study covariance matrix for each study.

Constant-only model: Multivariate meta-analysis

If we were to perform two separate univariate meta-analyses for outcomes y1 and y2, we would be ignoring the dependence among the two outcomes, which may lead to incorrect inference. We use the command meta mvregress to conduct a bivariate meta-analysis as follows:

. meta mvregress y1 y2, wcovvariables(v11 v12 v22)

Performing EM optimization ...

Performing gradient-based optimization:
Iteration 0:   log restricted-likelihood =  2.0594015
Iteration 1:   log restricted-likelihood =  2.0822925
Iteration 2:   log restricted-likelihood =  2.0823276
Iteration 3:   log restricted-likelihood =  2.0823276

Multivariate random-effects meta-analysis       Number of obs     =         10
Method: REML                                    Number of studies =          5
                                                Obs per study:
                                                              min =          2
                                                              avg =        2.0
                                                              max =          2
                                                Wald chi2(0)      =          .
Log restricted-likelihood = 2.0823276           Prob > chi2       =          .

Coefficient Std. err. z P>|z| [95% conf. interval]
_cons .3534282 .0588486 6.01 0.000 .238087 .4687694
_cons -.3392152 .0879051 -3.86 0.000 -.5115061 -.1669243
Test of homogeneity: Q_M = chi2(8) = 128.23 Prob > Q_M = 0.0000
Random-effects parameters Estimate
sd(y1) .1083191
sd(y2) .1806968
corr(y1,y2) .6087987

The order in which you specify variables within wcovvariables() is important (see wcovvariables() in [META] meta mvregress for details).

The first table displays the regression (fixed-effects) coefficient estimates from the bivariate meta-analysis. These estimates correspond to the overall bivariate effect sizes. The overall improvement in probing depth is roughly 0.35 mm, and the overall attachment level was reduced by 0.34 mm.

The multivariate homogeneity test, which tests whether the bivariate effect sizes (\( \theta_{1j}, \theta_{2j} \)) are constant across studies, is rejected (p < 0.0001).

The second table displays the standard deviations of the random effects corresponding to outcomes y1 and y2, as well as their correlation.

We could have performed a fixed-effects multivariate meta-analysis by specifying option fixed:

. meta mvregress y1 y2, wcovvariables(v11 v12 v22) fixed
(output omitted)

By performing a fixed-effects multivariate meta-analysis, we assume that study-specific bivariate effect sizes are the same across the studies and that the observed variability is due to sampling error. This assumption is often not satisfied in practice.

Incorporating moderators: Multivariate meta-regression

Berkey et al. (1998) argued that as the surgical experience accumulates, the surgical procedure will become more efficient, so the most recent studies may show greater surgical benefits. We will include variable pubyear, a surrogate for the time when the trial was performed, as a moderator to explain a portion of the heterogeneity highlighted in the previous section.

. meta mvregress y1 y2 = pubyear, wcovvariables(v*)

Performing EM optimization ...

Performing gradient-based optimization:
Iteration 0:   log restricted-likelihood = -3.5544446
Iteration 1:   log restricted-likelihood = -3.5402141
Iteration 2:   log restricted-likelihood = -3.5399568
Iteration 3:   log restricted-likelihood = -3.5399567

Multivariate random-effects meta-regression     Number of obs     =         10
Method: REML                                    Number of studies =          5
                                                Obs per study:
                                                              min =          2
                                                              avg =        2.0
                                                              max =          2
                                                Wald chi2(2)      =       0.40
Log restricted-likelihood = -3.5399567          Prob > chi2       =     0.8197

Coefficient Std. err. z P>|z| [95% conf. interval]
pubyear .0048615 .0218511 0.22 0.824 -.0379658 .0476888
_cons .3587569 .07345 4.88 0.000 .2147975 .5027163
pubyear -.0115367 .0299635 -0.39 0.700 -.070264 .0471907
_cons -.3357368 .0979979 -3.43 0.001 -.5278091 -.1436645
Test of homogeneity: Q_M = chi2(6) = 125.76 Prob > Q_M = 0.0000
Random-effects parameters Estimate
sd(y1) .1429917
sd(y2) .2021314
corr(y1,y2) .561385

In the wcovvariables() option, we used the stub notation v* to refer to variables v11, v12, and v22.

The estimates of the regression coefficients of variable pubyear are 0.0049 with a 95% CI of [\(−\)0.0380, 0.0477] for outcome y1 and \(-\)0.0115 with a 95% CI of [\(−\)0.0703, 0.0472] for outcome y2. The coefficients are not significant according to the z tests, with the respective p-values, p = 0.824 and p = 0.7. It appears that pubyear did not explain the heterogeneity among the effect sizes.

Random-effects estimation methods and between-study covariance structure

Here we modify the default REML estimation method and use the ML estimation instead. We also use an independent covariance structure for the random effects. This may be done by specifying random(mle, covariance(independent)).

. meta mvregress y*, wcovvariables(v*) random(mle, covariance(independent))

Performing EM optimization ...

Performing gradient-based optimization:
Iteration 0:   log likelihood =  5.1641932
Iteration 1:   log likelihood =  5.1654142
Iteration 2:   log likelihood =  5.1654153
Iteration 3:   log likelihood =  5.1654153

Multivariate random-effects meta-analysis       Number of obs     =         10
Method: ML                                      Number of studies =          5
                                                Obs per study:
                                                              min =          2
                                                              avg =        2.0
                                                              max =          2
                                                Wald chi2(0)      =          .
Log likelihood = 5.1654153                      Prob > chi2       =          .

Coefficient Std. err. z P>|z| [95% conf. interval]
_cons .3572553 .0499616 7.15 0.000 .2593323 .4551782
_cons -.3538886 .0788344 -4.49 0.000 -.5084013 -.199376
Test of homogeneity: Q_M = chi2(8) = 128.23 Prob > Q_M = 0.0000
Random-effects parameters Estimate
sd(y1) .0845313
sd(y2) .1596039

The random-effects parameter table now reports two terms, sd(y1) and sd(y2), because the correlation is assumed to be 0 under the independent covariance structure assumption.

Also see [META] meta mvregress.

estat heterogeneity

After you fit your multivariate meta-analysis model, you should quantify the amount of heterogeneity among the studies that was not accounted for by the model. You can use estat heterogeneity to do that.

. estat heterogeneity

Method: Cochran
  I2 (%) = 93.76
      H2 = 16.03

Method: Jackson–White–Riley
  I2 (%) = 67.29
       R =  1.75

  I2 (%) = 94.40
       R =  4.23

  I2 (%) = 87.49
       R =  2.83

This command produces heterogeneity statistics that extend the concept of univariate heterogeneity statistics, such as \( Q \) and \( I^2 \), to the multivariate setting.

For instance, Cochran's \( I^2 = \) 93.76% means that 93.76% of the heterogeneity is due to true heterogeneity between the studies as opposed to the sampling variability.

One potential shortcoming of the Cochran statistics is that they quantify the amount of heterogeneity jointly for all outcomes. The Jackson–White–Riley statistics provide ways to assess the contribution of each outcome to the total heterogeneity, in addition to their joint contribution.

For example, we can see that there is more heterogeneity among the effect sizes of outcome y2 (\( I^2 = \) 94.40%) than among the effect sizes of y1 (\( I^2 = \) 67.29%)

Also see [META] estat heterogeneity.

Postestimation: Predicting random effects

. predict double u*, reffects reses(se_u*)

. list trial u* se_u*

  trial u1 u2 se_u1 se_u2
1. Philstrom et al. (1983) .05452276 .00844496 .05419716 .12705759
2. Lindhe et al. (1982) -.0829853 -.22848358 .05668313 .13754507
3. Knowles et al. (1979) .02838328 .21906828 .06359432 .13645361
4. Ramfjord et al. (1987) -.07043129 .04982566 .06192739 .13633583
5. Becker et al. (1988) .07051055 -.04885532 .04624526 .10346863

We listed the random-effects variables u1 and u2 with their corresponding standard-error variables se_u1 and se_u2. The random effects are study-specific deviations from the overall mean effect size. For example, for study 1 and outcome y1, the predicted mean improvement in probing depth is about 0.05 mm higher than the overall mean improvement in probing depth, \( \hat{\theta}_1 = \) 0.357.

For more detail about other postestimation tools, see [META] meta mvregress postestimation.


Antczak‐Bouckoms, A., K. Joshipura, E. Burdick, and J. F. Camilla Tulloch. 1993. Meta‐analysis of surgical versus non‐surgical methods of treatment for periodontal disease. Journal of Clinical Periodontology 20: 259–268.

Berkey, C. S., D. C. Hoaglin, A. Antczak‐Bouckoms, F. Mosteller, and G. A. Colditz. 1998. Meta‐analysis of multiple outcomes by regression with random effects. Statistics in Medicine 17: 2537–2550.

Additional resources

Learn more about Stata's meta-analysis features.

Read more about meta-analysis in the Stata Meta-Analysis Reference Manual.

See Tour of meta-analysis commands in [META] meta.





The Stata Blog: Not Elsewhere Classified Find us on Facebook Follow us on Twitter LinkedIn YouTube Instagram
© Copyright 1996–2023 StataCorp LLC   •   Terms of use   •   Privacy   •   Contact us