Network metaanalysis of multiple outcome measures accounting for borrowing of information across outcomes
 Felix A Achana^{1}Email author,
 Nicola J Cooper^{1},
 Sylwia Bujkiewicz^{1},
 Stephanie J Hubbard^{1},
 Denise Kendrick^{2},
 David R Jones^{1} and
 Alex J Sutton^{1}
DOI: 10.1186/147122881492
© Achana et al.; licensee BioMed Central Ltd. 2014
Received: 8 January 2014
Accepted: 2 July 2014
Published: 21 July 2014
Abstract
Background
Network metaanalysis (NMA) enables simultaneous comparison of multiple treatments while preserving randomisation. When summarising evidence to inform an economic evaluation, it is important that the analysis accurately reflects the dependency structure within the data, as correlations between outcomes may have implication for estimating the net benefit associated with treatment. A multivariate NMA offers a framework for evaluating multiple treatments across multiple outcome measures while accounting for the correlation structure between outcomes.
Methods
The standard NMA model is extended to multiple outcome settings in two stages. In the first stage, information is borrowed across outcomes as well across studies through modelling the withinstudy and betweenstudy correlation structure. In the second stage, we make use of the additional assumption that intervention effects are exchangeable between outcomes to predict effect estimates for all outcomes, including effect estimates on outcomes where evidence is either sparse or the treatment had not been considered by any one of the studies included in the analysis. We apply the methods to binary outcome data from a systematic review evaluating the effectiveness of nine home safety interventions on uptake of three poisoning prevention practices (safe storage of medicines, safe storage of other household products, and possession of poison centre control telephone number) in households with children. Analyses are conducted in WinBUGS using Markov Chain Monte Carlo (MCMC) simulations.
Results
Univariate and the first stage multivariate models produced broadly similar point estimates of intervention effects but the uncertainty around the multivariate estimates varied depending on the prior distribution specified for the betweenstudy covariance structure. The second stage multivariate analyses produced more precise effect estimates while enabling intervention effects to be predicted for all outcomes, including intervention effects on outcomes not directly considered by the studies included in the analysis.
Conclusions
Accounting for the dependency between outcomes in a multivariate metaanalysis may or may not improve the precision of effect estimates from a network metaanalysis compared to analysing each outcome separately.
Keywords
Network metaanalysis Mixed treatment comparisons Multiple outcomes Multivariate WinBUGSBackground
Metaanalysis or the quantitative synthesis of evidence, usually from systematic reviews, has become a popular tool in healthcare evaluations [1, 2]. Largely driven by a desire for more realistic synthesis of complex healthcare evidence, increasingly sophisticated methodology has been developed. One area of metaanalysis that has seen significant methodological development is the application of multivariate statistical methods for the comparison of treatments on two or more endpoints (usually known as multivariate metaanalysis) [3–8]. These methods are appealing because many studies and systematic reviews focus on broad health effects and therefore typically report several outcome measures [4, 6, 9]. In such instances, the multivariate approach offers some advantages over separate univariate analyses including the ability to account for the interrelationship between outcomes and borrow strength across studies as well as across outcomes [10] through modelling the correlation structure [7, 11]. This can potentially reduce outcome reporting bias [12] and the uncertainty with which intervention effects are estimated. Additionally, in a decision making context where the synthesis is meant to inform a health economic evaluation, accounting for the correlations between effect estimates on different outcomes is important as the dependence between outcomes may have implication for estimating quality of life or economic consequences associated with treatment [13]. An example is the situation where a particularly effective treatment for a disease condition is associated with a large side effect profile. Ignoring information about the interrelationships between beneficial and ‘side effect’ endpoints in such instances may have implications for quantifying the benefits associated with treatment.
When summarising effectiveness evidence, correlations between the effectiveness estimates typically arise at either withinstudy and or betweenstudy levels. At the withinstudy level, correlations arise mainly due to differences in patientlevel characteristics. They are rarely reported in the published literature and usually have to be estimated from external sources such as individual patient level data if available or elicited from expert opinion [8, 11, 14, 15]. At the betweenstudy level, correlations arise from i) differences in the distribution of patientlevel characteristics across studies, in which case they will be related to the withinstudy correlations and/or ii) differences in the distribution of other studylevel characteristics such as study design, population and baseline disease severity [16]. The withinstudy correlations thus give an indication of the association between multiple endpoints within a study while the betweenstudy correlations indicate how the underlying true studyspecific effects on different outcomes vary jointly across studies.
A second area of rapid methodological development is network metaanalysis (NMA) [17], also known as mixed treatment comparison metaanalysis [18–20] or multiple treatment metaanalysis [21–23]. NMA methods extend standard pairwise metaanalysis to enable simultaneous comparison of multiple treatments while maintaining randomisation of individual studies [18]. The method enables ‘direct’ evidence (i.e. evidence from studies directly comparing two interventions of interest) and ‘indirect’ evidence (i.e. evidence from studies that do not compare the two interventions directly) to be pooled under the assumption of evidence consistency [24]. Estimates of intervention effects can then be obtained, including effects between treatments not directly compared within any one individual study [19]. NMA methods thus provide a coherent framework for appraising all available evidence relevant to a specific decision problem. The results from such analyses are increasingly being used to inform economic evaluations in healthcare decision making where coherent decisions (about judicious use of scarce resource) need to be made based on sound appraisal of all available evidence.
Approaches to extend NMA methodology to multiple outcome settings have been proposed in the literature [13, 25–27], initially focusing on mutually exclusive competing risk outcomes [13] or a single outcome measured at multiple time points [26, 28]. More recently, Efthimiou et al.[14] proposed a method for modelling multiple correlated outcomes in networks of evidence with binary outcome measures. The proposed method accounts for both the withinstudy and betweenstudy correlation structure and includes a strategy for eliciting expert opinion to inform the withinstudy correlations. This paper contributes to the growing literature on the simultaneous evaluation of correlated outcomes. We do this in two stages. In the first stage (labelled as model 2 in the remainder of the paper), information is borrowed across studies as well as across outcomes through modelling the correlations between effectiveness estimates on different outcomes. In the second stage (labelled as model 3 in the remainder of the paper), additional information is borrowed across outcomes based on ideas for combining evidence across human and animal studies originally proposed by DuMouchel and Harris [29] and also revisited by Jones et al.[30]. The proposed second stage analysis methods allows: i) disconnected treatments to be incorporated as nodes in a network of evidence and ii) prediction of intervention effects for outcomes where evidence from primary studies is either sparse or not directly available from any one study included in the analysis. The motivating application area is injury prevention in children where a broad array of outcomes and intervention packages have been evaluated with the aim of increasing safety practices around the home (to ultimately reduce household injuries).
The remainder of this paper is structured as follows: the example dataset is first described followed by a Methods section describing the statistical models developed and implementation of the models. These are followed by sections presenting the results of applying the methods to the motivating dataset and a discussion.
Dataset
 a)
Safe storage of medicines
 b)
Safe storage of other household products (e.g. cleaning products) and
 c)
Possession of a poison control centre (PCC) telephone number.
Summary of the available evidence
Outcome information (no. of events/no. of households in control versus (vs.) treatment arm)  

Comparison  First author and year of publication  IPD  Safe storage of medicines  Safe storage of other household products  Possession of a PCC number 
Usual care (1) vs. Education (2)  Gielen 2007  Yes  178/271 vs. 188/249^{Ɨ}  44/62 vs. 57/73^{Ɨ}  
Nansel 2002  Yes  83/89 vs. 79/85  65/89 vs. 66/85  59/89 vs. 63/85  
Nansel 2008  Yes  72/74 vs. 140/144^{†}  59/73 vs. 117/144^{†}  50/59 vs. 90/119^{†}  
Kelly B 1987  No  54/54 vs. 55/55  43/54 vs. 49/55  
McDonald 2005  No  4/57 vs. 6/60  3/57 vs. 6/61  
Kelly N 2003  No  45.56/136.68 vs. 112.95/137.63^{*}  
Usual care (1) vs. Education + free/low cost safety equipment (3)  Clamp 1998  Yes  68/82 vs. 79/83  49/82 vs. 59/83  
Woolf 1987  No  29/143 vs. 47/119  
Woolf 1992  No  60/151 vs. 89/150  59/151 vs. 117/150  
Usual care (1) vs. Education + equipment (3) vs. Education + equipment + home safety inspection (4)  Babul 2007  Yes  147/149 vs. 171/173 vs. 160/163  
Usual care (1) vs. Education + equipment + home safety inspection (4)  Hendrickson 2002  Yes  14/40 vs. 34/38  8/40 vs. 34/38  
Swart 2008  No  70.26/79.58 vs. 74.07/80^{*}  46.86/57.96 vs. 50.87/58.27^{*}  
Kendrick 1999  Yes  317/367 vs. 322/363  
Usual care (1) vs. Education + equipment + fitting (5)  Watson 2005  Yes  683/738 vs. 712/762  327/669 vs. 368/693  
Usual care (1) vs. Education + home safety inspection (6)  Petridou 1997  No  67.26/100.12 vs. 71.08/97.83^{*}  
Usual care (1) vs. Education + equipment + home safety inspection + fitting (7)  Schwarz D 1993  No  88.42/248.37 vs. 128.16/248.37^{*}  
Phelan 2011  No  16/138 vs. 71/139  
Usual care (1) vs. Home visit (8)  Johnson 2006  No  82/91 vs. 222/232^{†}  
Education (2) vs. Education + equipment (3)  Posner 2004  Yes  14/47 vs. 19/49  22/47 vs. 34/49  27/47 vs. 35/49 
Education (2) vs. Education + equipment + fitting (5)  Sznajder 2003  Yes  44/49 vs. 43/45  32/41 vs. 40/48  
Education + equipment + home safety inspection (4) vs. Education + equipment + home safety inspection + fitting (7)  King J 2001  No  261/469 vs. 273/482  
Education + equipment (3) vs. Equipment (9)  Dershewitz 1979  No  22/101 vs. 20/104  1/101 vs. 0/104 
Thirteen of the 22 studies considered at least two of the three outcomes. Of these, 8 considered storage of medicines and storage of other household products, 2 considered storage of other household products and possession of a PCC telephone number, and 3 considered all three outcome measures. Individual patient data (IPD) were available for 8 of the 13 studies, of which 7 were in a format suitable for the analysis reported here as explained by the footnotes in Table 1.
 1)
Usual care (UC)
 2)
Education (E)
 3)
Education + provision of free/low cost equipment (E + FE)
 4)
Education + provision of free/low cost equipment + home safety inspection (E + FE + HSI)
 5)
Education + provision of free/low cost equipment + fitting of equipment (E + FE + F)
 6)
Education + home safety inspection (E + HSI)
 7)
Education + provision of free/low cost equipment + home safety inspection + fitting of equipment (E + FE + HSI + F)
 8)
Education + home visit (E + HV)
 9)
Provision of free/low cost equipment (FE).
Methods
In this section, we first present the NMA statistical model for one binary outcome measure and then extend it to compare multiple interventions across multiple outcomes. Throughout the paper, we refer to these single and multiple outcome models as univariate and multivariate NMAs respectively. Where studies report multiple outcomes, these will not be independent as each household provides information on the different outcome measures within intervention arms. The multivariate model takes this correlation structure into account by allowing the intervention effects measured by one outcome to be correlated with the intervention effects measured by other outcomes.
Model 1: Univariate NMA
The analysis is conducted within a Bayesian framework requiring prior distributions to be specified for all model parameters. Accordingly, we specified minimally informative prior distributions corresponding to a Normal (0,10^{3}) prior distribution for the pooled mean effects relative to usual care, d _{(Ak)} and the studyspecific baseline effects, μ _{ ib } and a Uniform (0,2) prior distribution for the betweenstudy standard deviation on the log odds ratio scale σ[34].
Model 2: Multivariate NMA
We extend the univariate NMA model defined above to the multiple outcomes settings in order to account for correlations between intervention effects on different outcomes. The method presented here follows from Ades et al. (2010) NMA with competing risks model [13] where only the withinstudy correlations are taken into account. We extend their method to account for the betweenstudy correlation as well.
where (μ _{ ib1}, μ _{ ib2}, ⋯, μ _{ ibM }) and (δ _{ i(bk)1}, δ _{ i(bk)2}, ⋯, δ _{ i(bk)M }) represent vectors of ‘true’ baseline and studyspecific effects in study i with baseline treatment b respectively. The quantities $\left({\widehat{\theta}}_{\mathit{ik}1},\phantom{\rule{0.25em}{0ex}}{\widehat{\theta}}_{\mathit{ik}2},\cdots ,\phantom{\rule{0.5em}{0ex}}{\theta}_{\mathit{ikM}}\right)$ and (θ _{ ik1}, θ _{ ik2}, ⋯, θ _{ ikM }) represent vectors of observed and ‘true’ logodds of response in arm k of study i and S _{ ik } is the associated withinstudy covariance matrix usually assumed known but estimated in practice from the data [35]. We calculated $\left({\widehat{\theta}}_{\mathit{ik}1},\phantom{\rule{0.5em}{0ex}}{\widehat{\theta}}_{\mathit{ik}2}\phantom{\rule{0.25em}{0ex}},\cdots ,\phantom{\rule{0.25em}{0ex}}{\widehat{\theta}}_{\mathit{ikM}}\right)$ and the diagonal elements of S _{ ik } using standard formulae for logodds and variance of the logodds [2]. We applied a continuity correction by adding 0.5 to the numerators and 1 to the denominators of studies with 0% or 100% event rate in one of the treatment arms [36, 37]. The offdiagonals of S _{ ik } were calculated from estimates of withinstudy correlations ${r}_{\mathit{ik}}^{\mathit{mn}}$ between outcomes m and n (m ≠ n) in arm k of study i obtained from studies with IPD (see Additional file 1: Table S1). The method used to estimate the correlations from the IPD is described in the implementation section below.
When summarising evidence across multiple endpoints, it is common to encounter instances where some studies do not report information for all outcomes of interest leading to incomplete vectors with missing studyspecific effects for the outcomes not reported [5, 10]. Such studies can be included in our model under the assumption that the effects for outcomes not reported are missing at random. When implemented using the WinBUGS software, the missing study effects and standard errors are coded as NA in the data, a strategy previously outlined in Bujkiewicz et al.[10] and Dakin et al.[28]. This enables WinBUGS to automatically ‘impute’ values for the missing information under missing at random assumption with predicted distributions.
where, as in the univariate case, σ _{ m } represent the common betweenstudy standard deviation or heterogeneity parameter specific to outcome m.
where the offdiagonal elements of R contain correlation terms and diagonal elements equal 1. Lu and Ades [43] and also Wei and Higgins [39] showed that R can be written as R = L ^{ T } L using Cholesky decomposition where L is an upper triangular matrix. The spherical parameterization technique [39, 43] can be used to express R in terms of sine and cosine functions of the elements in L. Using this later technique, we specified Uniform (0,π) prior distributions for the spherical coordinate φ _{ mn } in our model to ensure that elements of the correlation matrix R lie in the interval (−1,1). Finally, the elements of V ^{1/2} correspond to the betweenstudy standard deviation terms in Σ _{(M × M)} and are given independent Uniform (0,2) prior distributions as in the univariate case (model 1).
Model 3: Borrowing strength across interventions and outcomes
Equation (10) thus embodies an assumption of equal or constant relative potency of treatments across outcomes which imply exchangeability of the relative effects between the nonreference/baseline treatments indicated by equation (11). For our example dataset, this implies that missing intervention effects for comparisons with the usual care intervention can be predicted directly from equation (10) as a linear combination of γ _{ m } and α _{ k } assuming that each treatment effect relative to usual care is reported on at least one outcome. The missing intervention effects between non reference/baseline treatments if required can similarly be predicted directly from the model as linear combinations of the intervention effects relative to usual care.
The parameter τ controls the accuracy of the constant relative potency assumption. Values of τ close to zero would thus indicate a high degree of confidence (and support from the data) in the parallelism of effect profiles across outcomes and the constant relative potency assumption. Conversely, larger values of τ would indicate otherwise.
 i.
Normal prior distribution centred on 0 with large variance and constrained to be positive, τ ~ N(0, 10^{2}), τ ≥ 0
 ii.
Gamma prior distribution placed on the precision: τ ^{− 2} ~ Gamma(0.001, 0.001).
The results of the sensitivity analyses are presented in Additional file 1: Figure S1.
There is a limitation to the number of data (i.e. intervention effects relative to the usual care) on outcomes allowed to be missing for the model hyperparameters to be identifiable. For K interventions and M outcomes, there will be (K − 1) × M equation (10) that are used to estimate a total of (K − 1) + M hyperparameters (i.e. (K − 1) of α _{ k } and M of γ _{ m } hyperparameters). Therefore no more than ((K − 1) × M) − ((K − 1) + M) missing values in total are allowed. For example, for K = 3 treatments and M = 2 outcomes, data has to be available on both outcomes for both treatment comparisons with the baseline when the prior distributions are noninformative. When large number of data on outcomes is missing, placing informative prior distributions on the hyperparameters can improve convergence.
Implementation of the models
where ${r}_{\mathit{ik}}^{\mathit{mn}}$ is the withinstudy correlation between the outcomes m and n effects measured on the logodds scale in arm k of study i, and $\overline{{r}^{\mathit{mn}}}$ and var (r ^{ mn })are the mean and variance of the withinstudy correlation between outcomes m and n effects measured by IPD respectively.
We also assessed consistency of the evidence within each network using a method based on node splitting [24]. We found no evidence of conflict between the direct and indirect sources (Additional file 1: Table S2) in all three outcome networks.
We fitted all models described above in WinBUGS [47] using Markov Chain Monte Carlo (MCMC) simulations. The univariate models were fitted separately for each outcome using WinBUGS code available from Dias et al.[48]. The WinBUGS code for the multivariate models is provided in the appendices 1 and 2 of the Additional file 1. Convergence was assessed by examination of the trace and autocorrelation plots and the RubinGelman statistic after running 400 000 simulations and discarding the first 200 000 samples as ‘burn in samples’.
Results
Univariate and multivariate analyses
Posterior median and 95% credible intervals of the betweenstudy standard deviation and correlation parameters
Parameter  Description/prior distribution  Model 1: univariate  Model 2a: Multivariate using inverseWishart prior distribution for Σ_{(M × M)}  Model 2b: Multivariate using a separation strategy to specify priors for elements of Σ_{(M × M)}  Model 3: Multivariate with extrapolation of effects across outcomes 

σ _{1}  Betweenstudy standard deviation: safe storage of medicines  0.26 (0.03, 1.02)  0.58 (0.33, 1.18)  0.27 (0.01, 1.08)  0.23 (0.01, 0.80) 
σ _{2}  Betweenstudy standard deviation: safe storage of other household products  0.56 (0.13, 1.27)  0.62 (0.35, 1.15)  0.47 (0.04, 1.18)  0.31 (0.01, 0.81) 
σ _{3}  Betweenstudy standard deviation: PCC  1.16 (0.57, 1.93)  0.94 (0.53, 1.99)  1.18 (0.57, 1.93)  1.08 (0.58, 1.85) 
τ  Primary analysis: τ ~Uniform (0, 2)        0.10 (0.01, 0.53) 
τ  Sensitivity analysis: τ ~ Normal (0, 10^{2}), τ ≥ 0        0.11 (0.00, 0.56) 
τ  Sensitivity analysis: τ ^{2} ~ Inverse − Gamma (0.001, 0.001)        0.08 (0.02, 0.36) 
ρ ^{12}  Betweenstudy correlation [medicines, other household products]    0.03 (−0.73, 0.76)  0.05 (−1.00, 1.00)  0.45 (−0.99, 1.00) 
ρ ^{13}  Betweenstudy correlation [medicines, PCC]    0.06 (−0.80, 0.81)  0.20 (−1.00, 1.00)  0.50 (−0.98, 1.00) 
ρ ^{23}  Betweenstudy correlation [Other household products, PCC]    0.08 (−0.81, 0.83)  0.13 (−0.97, 0.98)  0.60 (−0.87, 0.99) 
Borrowing strength across outcomes
It can be seen from Figure 2 that the effect of E + HSI and E + HV relative to usual care intervention on safe storage of medicines and safe storage of other household products, and E + FE + F and FE on possession of a PCC telephone number were only estimated in model 3 as none of the studies had trialled these interventions on the respective outcomes. In this model, estimates of relative effects between non reference/baseline treatments were assumed to be exchangeable across outcomes, which enabled estimates to be obtained for all outcomes by predicting effects where the interventions have not been considered for the particular outcome of interest. For the intervention/outcome pair where data from trials were available, this extrapolation step had the additional effect of producing more precise estimates of the treatment effect in comparison to the models that do not assume exchangeability effects across outcomes.
The posterior median and 95% credible intervals of intervention effects relative to usual care were unaffected by placing alternative minimally informative prior distributions on τ (Additional file 1: Appendix 2). The posterior median and credible intervals for τ (Table 2) were similarly not sensitive to the choice of prior distribution placed on τ in the primary and sensitivity analyses. The posterior median estimates of τ were all close to zero, which suggest that assumptions about the parallelism of effect profiles across outcomes is supported by the data. The estimates of τ would thus suggest with 95% probability, that based on the information in our example dataset, the extrapolation model could be accurate to within a factor of about e ^{(2 × 0.10)} = 1.24 (95% CrI 1.01 to 2.87).
Discussion
We have presented methods for simultaneous comparison of multiple treatments across multiple outcome measures while preserving the internal randomisation of individual studies. Our method may be viewed as an extension of Ades et al.’s (2010) NMA with competing risks paper [13] wherein only the withinstudy correlation is taken into account. We have extended their method to account for the dependency between outcome effects across studies as well as withinstudies.
In this particular application of the multivariate approach, accounting for the correlation between outcomes alone (models 2a and 2b) did not reduce the uncertainty around estimates of intervention effects compared to analysing each outcome separately (model 1). Assuming that intervention effects are exchangeable across outcome did however lead to a modest reduction in uncertainty around effectiveness estimates (model 3).
The betweenstudy correlations were estimated with considerable uncertainty (Table 2) and appear to have little impact on overall effect estimates. This may be because the betweenstudy correlation arises due to, among other things, differences in studylevel characteristics that also give rise to betweenstudy heterogeneity in a metaanalysis. Based on a criterion outlined in Spiegelhalter et al.[49] the posterior median estimates of the betweenstudy standard deviations, σ _{1} and σ _{2} on the log odds ratio scale (Table 2) could be interpreted as indicating evidence of low to moderate heterogeneity for storage of medicines and storage of other household products outcomes. Only the estimates for possession of poison control centre number exhibited a considerable degree of heterogeneity. Consequently, the posterior medians of the betweenstudy correlations were small. There was therefore very little gain (in terms of increasing the precision of estimates) from formulating the betweenstudy covariance structure described for the analysis presented here. Accounting for the betweenstudy correlation is likely to be beneficial in situations where the betweenstudy variance (heterogeneity) is large relative to withinstudy variances.
We opted to incorporate the withinstudy correlation through the armspecific effects (logodds) rather than the studyspecific treatment difference (logodds ratio) as is often done in multivariate metaanalysis [3, 4, 15, 38]. This approach greatly simplifies the likelihood for multiarm studies because treatment arms can be considered independent as a consequence of randomisation. Hence, there is no requirement to account for the additional correlations between effect estimates which share a common comparator treatment in the model likelihood [50]. The armbased approach is also likely to be useful when (as is typical with many practical application of multivariate metaanalysis) the withinstudy correlations are not available [10, 12, 15, 51] and have to be obtained from an external source such as expert opinion [14]. In such situations, formulating questions about correlations between outcomespecific event probabilities (which can be used directly in an armbased approach) is more likely to be intuitive and easily understood by nonstatistician healthcare experts than questions about correlations between intervention effects. It is acknowledged however, that the correlations between the intervention effects if required can easily be obtained from the correlations between the outcomes [14, 51].
At the betweenstudy level, we assumed a common correlation structure across treatments in addition to the common variance assumption underlying most practical application of NMA methods. The common correlation assumption implies that if several separate multivariate metaanalyses were conducted with the same outcomes, each with a different set of k versus b comparison, the assumption is that the betweenstudy correlations would be the same across the different sets of bk comparisons. We suggested this structure to simplify the covariance structure and reduce the number of parameters in the model. Appropriateness of such modelling assumptions would need to be considered carefully and assessed when it is feasible to do so.
Initially we specified an inverseWishart distribution for the betweenstudy covariance matrix Σ _{(M × M)}. However, we believe this prior distribution to be influential due to the small number of studies in our example dataset relative to the number of outcomes. Under these conditions, the inverseWishart prior distribution produced upwardlybiased estimates of σ _{1} and σ _{2} and downward bias in the estimate for σ _{3} when compared to the corresponding estimates obtained from the univariate model (Table 2). These findings are consistent with observations in the univariate case where the use of a Gamma prior distribution (which is the univariate analogue of the InverseWishart prior distribution) can lead to an overestimation of the heterogeneity parameter when the true value is close to 0 [46, 52]. As an alternative to an inverseWishart prior distribution therefore, we followed the spherical decomposition technique suggested by Lu and Ades [43]. This parameterization offered greater flexibility in formulating independent prior distributions for the standard deviation and correlation terms in Σ _{(M × M)}.
An obvious limitation to implementation of the multivariate models presented in this paper is the limited availability of data including i) the problem of missing withinstudy correlations and ii) the requirement for a relatively large number of studies to estimate all model parameters. The problem of missing withinstudy correlations has traditionally hampered the widespread application of multivariate metaanalysis [7, 10, 15]. In our example, IPD was available from a proportion of the included studies and we have used correlations estimated from the IPD to formulate informative prior distributions for the withinstudy model. Alternative approaches to dealing with missing withinstudy correlations when IPD is not available include: i) using the observed correlation from the summary studyspecific effects [12], ii) eliciting information about the correlations from external sources such as clinical experts [14] and iii) specifying ‘vague’ prior distributions for analysis conducted within a Bayesian framework [6].
The second data issue concerns the number of studies needed to estimate the full unstructured betweenstudy covariance matrix presented in equation (6). We anticipate a large number of multiarm studies reporting across the three outcomes will be needed to identify Σ _{(bk)} and estimate all model parameters. This can be problematic considering the fact that most applications of network metaanalysis typically include mostly twoarm studies with very small numbers of multiarm studies. Even with the simplification of the betweenstudy covariance matrix given in equation (5), a relatively large number of studies in comparison to the total number of outcomes being considered may still be needed. We are unable to answer the question of how many studies should be considered large enough for a NMA with multiple outcomes. As a guide, Wei and Higgins [39] recently estimated 15, 27 and 42 studies as a minimum for multivariate pairwise metaanalysis with two, three and fouroutcomes respectively. Hence, we believe an even larger number of studies will be required for the NMA with multiple outcomes.
Another limitation of the multivariate models presented here is that they rely on the normal approximation to binomial distribution to incorporate the withinstudy correlations in the model. The normal approximation frequently fails and may not provide adequate fit to the data in the presence of studies with zero or a small number of events, necessitating use of continuity corrections. We were unable to use the exact binomial distribution as our primary interest was to develop models for summary binary data where outcomes are not mutually exclusive, and where it is not reasonable to assume that withinstudy correlations are zero so that the likelihood factorises easily as in Arends et al.[3]. Further methodological investigations into modelling multivariate summary data that is not normally distributed will therefore be useful. An example is provided in Chu et al.[53] where parameterization of the withinstudy model enabled the special case of diagnostic sensitivity and specificity to be jointly modelled with disease prevalence using a trivariate binomial likelihood. In the interim, an alternative formulation which bypasses the need for approximating normal distributions is to directly model the IPD where this is available. This will require extending Saramago et al.’s[54] NMA model with aggregate and individual participant level data from single outcome to multiple outcome settings.
We assessed the consistency of each outcome network separately using the method of node splitting [24]. We found no evidence of conflict between the direct and indirect sources on pairwise contrasts that have both sources of evidence in model 1. We did not assess the consistency of the multivariate estimates partly because we are unaware of current methods for carrying out this type of assessment. We are investigating extensions of the nodesplit method to multiple outcome networks and investigate the effect of jointly synthesising evidence across multiple endpoints on evidence consistency in a simulation study.
Our initial motivation for a multiple outcome NMA was to estimate intervention effects for all the outcomes, including effects of interventions on outcomes not considered by any of the studies included in the analysis. This requires the correlation structure between effects on multiple outcomes to be appropriately modelled and also ensuring the mechanism of “borrow strength” across outcomes through the assumption of exchangeability of the random effect across outcomes. This implies a priori assumption that outcomes are related but different and that there is no way of knowing the order of magnitude of effects on outcomes. If this assumption does not hold, it may potentially lead to worse or more biased effectiveness estimates. In our example, the outcomes are similar and measured on the same scale. It would be clearly inappropriate to assume that intervention effects are exchangeable across outcomes that are different in some important respects such as being measured on different scales (e.g. where one outcome reports a weighted mean difference and another outcome reports a logodds ratio) as such estimates will differ in terms of the precision with which they are estimated.
Conclusion
Our aim in this paper was to present methods for simultaneous comparison of multiple treatments across multiple outcome measures while preserving the internal randomisation of individual studies. Application of the method to the poison prevention data yielded similar point estimates of treatment effect to those obtained from a univariate NMA but the uncertainty around the multivariate estimates increased or decreased depending on the prior distribution specified for the betweenstudy covariance structure. The proposed method followed the usual hierarchical approach to multivariate metaanalysis where correlations between outcomes are modelled at the withinstudy and or betweenstudy levels.
Betweenstudy covariance for multiarm studies reporting multiple outcomes
Where δ _{ i(Ak)m } and d _{(Ak)m } are studyspecific and mean effect of treatment k relative A (reference treatment) on outcome m in study i respectively and Σ _{ FULL } is the full (K1) × (K1) blocks of M × M withintreatment betweenoutcome covariance matrix. The parameters in Σ _{ FULL } have the following interpretation:(${\sigma}_{\left(\mathit{Ak}\right)m}^{2}$) indicate the variance of the effect of treatment k (k = B, C, ⋯, K) relative to A on outcome m across studies.
${\rho}_{\left(\mathit{Ak},\mathit{Ak}\right)}^{\mathit{mn}}$indicate the correlation between δ _{ i(Ak)m } and δ _{ i(Ak)n } (i.e. the correlation between the effect of treatment k relative to A on outcome m and the effect of treatment k relative to A on outcome n (m ≠ n)) specific to the Ak comparison.
Abbreviations
 CrI:

Credible interval
 IPD:

Individual patient data
 MCMC:

Markov Chain Monte Carlo
 NMA:

Network metaanalysis
 PCC:

Poison centre control number.
Declarations
Acknowledgements
FAA is funded by the National Institute for Health Research (NIHR) under its Programme Grants for Applied Research programme (Reference Number RPPG 0407–10231) to pursue a PhD in Public Health Evaluation of Interventions.
SB is supported by the MRC Methodology Research Programme (grant number MR/L009854/1).
Authors’ Affiliations
References
 Cooper NJ, Sutton AJ, Ades AE, Paisley S, Jones DR, W. G. O. U. o. E. i. E. D: Models, Use of evidence in economic decision models: Practical issues and methodological challenges. Health Economics. 2007, 16 (12): 12771286.View ArticlePubMedGoogle Scholar
 Sutton A, Abrams K, Jones DR, Sheldon TA, Song F: Methods for Metaanalysis in Medical Research. 2000, Chichester: WileyGoogle Scholar
 Arends LR, Voko Z, Stijnen T: Combining multiple outcome measures in a metaanalysis: an application. Statistics in Medicine. 2003, 22 (8): 13351353.View ArticlePubMedGoogle Scholar
 Berkey CS, Hoaglin DC, AntczakBouckoms A, Mosteller F, Colditz GA: Metaanalysis of multiple outcomes by regression with random effects. Statistics in Medicine. 1998, 17 (22): 25372550.View ArticlePubMedGoogle Scholar
 Jackson D, Riley R, White IR: Multivariate metaanalysis: Potential and promise. Statistics in Medicine. 2011, 30 (20): 24812498.View ArticlePubMedPubMed CentralGoogle Scholar
 Nam IS, Mengersen K, Garthwaite P: Multivariate metaanalysis. Statistics in Medicine. 2003, 22 (14): 23092333.View ArticlePubMedGoogle Scholar
 Riley RD, Abrams KR, Sutton AJ, Lambert PC, Thompson JR: Bivariate randomeffects metaanalysis and the estimation of betweenstudy correlation. BMC Medical Research Methodology. 2007, 7:Google Scholar
 Riley RD, Thompson JR, Abrams KR: An alternative model for bivariate randomeffects metaanalysis when the withinstudy correlations are unknown. Biostatistics. 2008, 9 (1): 172186.View ArticlePubMedGoogle Scholar
 Kendrick D, Coupland C, Mulvaney C, Simpson J, Smith SJ, Sutton A, Watson M, Woods A: Home safety education and provision of safety equipment for injury prevention. Cochrane Database Systematic reviews. 2007, 1: CD005014Google Scholar
 Bujkiewicz S, Thompson JR, Sutton AJ, Cooper NJ, Harrison MJ, Symmons DPM, Abrams KR: Multivariate metaanalysis of mixed outcomes: a Bayesian approach. Statistics in Medicine. 2013, 32 (22): 39263943.View ArticlePubMedPubMed CentralGoogle Scholar
 Riley RD, Abrams KR, Lambert PC, Sutton AJ, Thompson JR: An evaluation of bivariate randomeffects metaanalysis for the joint synthesis of two correlated outcomes. Statistics in Medicine. 2007, 26 (1): 7897.View ArticlePubMedGoogle Scholar
 Kirkham JJ, Riley RD, Williamson PR: A multivariate metaanalysis approach for reducing the impact of outcome reporting bias in systematic reviews. Statistics in Medicine. 2012, 31 (20): 21792195.View ArticlePubMedGoogle Scholar
 Ades AE, Mavranezouli I, Dias S, Welton NJ, Whittington C, Kendall T: Network metaanalysis with competing risk outcomes. Value Health. 2010, 13 (8): 976983.View ArticlePubMedGoogle Scholar
 Efthimiou O, Mavridis D, Cipriani A, Leucht S, Bagos P, Salanti G: An approach for modelling multiple correlated outcomes in a network of interventions using odds ratios. Statistics in Medicine. 2014, 33 (13): 22752287.View ArticlePubMedGoogle Scholar
 Riley RD: Multivariate metaanalysis: the effect of ignoring withinstudy correlation. Journal of the Royal Statistical Society Stat Soc. 2009, 172: 789811.View ArticleGoogle Scholar
 Rodgers M, Sowden A, Petticrew M, Arai L, Roberts H, Britten N, Popay J: Testing Methodological Guidance on the Conduct of Narrative Synthesis in Systematic Reviews. Evaluation. 2009, 15 (1): 4973.View ArticleGoogle Scholar
 Lumley T: Network metaanalysis for indirect treatment comparisons. Statistics in Medicine. 2002, 21 (16): 23132324.View ArticlePubMedGoogle Scholar
 Caldwell DM, Ades AE, Higgins JP: Simultaneous comparison of multiple treatments: combining direct and indirect evidence. British Medical Journal. 2005, 331 (7521): 897900.View ArticlePubMedPubMed CentralGoogle Scholar
 Caldwell DM, Welton NJ, Ades AE: Mixed treatment comparison analysis provides internally coherent treatment effect estimates based on overviews of reviews and can reveal inconsistency. Journal of Clinical Epidemiology. 2010, 63 (8): 875882.View ArticlePubMedGoogle Scholar
 Lu G, Ades AE: Combination of direct and indirect evidence in mixed treatment comparisons. Statistics in Medicine. 2004, 23 (20): 31053124.View ArticlePubMedGoogle Scholar
 Salanti G, Ades AE, Ioannidis JP: Graphical methods and numerical summaries for presenting results from multipletreatment metaanalysis: an overview and tutorial. Journal of Clinical Epidemiology. 2010, 64 (2): 163171.View ArticlePubMedGoogle Scholar
 Salanti G, Higgins JP, Ades AE, Ioannidis JP: Evaluation of networks of randomized trials. Statistical Methods in Medical Research. 2008, 17 (3): 279301.View ArticlePubMedGoogle Scholar
 Salanti G, Marinho V, Higgins JP: A case study of multipletreatments metaanalysis demonstrates that covariates should be considered. Journal of Clinical Epidemiology. 2009, 62 (8): 857864.View ArticlePubMedGoogle Scholar
 Dias S, Welton NJ, Caldwell DM, Ades AE: Checking consistency in mixed treatment comparison metaanalysis. Statistics in Medicine. 2010, 29 (7–8): 932944.View ArticlePubMedGoogle Scholar
 Hong H, Carlin BP, Chu H, Shamliyan TA, Wang S, Kane RL: Methods Research Report: A Bayesian Missing Data Framework for Multiple Continuous Outcome Mixed Treatment Comparisons. 2013, Avaialable online from http://www.ncbi.nlm.nih.gov/books/NBK116689/pdf/TOC.pdf Google Scholar
 Lu G, Ades AE, Sutton AJ, Cooper NJ, Briggs AH, Caldwell DM: Metaanalysis of mixed treatment comparisons at multiple followup times. Statistics in Medicine. 2007, 26 (20): 36813699.View ArticlePubMedGoogle Scholar
 Welton NJ, Cooper NJ, Ades AE, Lu G, Sutton AJ: Mixed treatment comparison with multiple outcomes reported inconsistently across trials: evaluation of antivirals for treatment of influenza A and B. Statistics in Medicine. 2008, 27 (27): 56205639.View ArticlePubMedGoogle Scholar
 Dakin HA, Welton NJ, Ades AE, Collins S, Orme M, Kelly S: Mixed treatment comparison of repeated measurements of a continuous endpoint: an example using topical treatments for primary openangle glaucoma and ocular hypertension. Statistics in Medicine. 2011, 30 (20): 25112535.View ArticlePubMedGoogle Scholar
 DuMouchel WH, Harris JE: Bayes Methods for Combining the Results of Cancer Studies in Humans and Other Species. Journal of the American Statistical Association. 1983, 78 (382): 293308.View ArticleGoogle Scholar
 Jones DR, Peters JL, Rushton L, Sutton AJ, Abrams KR: Interspecies extrapolation in environmental exposure standard setting: A Bayesian synthesis approach. Regulatory Toxicology Pharmacology. 2009, 53 (3): 217225.View ArticlePubMedGoogle Scholar
 Kendrick D, Young B, MasonJones AJ, Ilyas N, Achana FA, Cooper NJ, Hubbard SJ, Sutton AJ, Smith S, Wynn P, Mulvaney C, Watson MC, Coupland C: Home safety education and provision of safety equipment for injury prevention. Cochrane Database Systematic reviews. 2012, 9Google Scholar
 Johnston BD, Huebner CE, Anderson ML, Tyll LT, Thompson RS: Healthy steps in an integrated delivery system  Child and parent outcomes at 30 months. Archives of Pediatrics and Adolescent Medicine. 2006, 160 (8): 793800.View ArticlePubMedGoogle Scholar
 Babul S, Olsen L, Janssen P, McIntee P, Raina P: A randomized trial to assess the effectiveness of an infant home safety programme. International Journal of Control and Safety Promotion. 2007, 14 (2): 109117.View ArticleGoogle Scholar
 Dias S, Sutton AJ, Ades AE, Welton NJ: A Generalized Linear Modeling Framework for Pairwise and Network Metaanalysis of Randomized Controlled Trials. Medical Decision Making. 2012, 33 (5): 60717.View ArticlePubMedGoogle Scholar
 van Houwelingen HC, Arends LR, Stijnen T: Advanced methods in metaanalysis: multivariate approach and metaregression. Statistics in Medicine. 2002, 21 (4): 589624.View ArticlePubMedGoogle Scholar
 Dershewitz RA, Williamson JW: Prevention of Childhood Household Injuries  Controlled ClinicalTrial. American Journal of Public Health. 1977, 67 (12): 11481153.View ArticlePubMedPubMed CentralGoogle Scholar
 Kelly B, Sein C, McCarthy PL: Safety education in a pediatric primary care setting. Pediatrics. 1987, 79 (5): 818824.PubMedGoogle Scholar
 Mavridis D, Salanti G: A practical introduction to multivariate metaanalysis. Statistical Methods in Medical Research. 2013, 22 (2): 133158.View ArticlePubMedGoogle Scholar
 Wei Y, Higgins JPT: Bayesian multivariate metaanalysis with multiple outcomes. Statistics in Medicine. 2013, 32 (17): 29112934.View ArticlePubMedGoogle Scholar
 Carter P, Achana F, Troughton J, Gray LJ, Khunti K, Davies MJ: A Mediterranean diet improves HbA1c but not fasting blood glucose compared to alternative dietary strategies: a network metaanalysis. Journal of Human Nutrition and Dietetics. 2014, 27 (3): 280297.View ArticlePubMedGoogle Scholar
 Cooper NJ, Sutton AJ, Lu G, Khunti K: Mixed comparison of stroke prevention treatments in individuals with nonrheumatic atrial fibrillation. Archives Internal Medicine. 2006, 166 (12): 12691275.View ArticleGoogle Scholar
 Spiegelhalter D, Thomas A, Best N, Lunn D: WinBUGS User Manual. Available from http://www.politicalbubbles.org/bayes_beach/manual14.pdf, 2003. Accessed 10^{th} December 2013
 Lu GB, Ades A: Modeling betweentrial variance structure in mixed treatment comparisons. Biostatistics. 2009, 10 (4): 792805.View ArticlePubMedGoogle Scholar
 Barnard J, McCulloch R, Meng XL: Modeling covariance matrices in terms of standard deviations and correlations, with application to shrinkage. Statistica Sinica. 2000, 10 (4): 12811311.Google Scholar
 DuMouchel W, Groer PG: A Bayesian Methodology for Scaling Radiation Studies from Animals to Man. Health Physics. 1989, 57: 411418.View ArticlePubMedGoogle Scholar
 Lambert PC, Sutton AJ, Burton PR, Abrams KR, Jones DR: How vague is vague? A simulation study of the impact of the use of vague prior distributions in MCMC using WinBUGS. Statistics in Medicine. 2005, 24 (15): 24012428.View ArticlePubMedGoogle Scholar
 Lunn DJ, Thomas A, Best N, Spiegelhalter D: WinBUGS  A Bayesian modelling framework: Concepts, structure, and extensibility. Statistical Computing. 2000, 10 (4): 325337.View ArticleGoogle Scholar
 Dias S, Welton NJ, Sutton AJ, Ades AE: NICE DSU Technical Support Document 2: A Generalised Linear Modelling Framework for Pairwise and Network Metaanalysis of Randomised Controlled Trials. 2011, Available from http://www.nicedsu.org.uk Google Scholar
 Spiegelhalter DJ, Abrams KR, Myles JP: Bayesian Approaches to Clinical Trials and HealthCare Evaluation. Statistics in Practice. 2004, Chichester: John Wiley & SonsGoogle Scholar
 Franchini AJ, Dias S, Ades AE, Jansen JP, Welton NJ: Accounting for correlation in network metaanalysis with multiarm trials. Research Synthesis Methodology. 2012, 3 (2): 142160.View ArticleGoogle Scholar
 Wei Y, Higgins JPT: Estimating withinstudy covariances in multivariate metaanalysis with multiple outcomes. Statistics in Medicine. 2013, 32 (7): 11911205.View ArticlePubMedGoogle Scholar
 Gelman A: Prior distributions for variance parameters in hierarchical models(Comment on an Article by Browne and Draper). Bayesian Analysis. 2006, 1 (3): 515533.Google Scholar
 Chu HT, Nie L, Cole SR, Poole C: Metaanalysis of diagnostic accuracy studies accounting for disease prevalence: Alternative parameterizations and model selection. Statistics in Medicine. 2009, 28 (18): 23842399.View ArticlePubMedGoogle Scholar
 Saramago P, Sutton AJ, Cooper NJ, Manca A: Mixed treatment comparisons using aggregate and individual participant level data. Statistics in Medicine. 2012, 31 (28): 35163536.View ArticlePubMedGoogle Scholar
 The prepublication history for this paper can be accessed here:http://www.biomedcentral.com/14712288/14/92/prepub
Prepublication history
Copyright
This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly credited. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.