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.

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.

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) .describeContains 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 | ||

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

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] | ||

y1 | ||

_cons | .3534282 .0588486 6.01 0.000 .238087 .4687694 | |

y2 | ||

_cons | -.3392152 .0879051 -3.86 0.000 -.5115061 -.1669243 | |

Random-effects parameters | Estimate | |

Unstructured: | ||

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.

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] | ||

y1 | ||

pubyear | .0048615 .0218511 0.22 0.824 -.0379658 .0476888 | |

_cons | .3587569 .07345 4.88 0.000 .2147975 .5027163 | |

y2 | ||

pubyear | -.0115367 .0299635 -0.39 0.700 -.070264 .0471907 | |

_cons | -.3357368 .0979979 -3.43 0.001 -.5278091 -.1436645 | |

Random-effects parameters | Estimate | |

Unstructured: | ||

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.

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] | ||

y1 | ||

_cons | .3572553 .0499616 7.15 0.000 .2593323 .4551782 | |

y2 | ||

_cons | -.3538886 .0788344 -4.49 0.000 -.5084013 -.199376 | |

Random-effects parameters | Estimate | |

Independent: | ||

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**.

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 heterogeneityMethod: Cochran Joint: I2 (%) = 93.76 H2 = 16.03 Method: Jackson–White–Riley y1: I2 (%) = 67.29 R = 1.75 y2: I2 (%) = 94.40 R = 4.23 Joint: 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**.

.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.

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**.