*! v4.0.2 21apr2022

Title

mvmeta -- Multivariate random-effects meta-analysis and meta-regression

Description

Standard meta-analysis combines estimates of one parameter over several studies.
Multivariate meta-analysis is an extension that can combine estimates of several related parameters.
For example, we may have estimates of treatment effects on two different outcomes in a set of randomised trials,
or we may have estimates of the difference in outcome between "high", "medium" and "low" levels of a covariate.
For further discussion of the multivariate meta-analysis model and its applications see Jackson et al (2011).

mvmeta performs multivariate random-effects meta-analysis (White, 2009)
and multivariate random-effects meta-regression (White, 2011)
on a data-set of point estimates, variances and (optionally) within-study covariances.

Demonstrations are available on
getting the data into mvmeta format
and running mvmeta.

If you are doing network meta-analysis, please see my network suite,
which includes commands to formulate and run mvmeta models for network meta-analysis.

Syntax

mvmeta b V [xvars]
[if] [in]
[,

vars(varlist)
wscorr(expression|riley)
bscovariance(covtype [matexp])
equations(yvar1:xvars1[, yvar2:xvars2[, ...]])
noconstant
commonparm

reml|ml|mm1|mm2|fixed
start(expression)
longparm
noposdef
psdcrit(#)
maximize_options
suppress([fe] [uv] [mm])

nouncertainv
eform[(name)]
dof(expression)
randfix[(suboptions)]

wt[(suboptions)]

pbest(min|max [if] [in], pbest_options)

print(string)
i2
i2fmt(%fmt)
ciscale(sd|logsd|logh)
testsigma

showstart
id(varname)
showall
warning(error|text|off)
bubble[(suboptions)
qscalar
pi[(suboptions)]]

where the data are arranged with one line per study:
the point estimates are held in variables b*,
the variance of bx is held in variable Vxx,
and the within-study covariance of bx and by is held in variable Vxy.
If you have individual participant data, you can use mvmeta_make
to produce a dataset of the required format.

Model options

vars(varlist) specifies which variables are to be used.
By default, all variables named b* are used, but any variable named b is ignored.
Note that vars(b1 b2 b3) and vars(b3 b2 b1)
specify the same model but a different parameterisation.

wscorr(expression) specifies the value of all within-study correlations.
This means that within-study covariance variables Vxy need not exist.
(For any that do exist, wscorr() is ignored.)
Alternatively, wscorr(riley) can be used when within-study correlations are unknown.
It uses the alternative model of Riley et al (2008) to estimate an overall correlation.
Riley (2009) discusses other ways to handle unknown within-study correlations.

bscovariance(covtype)
specifies the between-studies covariance structure:
see Between-studies covariance structures below.

equations(yvars1:xvars1 [,yvars2:xvars2[,...]])
allows different outcomes to have different regression models.
For example, equations(y1 y2:x1, y3:x2) regresses y1 and y2 on x1 and y3 on x2.

noconstant suppresses the constant in the regression models.

commonparm forces the coefficients to be the same for all outcomes.
This option requires all equations to contain the same number of variables.
It can be useful in network meta-analysis and is illustrated in White et al (2012).

Estimation options

reml specifies that restricted maximum likelihood is to be used for estimation (the default).

ml specifies that maximum likelihood is to be used for estimation.
This is likely to underestimate the variance, so restricted maximum likelihood is usually preferred.

mm1 specifies that the multivariate method of moments procedure of
Jackson et al, 2010 is to be used for estimation.
This is a multivariate generalisation of the procedure of
DerSimonian and Laird (1986) and is faster than the likelihood-based methods.
However, the mm2 option is to be preferred.

mm2 specifies that the multivariate method of moments procedure of
Jackson et al, 2013 is to be used for estimation.
This is a matrix-based extension of the method given by the mm1 option and is to be preferred.
This method is used if mm is specified.

fixed specifies that the fixed-effect or common-effect model
(that is, the model allowing no heterogeneity between studies)
is to be used for estimation.

start(matrix|matrix expression|mm|#)
specifies a starting value for the between-studies variance.
The syntax depends on the between-studies covariance structure:
see between-studies covariance structures below.

longparm parameterises the model as one regression model for each outcome.
Without covariates, this is usually less convenient than the default
(all outcomes forming a single regression model), but is required if the pbest() option will be used.
With covariates, longparm is the default and cannot be changed.

noposdef does not immediately halt execution if a within-study variance-covariance matrix is declared not to be positive semidefinite. This option suppresses the fixed-effect analysis.

psdcrit(#) changes the criterion for judging a variance-covariance matrix not to be positive semidefinite.
A variance-covariance matrix is regarded as positive semidefinite if the ratio of the smallest to the largest ratio is greater than the negative of #.
Default is 1E-8.

maximize_options allows standard maximization options: see [R] maximize.

suppress([fe] [uv] [mm]) suppresses one of more of the following analyses
which are run by default:
a fixed-effect analysis (used by the randfix() option),
univariate analyses (used by the wt(rv) option),
and a method-of-moments analysis with unstructured between-studies covariance matrix (used by the i2 option).
This option is useful if these analyses fail.

Output options: fixed parameters

nouncertainv computes alternative (smaller) standard errors that ignore
the uncertainty in the estimated between-studies variance-covariance matrix, and therefore
agree with results produced by procedures such as SAS PROC MIXED (without the
ddfm=kr option) and metareg. (Note however that confidence intervals do not
agree, because mvmeta by default uses a Normal approximation whereas the other procedures
approximate the degrees of freedom of a t distribution.)

eform or eform(name) exponentiates the reported mean parameters and,
if the optional argument is used, reports them as name.

dof(expression) specifies the degrees of freedom for t-tests and confidence intervals.
The expression may include "n", the number of observations. The default is to use a normal distribution.

randfix or randfix(varlist) describes the impact of heterogeneity on the estimated coefficients in the models for the specified variables (which must be outcome variables in the fitted model).
If no varlist is specified then all outcome variables are used.

The method is described by Jackson et al (2012).
The estimated variance-covariance matrix of the parameters is compared with that estimated under a fixed-effect model:
the square roots of the determinants of these matrices are reported, followed by R, defined as the qth root of the ratio of these square roots, where q is the number of parameters involved.
R is a multivariate generalisation of the R statistic of Higgins and Thompson (2002), who note that R is often approximately equal to H, where I^2=(H^2-1)/H^2.

Output options: weights and borrowing of strength

The option wt or wt(suboptions)
reports study weights and/or borrowing of strength.
The methods are described in Jackson et al (2017).
The suboptions are:

sd (the default) reports study weights and borrowing of strength
using the score decomposition method.
The method computes, for each parameter, the percentage of the total information for that parameter
by study and source,
where source is direct (from that outcome), borrowed (from other outcomes) and total.
By default, the total for each study across sources
and the total for each source across studies are reported.

rv reports borrowing of strength using the relative variances method.
This compares the variance of each coefficient in the multivariate meta-analysis
with its variance in a univariate meta-analysis.
Unlike the sd method, the rv method takes account of changes
in the estimated between-studies heterogeneity matrix,
and hence may report negative values for borrowing of strength.

dpc reports data point coefficients. These show how each estimated parameter
is derived as a weighted sum of data points.

details outputs a table of the full score decomposition (sd method only)
or of the separate standard errors (rv method only).

nostudies suppresses listing of the study-specific weights (sd method only).

format(fmt) specifies the output format for all methods.

clear (sd method only) loads the data for the table into memory.

keepmat(name) saves the matrices.
With the sd method, three matrices,
nameborrowed, namedirect and nametotal, are saved.
For rv & dpc methods, one matrix, name, is saved.

unscaled (sd method only)
modifies the details and keepmat(name) options
to give unscaled, not scaled weights.

wide with the sd and details options causes the output to be in wide format.

Output options: probability-of-being-best

pbest(min|max [if] [in], [suboptions])
is useful in network meta-analysis.
It requests estimation of the probability that each linear predictor is the best
-- that is, the maximum or minimum, depending on the first argument of pbest().
Estimation is performed for each record in the current data that satisfies the if and in criteria.

The probability is estimated under a Bayesian model with flat priors,
assuming that the posterior distribution of the parameter estimates is approximated
by a Normal distribution with mean and variance
equal to the frequentist estimates and variance-covariance matrix of the parameters.
Rankings are constructed by drawing the coefficients multiple times
from their approximate posterior density.
For each draw, the linear predictor is evaluated for each study,
and the largest linear predictor is noted.

For models without covariates, pbest() is only available if longparm was specified when the model was fitted.

The suboptions of pbest() are:

reps(#) specifies the number of draws used. The default is 1000.

zero specifies that zero is to be considered as another linear predictor.

gen(string) specifies that the probabilities be saved in variables with prefix string.

seed(#) specifies the random number seed.

format(%fmt) specifies the output format.

id(varname) specifies an identifier for the output.

predict ranks the true effects in a future study with the same covariates,
thus allowing for heterogeneity as well as parameter uncertainty,
as in the calculation of prediction intervals (Higgins et al, 2009).
The default behaviour is instead to rank linear predictors and does not allow for heterogeneity.

bestonly requests only the probabilities of being the best treatment.
This can be misleading (Salanti et al, 2011)).
The default is to report the probabilities for all ranks.

saving(filename) writes the draws from the posterior distribution
(indexed by the identifier and the replication number) to filename,
and replace allows this file to be overwritten.

clear causes the summarised probabilities to be loaded into memory,
so that the user can produce their own tabulations or graphs. This disables the gen() option.

bar draws a bar graph of the probabilities.

line draws a line graph of the probabilities.

cumulative changes the bar or line graph to show cumualtive probabilties.
The rankogram of Salanti et al (2011)
is produced by specifying this together with the line suboption.

mcse adds the Monte Carlo standard errors to the tables.

meanrank adds the mean rank and the SUCRA (Salanti et al, 2011)
to the table.
The SUCRA is the rescaled mean rank: it
is 1 when a treatment is certain to be the best and 0 when a treatment is certain to be the worst.

tabdispoptions(string) specifies options valid for tabdisp which draws the results table.
For example, tabdispoptions(cellwidth(10)).

Output options: variance parameters

print(string) determines how the between-studies variance-covariance matrix is reported.
print(bscorr), the default, reports the between-studies standard deviations and correlation matrix.
print(bscov) reports the between-studies variance-covariance matrix.
print(bscov bscorr) reports both.

i2 reports the between-study variance and the I-squared statistic for each outcome,
together with their confidence intervals.
I-squared is computed as between / (within + between)
where "between" is the appropriate element of Sigma
and "within" is given by equation (9) of Higgins and Thompson (2002).
The method is described in section 3.6 of White (2011).

Confidence intervals are computed as follows.
With fixed or mm1 estimation, the method of Higgins and Thompson (2002) is used,
implemented by heterogi.
With reml or ml estimation, confidence intervals are computed by nlcom
on a scale specified by the ciscale option.
With mm2 estimation, I have implemented a new and even more ad hoc method:
the estimated between-study variance is used to reconstruct the diagonal elements of the Q matrix,
and the method of Higgins and Thompson (2002)
is used as if this were the observed Q matrix.

i2fmt(%fmt) specifies the output format for the I-squared statistics.

ncchi2 uses this option of heterogi in computing confidence intervals. It is only relevant after mm estimation.

ciscale(sd|logsd|logh) determines the scale on which confidence intervals for the between-study variance and I-squared are computed after reml or ml estimation:
tau, log(tau) or log(H), where tau^2 is the "between" variance above and H^2 = (1 + between / within).
The default is ciscale(sd).

testsigma performs a likelihood ratio test of Sigma=0, if reml or ml estimation was used.

qscalar performs a global test for heterogeneity using the scalar Q statistic.

pi or pi(suboptions) reports prediction intervals for the mean of each outcome. Confidence intervals are also reported. Suboptions are

format(string) specifies the numerical format of the values.

level(#) specifies the confidence level to be used.

xvar(varname) specifies the variable in a meta-regression whose coefficient is to be used.
By default the intercept (constant) is used.
This option is only useful when the constant is represented by a term other than _cons.

Output options: bubble plot

bubble or bubble(suboptions) draws a bubble plot of the data for two of the outcomes of the last fitted mvmeta command.
The two outcomes are plotted on the two axes and each study is represented by the point estimate and a confidence ellipse reflecting the within-study variances.
The univariate and multivariate pooled estimates are similarly plotted.
Suboptions are:

variables(varlist) specifies the two outcome variables to be plotted. The default is the first two outcome variables.

nomv omits the univariate pooled estimates from the bubble plot.

nouv omits the multivariate pooled estimates from the bubble plot.

studylabel(text) specifies the label to appear in the legend for the study estimates.

uvlabel(text) specifies the label to appear in the legend for the univariate pooled estimates.

mvlabel(text) specifies the label to appear in the legend for the multivariate pooled estimates.

clear replaces the data in memory with the data set used to draw the bubble plot. The bubble plot command may be retrieved with F9 and may be edited.

mlabel labels each point estimate with its identifier.

pct(numlist) specifies one or more percentages for confidence ellipses. The default is pct(50).

n(#) specifies the number of points to be plotted in each confidence ellipse. The default is 36.

missval(# [#]) specifies the axis values at which missing values are plotted.
If one value is given, it is for both axes.
If two values are given, they are for the x and y axes respectively.
The default is missval(0).

Whereas a study with both outcomes observed is represented by confidence ellipse,
a study with a missing value on one outcome is represented by a confidence interval for the non-missing outcome.
A study with missing values on both outcomes is omitted from the bubble plot.

stagger(#) specifies that the axis values in missval() should be separated by an amount # if there are multiple studies with missing values.

colors(text) specifies the colours used for the studies, univariate summary and multivariate summary respectively.

Other graph options are likely to work, e.g. lpatterns(dash).

Output options: miscellaneous

showstart reports the starting values used.

id(varname) specifies an identifier for the output.
This affects wt(), pbest() and bubble(). The default identifier is the observation number.

showall reports the estimated values of the basic parameters underlying the between-studies
variance matrix.

warning(error|text|off) controls whether warning messages are shown as errors, shown as text, or not shown.
warning(error) is the default, and is strongly recommended.

Between-studies covariance structures

The between-studies variance-covariance matrix Sigma may be modelled in various ways.
Each option has a different way to specify the starting values for Sigma.
In each case, the starting value for the fixed parameters is derived from the starting value of Sigma.

bscovariance(unstructured) estimates an unstructured Sigma, and is the default.
Starting values for Sigma may be specified explicitly by start(matrix_expression).
start(mm) (the default) specifies that the starting value is computed by the mm method.
start(0) uses a starting value of 0.001 times the default.

bscovariance(proportional matexp) models Sigma = tau^2*matexp, where tau is an unknown parameter and matexp is a known matrix expression (e.g. a matrix name or I(2)).
start(#) specifies the starting value for the scalar tau.

bscovariance(exchangeable #) is a shorthand for bscovariance(proportional P)
where P is a matrix with 1's on the diagonal and # off the diagonal.
bscovariance(exchangeable 0.5) is widely used in network meta-analysis.

bscovariance(exchangeable) models Sigma = tau^2*P,
where P is a matrix with 1's on the diagonal and a constant rho off the diagonal.
start(# #) specifies the starting values for tau and rho.

bscovariance(equals matexp) forces Sigma = matexp, where is a known matrix expression (e.g. a matrix name or I(2)).
start() is not relevant.

bscovariance(correlation matexp) models Sigma = D*matexp*D,
where matexp is a known matrix expression containing the between-study correlations,
and D is an unknown diagonal matrix containing the between-studies standard deviations.
start(rowvector) specifies the starting values for the diagonal of D.

Studies in which some outcomes are unestimated

mvmeta now deals naturally with cases where a study reports only a subset of outcomes:
that is, all computation methods are adapted to handle this case.
mvmeta ignores within-study variances and covariances specified for missing point estimates.
Conversely, it expects non-missing within-study variances and covariances to accompany non-missing point estimates.

mvmeta_make automatically fills in missing values using the augmentation algorithm described in White (2009).

Changes in mvmeta version 2.3

mvmeta version 2.3 was published in the SJ (White, 2011) with the following changes:

Meta-regression is allowed.
The simple syntax is mvmeta b V xvars.
The more flexible syntax uses the eq() option.
For example, for 2-dimensional b, mvmeta b V x is the same as mvmeta b V, eq(b1:x,b2:x).

mvmeta, typed without specifying b and V, redisplays the latest estimation results.
All the output options listed above may be used, except keepmat().

eform is correctly implemented, and is ignored if long parameterisation is used without covariates.

The starting values (for bscov(uns)) are produced by default by the method of moments.

Option showchol has been renamed showall.

The likelihood is coded using mata and appears on initial tests to be 2-5 times faster.

corr() has been renamed wscorr().

bscorr and bscov have been renamed print(bscorr) and print(bscov).

Changes from version 2.3 to version 4.0

The command has been modified to work with the new network suite for network meta-analysis.
Particular changes include the new commonparm option, and the shorthand bscov(exchangeable #).

The equations() option now allows lists of y-variables.

The matrix-based method of moments has been added as the mm2 option
and is the default if mm is specified.

The following new suboptions for pbest allow various graphical displays, saving of results and numerical summaries of the estimated ranks:
all
saving(filename)
clear
bar
line
cumulative
mcse
meanrank
tabdispoptions(string)

The new wt option displays study weights and borrowing of strength.

The new randfix option compares random-effects with fixed-effect results.

The program structure has been changed: by default, the fixed-effect model, all univariate models and the unstructured method of moments are fitted before the specified model.
The suppress() option can be used to suppress some or all of these analyses; in particular, it is used by network.

The augmentation procedures used by version 1 of mvmeta are unnecessary and are now undocumented.
The options augment, augquiet, missest(#) and missvar(#) remain available.

A bug in the estimation procedure for the wscorr(riley) method which led to wrong answers has been fixed.

The bubble and pi options have been added.

The bscov(exch) option has been added, with the ability to estimate the between-studies correlation.

Constrained estimation can now be done: see constraint.

A number of minor bugs have been fixed, and the output improved.

Examples

First stage, starting with individual participant data (fg has levels 1-5):

. xi: mvmeta_make stcox ages i.fg, strata(sex tr) nohr saving(fscstage1) replace by(cohort) usevars(i.fg) names(b V) esave(N)

The individual participant data are not publicly available, but you can get the summary data produced by the above command using

. net get st0156, from(http://www.stata-journal.com/software/sj9-1)

Second stage:

. use fscstage1, clear

. mvmeta b V

For more examples, please see the demonstration.

Details

The reml and ml methods use Newton-Raphson maximisation of the likelihood
or restricted likelihood using
ml.
The between-studies variance matrix (using bscov(uns)) is parameterised
via its Cholesky decomposition in
order to ensure that it is non-negative definite.

Parts of mvmeta require the additional program sencode.

mvmeta aims to work under all Stata versions 9 and later, but I have been unable to test it on versions 9-11.

Known problems

With methods of moments, none of wscorr(riley), bscov(prop) and bscov(corr) has been implemented: in fact we don't yet have methods.

Please report any other problems to ian.white@ucl.ac.uk.

References

DerSimonian R, Laird N.
Meta-analysis in clinical trials.
Controlled Clinical Trials 1986; 7: 177-188.
http://www.sciencedirect.com/science/article/pii/0197245686900462

Higgins JPT, Thompson SG.
Quantifying heterogeneity in a meta-analysis. Statistics in Medicine 2002; 21: 1539-58.
http://onlinelibrary.wiley.com/doi/10.1002/sim.1186/abstract

Higgins JPT, Thompson SG, Spiegelhalter DJ.
A re-evaluation of random-effects meta-analysis.
Journal of the Royal Statistical Society (A) 2009; 172: 137-159.
http://onlinelibrary.wiley.com/doi/10.1111/j.1467-985X.2008.00552.x/abstract

Jackson D, White IR, Thompson SG.
Extending DerSimonian and Laird's methodology to perform multivariate random effects meta-analyses.
Statistics in Medicine 2010; 29: 1282-1297.
http://onlinelibrary.wiley.com/doi/10.1002/sim.3602/abstract

Jackson D, Riley R, White IR.
Multivariate meta-analysis: potential and promise.
Statistics in Medicine 2011; 30: 2481-2498.
http://onlinelibrary.wiley.com/doi/10.1002/sim.4172/abstract

Jackson D, White IR, Riley R.
Quantifying the impact of between-study heterogeneity in multivariate meta-analyses
Statistics in Medicine 2012; 31: 3805-3820.
http://onlinelibrary.wiley.com/doi/10.1002/sim.5453/abstract

Jackson D, White IR, Riley R.
A matrix based method of moments for fitting the random effects model
for meta-analysis and meta-regression.
Biometrical Journal 2013; 55: 231-245.
http://onlinelibrary.wiley.com/doi/10.1002/bimj.201200152/abstract

Jackson D, White IR, Price M, Copas J, Riley RD.
Borrowing of strength and study weights in multivariate and network meta-analysis.
Statistical Methods in Medical Research 2017; 26(6), 2853–2868.
https://doi.org/10.1177/0962280215611702

Riley RD.
Multivariate meta-analysis: the effect of ignoring within-study correlation. Journal of the Royal Statistical Society (A) 2009; 172: 789-811.
http://onlinelibrary.wiley.com/doi/10.1111/j.1467-985X.2008.00593.x/abstract

Riley RD, Thompson JR, Abrams KR.
An alternative model for bivariate random-effects meta-analysis when the within-study correlations are unknown. Biostatistics 2008; 9: 172-186.
http://biostatistics.oxfordjournals.org/content/9/1/172.short

Salanti G, Ades A, Ioannidis J.
Graphical methods and numerical summaries for presenting results from multiple-treatment meta-analysis:
an overview and tutorial.
Journal of Clinical Epidemiology 2011; 64: 163-171.
http://www.ncbi.nlm.nih.gov/pubmed/20688472

* White IR.
Multivariate random-effects meta-analysis.
Stata Journal 2009; 9: 40-56.
http://www.stata-journal.com/article.html?article=st0156

* White IR.
Multivariate random-effects meta-regression: Updates to mvmeta.
Stata Journal 2011; 11: 255-270.
http://www.stata-journal.com/article.html?article=st0156_1

White IR, Barrett JK, Jackson D, Higgins JPT.
Consistency and inconsistency in network meta-analysis:
model estimation using multivariate meta-regression.
Research Synthesis Methods 2012; 3: 111-125.
http://onlinelibrary.wiley.com/doi/10.1002/jrsm.1045/abstract

* Please use these references to cite this program.

Author and updates

Ian White, MRC Clinical Trials Unit at UCL, London, UK.
Email ian.white@ucl.ac.uk.

You can get the latest version of this package by visiting https://github.com/UCL/mvmeta
or within Stata by running
net from https://raw.githubusercontent.com/UCL/mvmeta/master/package/.