## In the spotlight: Multivariate meta-analysis

Multivariate meta-analysis is the science of jointly synthesizing correlated effect sizes from multiple studies addressing a common research question. For example, in hypertension trials, effect sizes for both systolic blood pressure and diastolic blood pressure are measured on the same set of patients. Because each subject contributes two observations within a study, the two effect sizes are correlated and any meta-analysis that aims at pooling these effect sizes across studies must account for this within-study correlation. Traditionally, the dependence among the effect sizes has been ignored and a separate univariate meta-analysis was conducted for each effect size. This typically leads to biased pooled effects and overestimated standard errors.

Stata 17 has added a new command, **meta mvregress**, to the
**meta** suite to perform multivariate meta-analysis and
meta-regression.

## Teacher burnout data

Many studies suggest a strong link between teacher burnout and teacher attrition. Understanding the driving factors behind teacher burnout has therefore been a main research topic in educational psychology. Consider the meta-analysis conducted by Aloe, Amo, and Shanahan (2014) to investigate the effect of classroom management self-efficacy (CMSE) on teacher burnout. The burnout is measured by the Maslach Burnout Inventory: Educators Survey (MBI-ES). The MBI-ES is an instrument based on three variables: emotional exhaustion (EE), depersonalization (DP), and lowered personal accomplishment (LPA).

According to Aloe, Amo, and Shanahan (2014), ''CMSE is the extent to which a teacher feels that they are capable of gaining and maintaining students’ attention, and dealing with disruption and misbehaving students (Emmer 1990; Emmer and Hickman 1991)'' and can be seen as a protective factor against burnout.

Each of the 16 studies in the meta-analysis reported 3 (correlated) effect
sizes that are stored in variables **y1**, **y2**, and **y3**. These
variables represent the correlation between the CMSE and each of the three
dimensions of burnout.

.describeContains data from aloe.dta Observations: 16 Variables: 19 14 Oct 2021 00:01

Variable Storage Display Value name type format label Variable label | ||

study str26 %26s Study y1 float %9.0g Emotional exhaustion y2 float %9.0g Depersonalization y3 float %9.0g Lowered personal accomplishment v11 float %9.0g Variance of y1 v12 float %9.0g Covariance between y1 and y2 v13 float %9.0g Covariance between y1 and y3 v22 float %9.0g Variance of y2 v23 float %9.0g Covariance between y2 and y3 v33 float %9.0g Variance of y3 country str11 %11s Country of origin female_p float %9.0g Percentage of females experience float %9.0g Years of experience s1 double %10.0g Std. err. of y1 s2 double %10.0g Std. err. of y2 s3 double %10.0g Std. err. of y3 pub_type long %12.0g pubtype Publication type school long %14.0g schoollvl School level pub_year int %8.0g Publication year - 2000 | ||

Traditionally, a separate (univariate) meta-analysis is performed for each
outcome. Below, we use univariate analysis as a descriptive tool to gain
some insight about each dimension of burnout separately.
Frames
can facilitate the task of performing multiple univariate meta-analyses, and
their usage with the **meta** suite is illustrated next.

We start by copying the dataset into three frames, **fr_1**, **fr_2**,
and **fr_3**, where **fr_***i* corresponds to outcome **y***i*,
where *i*=1, 2, or 3. Then, within each frame, we declare our data
(corresponding with one of the dimensions of burnout) as **meta** data
by using the **meta set** command.

.forvalues i = 1/3 { frame copy `c(frame)' fr_`i' frame fr_`i': quietly meta set y`i' s`i', studylabel(study) }

By using frames, we avoid having to specify **meta set**
each time we want to analyze a different outcome. Rather, we
have three meta-analyses already declared in separate frames corresponding
to each outcome of interest. Subsequently, for any outcome, a meta-analysis
technique may be performed by specifying the appropriate command within the
corresponding frame. For example, we could construct a forest plot for each
outcome by specifying the **meta forestplot** command within each frame.
The forest plot in frame **fr_1** corresponds to effect sizes stored in
variable **y1**, and so on.

.local forest_opts nullrefline esrefline.frame fr_1: meta forestplot, name(fp1, replace) `forest_opts'.frame fr_2: meta forestplot, name(fp2, replace) `forest_opts'.frame fr_3: meta forestplot, name(fp3, replace) `forest_opts'

The overall (mean) correlations between CMSE and EE, DP, and LPA are –0.28, –0.33, and 0.43, respectively. The negative correlation between CMSE and both EE and DP suggests that a positive classroom experience leads to less EE and DP feelings. The scenario is reversed between CMSE and LPA (they are positively correlated), suggesting that a positive classroom experience will also lead to feeling more accomplished. From the plots, it seems there is considerable heterogeneity among studies, especially in the correlation between CMSE and the third dimension of burnout, LPA, where 10 out of the 16 studies had CIs that did not cross the overall effect vertical line.

You could also conduct univariate meta-regression of, say, **y2**
(correlation between CMSE and DP) on moderators such as percentage of
female teachers (**female_p**) and publication type (**pub_type**)
by typing

.frame fr_2: meta regress female_p i.pub_type(output omitted)

Or perhaps construct a funnel plot to assess the small-study effect for variable
**y3**:

.frame fr_3: meta funnelplot(output omitted)

## Performing a multivariate meta-analysis

The three separate meta-analyses do not account for the dependence between
**y1**, **y2**, and **y3**. We will use the **meta mvregress**
command to conduct a multivariate meta-analysis to pool the
three-dimensional effect of burnout across studies while accounting
for the dependence among the three outcomes. There are two
components that contribute to the correlation among the effect sizes:
the within-study correlation, which is accounted for by the six
variables **v11**, **v12**, **v13**, **v22**, **v23**, and
**v33**; and the between-study correlation, which is accounted for by
the random-effects model (more specifically, by the covariance matrix of
the random effects) that is assumed by default.

.meta mvregress y1 y2 y3, wcovvariables(v11 v12 v13 v22 v23 v33)Performing EM optimization ... Performing gradient-based optimization: Iteration 0: log restricted-likelihood = 37.738248 Iteration 1: log restricted-likelihood = 40.417156 Iteration 2: log restricted-likelihood = 40.56626 Iteration 3: log restricted-likelihood = 40.567178 Iteration 4: log restricted-likelihood = 40.567179 Multivariate random-effects meta-analysis Number of obs = 48 Method: REML Number of studies = 16 Obs per study: min = 3 avg = 3.0 max = 3 Wald chi2(0) = . Log restricted-likelihood = 40.567179 Prob > chi2 = .

Coefficient Std. err. z P>|z| [95% conf. interval] | ||

y1 | ||

_cons | -.2773755 .0302203 -9.18 0.000 -.3366062 -.2181447 | |

y2 | ||

_cons | -.3285461 .0284353 -11.55 0.000 -.3842783 -.272814 | |

y3 | ||

_cons | .4331267 .0449928 9.63 0.000 .3449425 .521311 | |

Random-effects parameters | Estimate | |

Unstructured: | ||

sd(y1) | .1066603 | |

sd(y2) | .0992441 | |

sd(y3) | .1708421 | |

corr(y1,y2) | .8779874 | |

corr(y1,y3) | -.9685017 | |

corr(y2,y3) | -.876861 | |

More conveniently, you may express the above specification by using the stub shortcut notation:

.meta mvregress y*, wcovvariables(v*)(output omitted)

The first table displays the overall trivariate correlation from the multivariate model. The estimates are very similar to the univariate ones reported on the forest plots. This is often the case, especially when no missing data are reported by the studies (which is our case). Whenever studies report missing values on some outcome of interest, multivariate meta-analysis can, by observing related outcomes, borrow strength from these outcomes to produce more efficient overall estimates for all outcomes. The three overall correlations are significant, suggesting that EE, DP, and LPA all relate to CMSE.

The multivariate homogeneity test, which tests whether the trivariate
correlation is constant across studies, is rejected (*p* < 0.00001).
The second table shows the standard deviations and correlations of
the random effects.

The default model assumed by **meta mvregress** is a random-effects
model with an unstructured covariance matrix for the random effects
and is fit via restricted maximum-likelihood estimation (REML).
You may use the **random()** option to request other
settings.
For example, below we specify a random-effects model with an exchangeable
covariance structure for the random effects estimated via maximum likelihood
estimation (MLE).

.meta mvregress y*, wcovvariables(v*) random(mle, covariance(exchangeable))Performing EM optimization ... Performing gradient-based optimization: Iteration 0: log likelihood = 35.789919 Iteration 1: log likelihood = 36.294019 Iteration 2: log likelihood = 36.298873 Iteration 3: log likelihood = 36.298874 Multivariate random-effects meta-analysis Number of obs = 48 Method: ML Number of studies = 16 Obs per study: min = 3 avg = 3.0 max = 3 Wald chi2(0) = . Log likelihood = 36.298874 Prob > chi2 = .

Coefficient Std. err. z P>|z| [95% conf. interval] | ||

y1 | ||

_cons | -.2812558 .0312227 -9.01 0.000 -.3424512 -.2200604 | |

y2 | ||

_cons | -.2812558 .0312227 -9.01 0.000 -.3424512 -.2200604 | |

y3 | ||

_cons | .4406485 .0308653 14.28 0.000 .3801536 .5011434 | |

Random-effects parameters | Estimate | |

Exchangeable: | ||

sd(y1 y2 y3) | .1103228 | |

corr(y1 y2 y3) | -.4598057 | |

The first table shows similar values to the overall trivariate correlation from the previous REML model. The second table displays the assumed common value of the random-effects standard deviation (0.1103) and the assumed common correlation among the random effects (–0.4598).

You could also specify **random(jwriley)** to request the
Jackson—White—Riley noniterative method for estimating the
between-study random-effects covariance matrix. The method can be seen
as a multivariate extension of the popular univariate method of
DerSimonian and Laird.

If you wish to assume a fixed-effects multivariate MA model, you may
specify option **fixed**:

.meta mvregress y*, wcovvariables(v*) fixed(output omitted)

## Assessing multivariate heterogeneity

Upon successfully fitting your multivariate model, you may explore the
heterogeneity among the correlated effect sizes across the studies via
the **estat heterogeneity** postestimation command:

.quietly meta mvregress y*, wcovvariables(v*).estat heterogeneityMethod: Cochran Joint: I2 (%) = 82.47 H2 = 5.71 Method: Jackson–White–Riley y1: I2 (%) = 82.37 R = 2.38 y2: I2 (%) = 83.92 R = 2.49 y3: I2 (%) = 94.52 R = 4.27 Joint: I2 (%) = 77.06 R = 2.09

The Cochran heterogeneity statistics quantify the amount of
heterogeneity jointly for all outcomes. \( I^2\) = 82.47% means that
roughly 82% of the variability among the effect sizes across the
studies is due to true heterogeneity between the studies as opposed to
the sampling variability. The Jackson—White—Riley statistics quantify
the amount of heterogeneity separately for each outcome variable and
jointly for all outcomes. For instance, we saw from the forest plot of
variable **y3** (the correlation between CMSE and LPA) that there was
considerable heterogeneity because of the many nonoverlapping study CIs.
This is confirmed by the highest \( I^2\) value (roughly 95%)
corresponding to outcome **y3**.

.estat heterogeneity, jwriley(y1 y2)Method: Jackson–White–Riley y1 y2: I2 (%) = 76.23 R = 2.05

The value of \( I^2\) is roughly 76%, which suggests that, when considered as a
pair, the effect sizes in variables **y1** and **y2** are less variable
than the individual effect sizes, albeit with considerable heterogeneity
still present.

## Predicting random effects and more

Various standard postestimation tools—such as
predict,
margins, and
test—are also
available after you fit a multivariate meta-analysis
model. Below, we predict random effects by using **predict,
reffects**, and we obtain their standard errors
by specifying the **reses()** option.

.predict double re*, reffects reses(re_se*).set linesize 90.list re*

re1 re2 re3 re_se1 re_se2 re_se3 | |

1. | -.10283588 -.03131996 .17620617 .09978001 .0901633 .16225359 |

2. | -.09379568 -.05495025 .12575747 .09988817 .09205702 .16277886 |

3. | -.02938791 -.01429996 .0458616 .09973658 .09082035 .16198752 |

4. | -.00674903 -.04699805 -.02200034 .098104 .08902432 .16011173 |

5. | -.15505505 -.14830599 .25088819 .09530048 .08396402 .15573488 |

6. | .02059561 .01613239 -.03490522 .09268252 .08009149 .15180802 |

7. | -.05697903 -.02027662 .12091516 .09937867 .09083816 .16232895 |

8. | .09958159 .04976988 -.16574508 .09476274 .08488424 .15428564 |

9. | .18321031 .19925303 -.29725953 .09895307 .09032536 .16101285 |

10. | .1637504 .14524732 -.27103958 .08737944 .07581568 .14315724 |

11. | .1196294 .08084109 -.18393054 .0898793 .07911838 .14714454 |

12. | -.03289968 -.0253095 .07421782 .09507486 .08363428 .15483541 |

13. | -.00577343 -.07779175 .04872011 .09989025 .09205789 .16304962 |

14. | -.12152979 -.11379705 .18690794 .09933887 .08977333 .1614794 |

15. | -.00072944 .00359832 -.02299691 .09546701 .08470254 .15506071 |

16. | .01896759 .0382071 -.03159726 .09648813 .08623672 .15749337 |

Each row of predicted random effects represents the study-specific
deviations from the mean trivariate correlation. Instead of the
**reses()** option, you could have specified the **revce()**
option to compute the full variance—covariance matrix of the
predicted random effects in lieu of standard errors only.

## Fitting a multivariate meta-regression (moderator analysis)

We can try to explain some of the heterogeneity highlighted in the
univariate forest plots but also in the multivariate homogeneity test
by performing a multivariate meta-regression using the variable
**pub_year** (centered around year 2000). Perhaps newer studies report
different magnitudes of the three outcomes compared with older studies,
which might explain the discrepancies between the trivariate correlations
across the studies. We simply incorporate our moderator by specifying
it on the right-hand side of the equals sign, as follows:

.set linesize 80.meta mvregress y* = pub_year, wcovvariables(v*) variancePerforming EM optimization ... Performing gradient-based optimization: Iteration 0: log restricted-likelihood = 26.721723 Iteration 1: log restricted-likelihood = 28.24878 Iteration 2: log restricted-likelihood = 28.36236 Iteration 3: log restricted-likelihood = 28.362971 Iteration 4: log restricted-likelihood = 28.362971 Multivariate random-effects meta-regression Number of obs = 48 Method: REML Number of studies = 16 Obs per study: min = 3 avg = 3.0 max = 3 Wald chi2(3) = 7.36 Log restricted-likelihood = 28.362971 Prob > chi2 = 0.0613

Coefficient Std. err. z P>|z| [95% conf. interval] | ||

y1 | ||

pub_year | -.0050658 .0044898 -1.13 0.259 -.0138657 .003734 | |

_cons | -.2495888 .0382269 -6.53 0.000 -.3245121 -.1746654 | |

y2 | ||

pub_year | -.0079188 .0034126 -2.32 0.020 -.0146074 -.0012302 | |

_cons | -.2903522 .0292901 -9.91 0.000 -.3477597 -.2329446 | |

y3 | ||

pub_year | .0085855 .0066232 1.30 0.195 -.0043957 .0215666 | |

_cons | .3873665 .056148 6.90 0.000 .2773184 .4974147 | |

Random-effects parameters | Estimate | |

Unstructured: | ||

var(y1) | .010966 | |

var(y2) | .0058382 | |

var(y3) | .0267623 | |

cov(y1,y2) | .007134 | |

cov(y1,y3) | -.0164038 | |

cov(y2,y3) | -.0107893 | |

We used option **variance** to report variances and covariances of random
effects instead of standard deviations and correlations. The slope of
**pub_year** for outcome **y2** (–0.008) is significant and
suggests that newer studies report smaller correlations between CMSE
and DP. If we compare the random-effects table with the one from
the constant-only model above, it seems that the moderator explained some
of the between-studies variability for outcome **y2**: **var(y2)** went
from 0.0098 (in the constant-only model) to 0.0058. However, the
residual homogeneity test (*p* < 0.0001) suggests there is still a
considerable amount of residual heterogeneity (remaining heterogeneity
not explained by moderators). This is also confirmed if we type

.estat heterogeneityMethod: Cochran Joint: I2 (%) = 79.39 H2 = 4.85 Method: Jackson–White–Riley y1: I2 (%) = 83.35 R = 2.45 y2: I2 (%) = 77.81 R = 2.12 y3: I2 (%) = 94.42 R = 4.23 Joint: I2 (%) = 75.67 R = 2.03

Notice that there was a reduction (from 83.92% for the constant-only model to
77.81%) in the Jackson—White—Riley \( I^2\) statistic
corresponding to outcome **y2** but the other values of the \( I^2\)
statistic corresponding to **y1** and **y3** remained nearly
the same.

## References

Aloe, A. M., L. C. Amo, and M. E. Shanahan. 2014. Classroom management
self-efficacy and burnout: A multivariate meta-analysis. *Educational
Psychology Review* 26: 101–126.

Emmer, E. T. 1990. A scale for measuring teacher efficacy in classroom management and discipline. Paper presented at the annual meeting of the American Educational Research Association, Boston, MA.

Emmer, E. T., and J. Hickman. 1991. Teacher efficacy in classroom management
and discipline. *Educational and Psychological Measurement* 51: 755–765.
doi:10.1177/0013164491513027.

— by Houssein Assaad

Principal Statistician and Software Developer