*! version 3.1 13jul2015

Demonstration for the mvmeta package in Stata: running mvmeta
Ian White

How to run this tutorial

You will need to install the mvmeta package and get the associated data and graph files.

Installing mvmeta

Click the following command

net from http://www.homepages.ucl.ac.uk/~rmjwiww/stata/

then click on meta, then on mvmeta, then on (click here to install). This installs the package.

Now click on (click here to return to the previous screen) and (click here to get).
This installs the associated data and graph files.

Getting the data in

A separate demonstration is available.

Berkey data

Data from Berkey et al (1998): treatment effects on
two outcomes (probing depth, y1; attachment level, y2)
in periodontal disease.
The within-trial variances and covariances were reported by the authors.

use berkey, clear

l, noo sep(0)

Now we'll draw a bubble plot of the data. This is an option of mvmeta.
We will use the noestimate option to avoid fitting the model yet:

mvmeta y V, noestimate bubble

We fit a bivariate meta-analysis:

mvmeta y V

The y1 results are very similar to those from univariate meta-analysis:

mvmeta y V, var(y1)

We compare the univariate mvmeta results with those from metan.

gen s1 = sqrt(V11)

metan y1 s1

These differ because by default (1) mvmeta fits the random-effects model,
whereas metan fits the fixed-effect model;
(2) mvmeta uses REML,
whereas metan uses the method of moments;
(3) mvmeta's standard error allows for uncertainty in estimating tau.
To get exact agreement:

metan y1 s1, random

mvmeta y V, var(y1) mm nouncertainv print(bscov)

Let's return to the bivariate setting and
explore a meta-regression on year of publication.

mvmeta y V pubyear

There's no evidence that either outcome is associated with year of publication.

p53 data

These data come from 6 observational studies in patients with squamous cell carcinoma of the oropharynx.
The presence of mutant p53 tumor suppressor gene is considered as a possible prognostic factor (Jackson et al, 2011).

The data are the estimated log hazard ratios (lnHR)
for mutant vs. normal p53 gene for two outcomes,
overall survival (OS) and disease-free survival (DFS),
together with their variances.

use p53, clear

l, noo sep(0) abb(12)

Now we'll draw a bubble plot of the data.
The data don't include the within-study correlations:
we assume they are all 0.7.

mvmeta lnHR VlnHR, bubble noestimate wscorr(0.7)

Note that the OS results in studies without
DFS results are much larger than those in studies with DFS results.
This suggests that the multivariate result for DFS may be substantially larger
than the univariate result - but also that we should be cautious of both results.

Let's fit the univariate meta-analysis for DFS:

mvmeta lnHR VlnHR, var(lnHRdfs)

Let's compare with the multivariate meta-analysis.

mvmeta lnHR VlnHR, wscorr(0.7)

Yes, the multivariate result for DFS is larger (lnHRdfs is less negative)
than the univariate result. It also has a larger between-studies variance,
and hence a larger standard error.

Various post-estimation options can be called after a mvmeta
fit without re-fitting the model.
The i2 option estimates a multivariate I-squared,
together with its confidence interval.
It also gives a confidence interval for the between-studies correlation:

mvmeta, i2

The eform option shows the exponentiated coefficients - here hazard ratios:

mvmeta, eform

We can also see the full parameterisation of the model,
including the Cholesky decomposition of the variance terms
which are usually hidden:

mvmeta, showall

Other post-estimation options include t(#) to specify
a t-distribution for inference.

Above we assumed the unknown within-study correlations were all 0.7.
We can avoid this assumption by using Riley's alternative model (Riley et al, 2008).
Usually this converges quickly,
but for these data on this occasion it can need more than 2000 iterations:

mvmeta lnHR VlnHR, wscorr(riley)

Adding the difficult option leads to convergence in a small number of iterations,
but a warning is printed because the overall correlation is close to 1.
It is better for these data not to rely on the alternative model.

mvmeta lnHR VlnHR, wscorr(riley) difficult

Fibrinogen Studies Collaboration 1

Fully and partly adjusted associations

The original data are from 31 studies relating plasma levels of fibrinogen,
a blood clotting factor, to time to a coronary heart disease (CHD) event (Fibrinogen Studies Collaboration, 2005).
In this example, we assume a linear association between fibrinogen and CHD,
and we wish to adjust for confounding. Some confounders are recorded in all studies, while
others are recorded in only 14 studies.
We therefore estimate a partly adjusted coefficient (log hazard ratio) in all 31 studies,
and a fully adjusted coefficient in the 14 studies.
We also estimate their (within-studies) correlation: in the paper we considered three methods,
but here we'll use the bootstrap method (Fibrinogen Studies Collaboration, 2009).

use fscfpama, clear

l, noo

In the variable names, "fa" refers to fully adjusted estimates, and "pa" refers to partly adjusted estimates.
"corrb" is the correlation between fully and partly adjusted estimates, estimated by the bootstrap method.

In order to do mvmeta we need to construct the variance-covariance matrices:

gen varfafa = sefa^2

gen varpapa = sepa^2

gen varfapa = corrb*sefa*sepa

The "standard" approach would be to analyse only the fully adjusted estimates:

mvmeta beta var, var(betafa)

Our new approach analyses the fully-adjusted estimates jointly with the partly-adjusted estimates
in order to gain precision:

mvmeta beta var

The standard error for betafa has decreased from 0.0389 to 0.0266 or 32%.
This represents a 53% decrease in variance.
Note the between-studies correlation is estimated as 1, so that the model is able to infer
fully-adjusted estimates quite precisely from partly-adjusted estimates.

The wt option estimates the borrowing of strength
- the degree to which results for one outcome gain precision
by the inclusion of the other outcome(s) in the analysis.

mvmeta, wt

The output shows a 53.3% borrowing of strength for the fully adjusted result (using the column headed "Overall: betafa").

The output also shows the relative contributions of the studies to the fully-adjusted result:
for example, study 12 makes the largest contribution (9.8% of the total weight).
A large study which is only partly adjusted may contribute more weight
than a small study which is fully adjusted.
For example, study 28 which is only partly adjusted contributes more weight (5.7% of the total weight)
than study 1 which is fully adjusted (0.5% of the total weight).
This is because study 28 is much more precise:

l cohort sepa if inlist(cohort,1,28)

Fibrinogen Studies Collaboration 2

Shape of exposure-outcome relationship

We now use the same original data to explore the shape of association
between fibrinogen and CHD, adjusting for complete confounders.
Each study has been analysed using a Cox model including fibrinogen
categorised into 5 groups and adjusting for confounders.
The "outcomes" of interest are therefore the 4 contrasts (log hazard ratios) of groups
2-5 with group 1.
Some studies (e.g. study 15)
have no participants or no events in group 1:
these have been handled by introducing ("augmenting") a very small amount of
data in group 1.
mvmeta_make has been used to automate the augmentation, fitting of the Cox models
and extraction of the point estimates, variances and covariances (White, 2009).

use fscstage1, clear

browse

mvmeta b V

estimates store fsc2full

That was a little slow.
We'll demonstrate some faster alternatives.
First and probably best, the method of moments. The original version is by Jackson et al (2010):

mvmeta b V, mm

and a matrix-based method of moments is by Jackson et al (2013):

mvmeta b V, mm2

The fixed-effect method (not recommended, because it ignores heterogeneity):

mvmeta b V, fixed

We could also assume that the between-studies variation
is captured by a random slope:

matrix B = (1,2,3,4)'*(1,2,3,4)

mvmeta b V, bscov(prop B)

lrtest fsc2full

Not significantly worse,
but I'd prefer to use the full model:

estimates replay fsc2full

Finally, here's a graph I made earlier:

graph use fsc2cigraph

References

Berkey CS, Hoaglin DC, Antczak-Bouckoms A, Mosteller F, Colditz GA (1998).
Meta-analysis of multiple outcomes by regression with random effects.
Statistics in Medicine 17: 2537-2550.
Link.

Jackson D, White IR, Thompson SG (2010).
Extending DerSimonian and Laird's methodology to perform multivariate random effects meta-analyses.
Statistics in Medicine 29: 1282-1297.
Link.

Jackson D, Riley R, White IR (2011).
Multivariate meta-analysis: potential and promise.
Statistics in Medicine 30: 2481-2498.
Link.

Jackson D, White IR, Riley R (2013).
A matrix based method of moments for fitting the random effects model
for meta-analysis and meta-regression.
Biometrical Journal 55: 231-245.
Link.

Fibrinogen Studies Collaboration (2005).
Plasma fibrinogen and the risk of major cardiovascular diseases and non-vascular mortality: meta-analysis of individual data for 154 211 adults in 31 prospective studies.
Journal of the American Medical Association 294: 1799-1809.
Link.

Fibrinogen Studies Collaboration (2009).
Systematically missing confounders in individual participant data meta-analysis of observational cohort studies.
Statistics in Medicine 28: 1218-1237.
Writing committee: Jackson D, White IR.
Link.

Riley RD, Thompson JR, Abrams KR (2008).
An alternative model for bivariate random-effects meta-analysis when the within-study correlations are unknown.
Biostatistics 9: 172-186.
Link.

White IR (2009).
Multivariate random-effects meta-analysis.
Stata Journal 9: 40-56.
Link.