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

Given arm-level binary data of the form presented in Table

1, a random effects NMA may be specified using the method of Lu and Ades [

20]. It is assumed that the occurrence of

*r*
_{
ik
} events out of a total of

*n*
_{
ik
} households in the

*k*th-arm (

*k = A,B,C,….,*) of the

*i*th-study follow a binomial distribution with underlying event probability

*p*
_{
ik
}:

where

*μ*
_{
ib
} is a study-specific baseline effect (i.e. the log-odds for the control group in study

*i* with baseline treatment

*b*),

_{
i(bk)
} is a study-specific log-odds ratio,

*d*
_{
(bk)
} is the pooled effect of treatment

*k* relative to treatment

*b* (a quantity usually of interest in a meta-analysis) and

is the between-study variance or heterogeneity parameter. Random effects NMA assumes that intervention effects are exchangeable across the network regardless of whether or not treatments

*b* and

*k* are included in study

*i*[

18]. This assumption implies that the pooled effects

*d*
_{
(bk)
}, can be expressed as functions of basic parameters with reference to a common comparator or baseline treatment (i.e.

*d*
_{(bk)} =

*d*
_{(Ak)} −

*d*
_{(Ab)}) [

24]. Throughout this paper, we take usual care (UC) intervention to be the reference or ‘baseline’ treatment (i.e. UC is taken as treatment

*A* of equation (

2) above). Multi-arm studies (i.e. studies with more than 2 treatment groups) present a special problem in network meta-analysis because they produce evidence on multiple treatment effects that are correlated through sharing a common reference or ‘baseline’ treatment. Under a homogenous variance assumption (

), the covariance between any two effects that share a common reference treatment is

[

20]. The homogeneous variance assumption allows for the distribution of effects (in a study with an arbitrary number of arms) to be expressed as a univariate marginal distribution and a series of univariate conditional distributions. Specifically, for the

*i*th-study with

*p + 1* arms and

*p* treatment effect estimates relative to the reference treatment, if

then the marginal and conditional univariate distributions for arm

*j*, given the previous 1, ⋯, (

*j* − 1) arms are:

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 study-specific baseline effects, *μ*
_{
ib
} and a Uniform (0,2) prior distribution for the between-study 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 within-study correlations are taken into account. We extend their method to account for the between-study correlation as well.

Note that in Ades

*et al.* (2010), a multinomial likelihood was appropriate as the three binary outcomes (relapse during treatment for Schizophrenia, discontinuation because of intolerable side effects, and discontinuation for other reasons) are mutually exclusive and event probabilities sum to 1 across outcomes. A multinomial likelihood will not be appropriate for our example dataset because each household can have one, two or all three outcomes simultaneously so that the event probabilities do not sum to 1 across outcomes. Instead, we use a Normal distribution on the log-odds scale to take account of the within-study correlations between outcomes. We assume that in each study

*i* and for each

*k*-th arm, the estimates

of the observed log-odds of an event on the

*m*th outcome (

*m* = 1, 2, ⋯,

*M*) jointly follow a multivariate normal distribution:

where (*μ*
_{
ib1}, *μ*
_{
ib2}, ⋯, *μ*
_{
ibM
}) and (*δ*
_{
i(bk)1}, *δ*
_{
i(bk)2}, ⋯, *δ*
_{
i(bk)M
}) represent vectors of *‘true’* baseline and study-specific effects in study *i* with baseline treatment *b* respectively. The quantities
and (*θ*
_{
ik1}, *θ*
_{
ik2}, ⋯, *θ*
_{
ikM
}) represent vectors of observed and ‘*true’* log-odds of response in arm *k* of study *i* and **S**
_{
ik
} is the associated within-study covariance matrix usually assumed known but estimated in practice from the data [35]. We calculated
and the diagonal elements of **S**
_{
ik
} using standard formulae for log-odds and variance of the log-odds [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 off-diagonals of **S**
_{
ik
} were calculated from estimates of within-study correlations
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 study-specific 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.

We refer to equation (

5) as the within-study model and the model describing the distribution of the

*‘true’* effects across studies (presented below) as the between-study model following standard terminology in multivariate meta-analysis [

5,

6,

10,

11,

38,

39]. For the network of two-arm trials, the between-study model for the

*i*th study is thus given by:

where the

*‘true’* effects

*δ*
_{
i(bk)m
} (

*m* = 1, 2, ⋯,

*M*) jointly follow a Normal distribution with mean effects

*d*
_{(bk)m
}. The parameters in equation (

6) have the same interpretation as in equation (

2) except that they are now specific to each outcome. The covariance matrix

*Σ*
_{(bk)} contains terms for the between-study variances,

for each outcome

*m* and the between-study correlations

between effects measured by outcome

*m* and

*n* (

*m ≠ n*) specific to each

*k* versus

*b* comparison. Fitting the full model would thus require a large number of possibly multi-arm studies in order to make

**Σ**
_{(bk)} identifiable [

5,

13]. The number of parameters in

*Σ*
_{(bk)}, can however be reduced if reasonable assumptions can be made about the covariance structure. In particular, most practical applications of NMA methods involve the assumption of a common between-study variance across treatment arms, often referred to as a homogenous variance assumption [

18,

40,

41]. Therefore, to simplify

**Σ**
_{(bk)} we make the additional assumption in this context of a common between-study correlation (

) leading to the following simplified between-study covariance structure for two-arm studies:

where, as in the univariate case, *σ*
_{
m
} represent the common between-study standard deviation or heterogeneity parameter specific to outcome *m*.

To include multi-arm studies in our model, we extend equations (

3) and (

4) to the multiple outcome setting. We show in Appendix A, that under evidence consistency and the homogenous between-study covariance structure (

), equation (

3) can be extended to the multiple outcome settings by formulating the distribution of effects in a multi-arm study

*i* with

*p + 1* arms reporting on

*m* = 1, 2, ⋯,

*M* outcomes as follows:

where

*p* is the number of treatment effect estimates

*.* The corresponding marginal and conditional distributions for arm

*j*, given the previous 1, ⋯, (

*j* − 1) arms are:

To complete model 2,

*μ*
_{
ibm
} and

*d*
_{(1k)m
} are given minimally informative prior distributions:

Prior distributions also need to be specified for

**Σ**
_{(M × M)} which, in general, is non-trivial because of the positive definite constraint. Initially we specified an Inverse-Wishart distribution [

42]:

where

**K** is

*M* ×

*M* scale matrix and

*M* is the total number of outcomes. Specifying minimally informative Inverse-Wishart prior distribution is, however, problematic, especially when the amount of data is small relative to the dimensions of

**Σ**
_{(M × M)} as is the case for our example data. Therefore, to allow for flexibility in formulating a prior distribution for

**Σ**
_{(M × M)}, we also followed a strategy outlined by Lu and Ades (2009) [

43] and more recently by Wei and Higgins (2013) [

39] to express

**Σ**
_{(M × M)} in terms of a diagonal matrix of standard deviations V

^{1/2} and squared positive semi-definite matrix of correlations

**R** based on a separation strategy (Barnard

*et al.*[

44]):

where the off-diagonal 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 between-study 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

From Table

1, it can be seen that none of the studies had considered the interventions E + HSI and E + HV for storage of medicines and storage of other household products. Similarly, interventions E + FE + F and FE were not trialled by any of the included studies on possession of a PCC number. To estimate the full set of 24 basic intervention effects relative to usual care from 9 interventions on 3 outcomes, we applied methods originally proposed by DuMouchel and Harris [

29] and revisited by DuMouchel and Groer [

45] and Jones

*et al.*[

30]. We assume that the pooled effects of treatment

*k* relative to usual care intervention

*d*
_{(Ak)m
}, can be expressed as a sum of treatment-specific effect

*α*
_{
k
} and outcome-specific effect

*γ*
_{
m
}. This assumption replaces the minimally informative prior distribution N(0, 10

^{3}) specified for

*d*
_{(Ak)m
} in model 2 with:

where

*K* is the total number of treatments being evaluated across

*M* outcomes, and for

*k* = 1 (i.e. reference treatment

*A*),

*d*
_{(Ak)m
} equal to zero. Note that on the logarithmic scale, this would imply that the ratio of any intervention effects is constant across outcomes as the

*γ*
_{
m
} cancel, i.e.

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 non-reference/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.

Multi-arm studies are included in model 3 based on equations (

8) and (

9) in the same way as in model 2. To complete model 3, the parameters

*α*
_{
k
},

*γ*
_{
m
} and

*τ* are given minimally informative prior distributions. For the mean effects, this is a normal distribution with zero mean and large variance:

We give

*τ* a Uniform (0, 2) prior distribution, considered to be minimally informative on the log-odds ratio scale. Sensitivity analyses were conducted to assess the impact of specifying alternative prior distributions for

*τ* that are also considered minimally informative [

46]:

- 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 hyper-parameters 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* hyper-parameters (i.e. (*K* − 1) of *α*
_{
k
} and *M* of *γ*
_{
m
} hyper-parameters). 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 non-informative. When large number of data on outcomes is missing, placing informative prior distributions on the hyper-parameters can improve convergence.

### Implementation of the models

We fitted a total of four models, models 1 and 3 as described above and two versions of model 2. In model 2a, we specified an inverse-Wishart prior distribution for the between-study covariance matrix

**Σ**
_{(M × M)} whilst in model 2b, we specified a prior distribution for

**Σ**
_{(M × M)} based on the separation strategy. All four models allowed for multi-arm trials to be included in the analysis. To fit the multivariate NMA models, the quantities

and the diagonals of

**S**
_{
ik
} were estimated using standard 2×2–

*table* formulae [

2]. Next, we obtained estimates of the within-study correlations from the IPD studies using the following three methods: i) Pearson correlation coefficient between the observed outcome events ii) Bootstrapping as described in Daniel and Hughes (1998), and iii) Generalised Estimating Equations (See details in Additional file

1). All three methods produced identical estimates of the correlations between pairs of outcome specific log-odds of event from each IPD study (Additional file

1: Table S1). Therefore, we formulated informative prior distributions for the correlation terms in

**S**
_{
ik
} of equation (

5) using the estimates of the correlations between the observed outcome events (Pearson correlation) as follows:

where
is the within-study correlation between the outcomes *m* and *n* effects measured on the log-odds scale in arm *k* of study *i*, and
and *var* (*r*
^{
mn
})are the mean and variance of the within-study 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 Rubin-Gelman statistic after running 400 000 simulations and discarding the first 200 000 samples as ‘burn in samples’.