A matrix-based method of moments for fitting multivariate network meta-analysis models with multiple outcomes and random inconsistency effects

TL;DR

This paper introduces a new multivariate network meta-analysis model that incorporates multiple outcomes and treatments, allowing for the estimation of heterogeneity and inconsistency effects within a unified framework.

Contribution

It extends existing univariate models to a multivariate network setting and provides an estimation procedure accommodating both heterogeneity and inconsistency.

Findings

Effective estimation of heterogeneity and inconsistency parameters.

Application to real medical data demonstrates model utility.

Extension of DerSimonian and Laird method to multivariate network meta-analysis.

Abstract

Random-effects meta-analyses are very commonly used in medical statistics. Recent methodological developments include multivariate (multiple outcomes) and network (multiple treatments) meta-analysis. Here we provide a new model and corresponding estimation procedure for multivariate network meta-analysis, so that multiple outcomes and treatments can be included in a single analysis. Our new multivariate model is a direct extension of a univariate model for network meta-analysis that has recently been proposed. We allow two types of unknown variance parameters in our model, which represent between-study heterogeneity and inconsistency. Inconsistency arises when different forms of direct and indirect evidence are not in agreement, even having taken between-study heterogeneity into account. However the consistency assumption is often assumed in practice and so we also explain how to fit a reduced model which makes this assumption. Our estimation method extends several other commonly used methods for meta-analysis, including the method proposed by DerSimonian and Laird (1986). We investigate the use of our proposed methods in the context of a real example.

Tables

Table 1. Table 1: Summary of the relapsing remitting multiple sclerosis dataset.

Study	Design	Outcomes
IFNB SG (1)	AB	All three outcomes measured
IFNB SG (2)	AB	All three outcomes measured
Jacobs/Simon	AC	All three outcomes measured
PRISMS (1)	AC	All three outcomes measured
PRISMS (2)	AC	All three outcomes measured
Johnson	AD	Relapse rate and disability progression only
Durelli	BC	Relapse rate and disability progression only
O’Connor (1)	BD	Relapse rate and disability progression only
O’Connor (2)	BD	Relapse rate and disability progression only
Mikol	CD	All three outcomes measured
FREEDOMS 1	AEF	All three outcomes measured
FREEDOMS 2	AEF	All three outcomes measured
TRANSFORMS	CEF	All three outcomes measured

Table 2. Table 2: Treatment effect estimates of each treatment relative to the reference treatment A (placebo).

model	estimate (se)
	AB	AC	AD	AE	AF
MRI ($y_1$)
univariate ($y_1$)	-0.95 (0.39)	-1.00 (0.21)	-0.68 (0.50)	-1.38 (0.26)	-1.52 (0.26)
bivariate ($y_1,y_2$)	-0.94 (0.39)	-1.00 (0.21)	-0.68 (0.50)	-1.39 (0.26)	-1.53 (0.26)
bivariate ($y_1,y_3$)	-0.96 (0.39)	-0.98 (0.22)	-0.66 (0.50)	-1.38 (0.26)	-1.51 (0.26)
trivariate ($y_1,y_2,y_3$)	-0.96 (0.39)	-0.97 (0.22)	-0.67 (0.50)	-1.38 (0.26)	-1.51 (0.26)
Relapse rate ($y_2$)
univariate ($y_2$)	-0.35 (0.10)	-0.25 (0.09)	-0.34 (0.11)	-0.81 (0.12)	-0.78 (0.12)
bivariate ($y_1,y_2$)	-0.35 (0.10)	-0.25 (0.09)	-0.34 (0.11)	-0.81 (0.12)	-0.78 (0.12)
bivariate ($y_2,y_3$)	-0.36 (0.11)	-0.23 (0.10)	-0.33 (0.12)	-0.80 (0.13)	-0.77 (0.13)
trivariate ($y_1,y_2,y_3$)	-0.36 (0.11)	-0.23 (0.10)	-0.33 (0.12)	-0.80 (0.13)	-0.77 (0.13)
Disability progression ($y_3$)
univariate ($y_3$)	-0.46 (0.25)	-0.11 (0.21)	-0.42 (0.25)	-0.33 (0.25)	-0.37 (0.24)
bivariate ($y_2,y_3$)	-0.47 (0.25)	-0.10 (0.21)	-0.43 (0.25)	-0.37 (0.25)	-0.37 (0.25)
bivariate ($y_1,y_3$)	-0.46 (0.25)	-0.11 (0.21)	-0.42 (0.25)	-0.34 (0.25)	-0.38 (0.25)
trivariate ($y_1,y_2,y_3$)	-0.47 (0.25)	-0.10 (0.21)	-0.43 (0.25)	-0.37 (0.25)	-0.37 (0.25)

Table 3. Table 3: Inconsistency and heterogeneity covariance matrices estimates.

model	${\Sigma }_{\omega 11}$	${\Sigma }_{\omega 12}$	${\Sigma }_{\omega 13}$	${\Sigma }_{\omega 22}$	${\Sigma }_{\omega 23}$	${\Sigma }_{\omega 33}$
univariate ($y_1$)	0.0000
univariate ($y_2$)				0.0115
univariate ($y_3$)						0.0713
bivariate ($y_1,y_2$)	0.0002	0.0017		0.0125
bivariate ($y_1,y_3$)	0.0018		0.0116			0.0741
bivariate ($y_2,y_3$)				0.0161	0.0344	0.0735
trivariate ($y_1,y_2,y_3$)	0.0027	0.0066	0.0143	0.0161	0.0349	0.0756
	${\Sigma }_{\beta 11}$	${\Sigma }_{\beta 12}$	${\Sigma }_{\beta 13}$	${\Sigma }_{\beta 22}$	${\Sigma }_{\beta 23}$	${\Sigma }_{\beta 33}$
univariate ($y_1$)	0.1508
univariate ($y_2$)				0.0043
univariate ($y_3$)						0.0000
bivariate ($y_1,y_2$)	0.1523	-0.0110		0.0047
bivariate ($y_1,y_3$)	0.1526		-0.0191			0.0024
bivariate ($y_2,y_3$)				0.0061	0.0015	0.0004
trivariate ($y_1,y_2,y_3$)	0.1538	-0.0116	-0.0195	0.0059	0.0024	0.0027

Equations

$${\mathbf{Y}}_{di}={\mathbf{Z}}_{(d)}\mathbf{\delta }+{\mathbf{\Theta }}_{di}+{\mathbf{\Omega }}_d+{\mathbf{ϵ}}_{di}$$

$$\mathbf{Y}\sim N(\mathbf{Z}\mathbf{\delta },{\tau }_{\beta }^2{\mathbf{M}}_1+{\tau }_{\omega }^2{\mathbf{M}}_2+\mathbf{S})$$

$${\mathbf{Y}}_{di}={\mathbf{X}}_{(d)}\mathbf{\delta }+{\mathbf{\Theta }}_{di}+{\mathbf{\Omega }}_d+{\mathbf{ϵ}}_{di}$$

$$\mathbf{Y}\sim N(\mathbf{X}\mathbf{\delta },{\mathbf{M}}_1\otimes {\mathbf{\Sigma }}_{\beta }+{\mathbf{M}}_2\otimes {\mathbf{\Sigma }}_{\omega }+\mathbf{S})$$

$$\mathbf{Q}=\mathbf{W}\{\mathbf{R}(\mathbf{Y}-\hat{\mathbf{Y}})\}\{\mathbf{R}(\mathbf{Y}-\hat{\mathbf{Y}}){\}}^T=\mathbf{W}(\mathbf{Y}-\hat{\mathbf{Y}})(\mathbf{Y}-\hat{\mathbf{Y}}{)}^T\mathbf{R}$$

$${\mathbf{Y}}_d\sim N({\mathbf{X}}_d{\mathbf{\beta }}_d,{\mathbf{M}}_1^{\mathbf{d}}\otimes {\mathbf{\Sigma }}_{\beta }+{\mathbf{S}}_d)$$

$${\mathbf{Q}}_d={\mathbf{W}}_d({\mathbf{Y}}_d-{\hat{\mathbf{Y}}}_d)({\mathbf{Y}}_d-{\hat{\mathbf{Y}}}_d{)}^T{\mathbf{R}}_d$$

$$\text{mbox}E[\text{mbox}btr(\mathbf{Q})]=\sum _{i=1}^n\sum _{j=1}^n\sum _{k=1}^n{m}_{1ij}{\mathbf{A}}_{k,i}{\mathbf{\Sigma }}_{\beta }{\mathbf{B}}_{j,k}+$$

$$\sum _{i=1}^n\sum _{j=1}^n\sum _{k=1}^n{m}_{2ij}{\mathbf{A}}_{k,i}{\mathbf{\Sigma }}_{\omega }{\mathbf{B}}_{j,k}+\text{mbox}btr(\mathbf{B}).$$

$$\text{mbox}vec(\text{mbox}E[\text{mbox}btr(\mathbf{Q})])=\mathbf{C}\text{mbox}vec({\mathbf{\Sigma }}_{\beta })+\mathbf{D}\text{mbox}vec({\mathbf{\Sigma }}_{\omega })+\mathbf{E}$$

$$\mathbf{C}=\sum _{i=1}^n\sum _{j=1}^n\sum _{k=1}^n{m}_{1ij}{\mathbf{B}}_{j,k}^T\otimes {\mathbf{A}}_{k,i}$$

$$\mathbf{D}=\sum _{i=1}^n\sum _{j=1}^n\sum _{k=1}^n{m}_{2ij}{\mathbf{B}}_{j,k}^T\otimes {\mathbf{A}}_{k,i}$$

$$\mathbf{E}=\text{mbox}vec(\text{mbox}btr(\mathbf{B})).$$

$$\text{mbox}vec\left(\text{mbox}E[\sum _{d=1}^D\text{mbox}btr({\mathbf{Q}}_d)]\right)=\left(\sum _{d=1}^D{\mathbf{C}}_d\right)\text{mbox}vec({\mathbf{\Sigma }}_{\beta })+\sum _{d=1}^D{\mathbf{E}}_d$$

$${\mathbf{C}}_d=\sum _{i=1}^{n_d}\sum _{j=1}^{n_d}\sum _{k=1}^{n_d}{m}_{1ij}^d{\mathbf{B}}_{d,j,k}^T\otimes {\mathbf{A}}_{d,k,i}$$

$${\mathbf{E}}_d=\text{mbox}vec(\text{mbox}btr({\mathbf{B}}_d)).$$

$$(\mathbf{ACB}{)}_{k,l}=\sum _{i=1}^n\sum _{j=1}^n{\mathbf{A}}_{k,i}{\mathbf{C}}_{i,j}{\mathbf{B}}_{j,l}$$

$$(\mathbf{A}(\mathbf{M}\otimes \mathbf{\Sigma })\mathbf{B}{)}_{k,l}=\sum _{i=1}^n\sum _{j=1}^n{m}_{ij}{\mathbf{A}}_{k,i}\mathbf{\Sigma }{\mathbf{B}}_{j,l}$$

$$\text{mbox}btr(\mathbf{A}(\mathbf{M}\otimes \mathbf{\Sigma })\mathbf{B})=\sum _{i=1}^n\sum _{j=1}^n\sum _{k=1}^n{m}_{ij}{\mathbf{A}}_{k,i}\mathbf{\Sigma }{\mathbf{B}}_{j,k}$$

$$\text{mbox}E[\mathbf{Q}]=\mathbf{A}({\mathbf{M}}_1\otimes {\mathbf{\Sigma }}_{\beta }+{\mathbf{M}}_2\otimes {\mathbf{\Sigma }}_{\omega })\mathbf{B}+\mathbf{B}$$

$$\mathbf{A}=({\mathbf{I}}_{np}-\mathbf{H}{)}^T\mathbf{W}$$

$$\mathbf{B}=({\mathbf{I}}_{np}-\mathbf{H}{)}^T\mathbf{R}$$

$$\text{mbox}E[\text{mbox}btr(\mathbf{Q})]=\sum _{i=1}^n\sum _{j=1}^n\sum _{k=1}^n{m}_{1ij}{\mathbf{A}}_{k,i}{\mathbf{\Sigma }}_{\beta }{\mathbf{B}}_{j,k}+$$

$$\sum _{i=1}^n\sum _{j=1}^n\sum _{k=1}^n{m}_{2ij}{\mathbf{A}}_{k,i}{\mathbf{\Sigma }}_{\omega }{\mathbf{B}}_{j,k}+\text{mbox}btr(\mathbf{B})$$

$${\mathbf{H}}_d={\mathbf{X}}_d({\mathbf{X}}_d^T{\mathbf{W}}_d{\mathbf{X}}_d{)}^{-1}{\mathbf{X}}_d^T{\mathbf{W}}_d$$

$${\mathbf{A}}_d=({\mathbf{I}}_{p{n}_d}-{\mathbf{H}}_d{)}^T{\mathbf{W}}_d$$

$${\mathbf{B}}_d=({\mathbf{I}}_{p{n}_d}-{\mathbf{H}}_d{)}^T{\mathbf{R}}_d$$

$$\text{mbox}vec(\text{mbox}E[\text{mbox}btr({\mathbf{Q}}_d)])={\mathbf{C}}_d\text{mbox}vec({\mathbf{\Sigma }}_{\beta })+{\mathbf{E}}_d$$

$${\mathbf{C}}_d=\sum _{i=1}^{n_d}\sum _{j=1}^{n_d}\sum _{k=1}^{n_d}{m}_{1ij}^d{\mathbf{B}}_{d,j,k}^T\otimes {\mathbf{A}}_{d,k,i}$$

$${\mathbf{E}}_d=\text{mbox}vec(\text{mbox}btr({\mathbf{B}}_d))$$

Full text

A matrix-based method of moments for fitting multivariate network meta-analysis models

with multiple outcomes

and random inconsistency effects

Dan Jackson

MRC Biostatistics Unit, Cambridge, UK

email: daniel.jackson@mrc-bsu.cam.ac.uk

Sylwia Bujkiewicz

Biostatistics Research Group, Department of Health Sciences, University of Leicester, UK

Martin Law

MRC Biostatistics Unit, Cambridge, UK

Richard D Riley

Centre for Prognosis Research, Research Institute for Primary Care and Health Sciences, University of Keele, UK

Ian White

MRC Biostatistics Unit, Cambridge, UK

MRC Clinical Trials Unit at University College London, UK

Random-effects meta-analyses are very commonly used in medical statistics. Recent methodological developments include multivariate (multiple outcomes) and network (multiple treatments) meta-analysis. Here we provide a new model and corresponding estimation procedure for multivariate network meta-analysis, so that multiple outcomes and treatments can be included in a single analysis. Our new multivariate model is a direct extension of a univariate model for network meta-analysis that has recently been proposed. We allow two types of unknown variance parameters in our model, which represent between-study heterogeneity and inconsistency. Inconsistency arises when different forms of direct and indirect evidence are not in agreement, even having taken between-study heterogeneity into account. However the consistency assumption is often assumed in practice and so we also explain how to fit a reduced model which makes this assumption. Our estimation method extends several other commonly used methods for meta-analysis, including the method proposed by DerSimonian and Laird (1986). We investigate the use of our proposed methods in the context of a real example.

1 Introduction

Meta-analysis, the statistical process of pooling the results from separate studies, is commonly used in medical statistics and now requires little introduction. The univariate random-effects model is often used for this purpose. This model has recently been extended to the multivariate (multiple outcomes; Jackson et al., 2011) and network (multiple treatments; Lu and Ades, 2004) meta-analysis settings. In a network meta-analysis, more than two treatments are included in the same analysis. The main advantage of network meta-analysis is that, by using indirect information contained in the network, more precise and coherent inference is possible, especially when direct evidence for particular treatment comparisons is limited. Here we describe a new model that extends the random-effects modelling framework to the multivariate network meta-analysis setting, so that both multiple outcomes and multiple treatments may be included in the same analysis.

Other multivariate extensions of univariate methods for network meta-analysis have previously been proposed. For example, Achana et al. (2014) analyse multiple correlated outcomes in multi-arm studies in public health. Efthimiou et al. (2014) propose a model for the joint modelling of odds ratios on multiple endpoints. Efthimiou et al. (2015) develop another model that is a network extension of an alternative multivariate meta-analytic model that was originally proposed by Riley et al. (2008). A network meta-analysis of multiple outcomes with individual patient data has also been proposed by Hong et al. (2015) under both contrast-based and arm-based parameterizations, and Hong et al. (2016) develop a Bayesian framework for multivariate network meta-analysis. These multivariate network meta-analysis models are based on the assumption of consistency in the network, extending the approach introduced by Lu and Ades (2004). In contrast to these previously developed methods, the method proposed here relaxes the consistency assumption. This assumption is sometimes found to be false across the entire network (Veroniki et al., 2013). We model the inconsistency using a design-by-treatment interaction, so that different forms of direct and indirect evidence may not agree, even after taking between-study heterogeneity into account. However we assume that the design-by-interaction terms follow normal distributions, and so conceptualise inconsistency as another source of random variation. This allows us to achieve the dual aim of estimating meaningful treatment effects whilst also allowing for inconsistency in the network.

Although we allow inconsistency in the network, we propose a relatively simple model. Our preference for a simple model is because the between-study covariance structure is typically hard to identify accurately in multivariate meta-analyses (Jackson et al., 2011) and also because network meta-analysis datasets are usually small (Nikolakopoulou et al., 2014). The new model that we propose for multivariate network meta-analysis is a direct generalisation of the univariate network meta-analysis model proposed by Jackson et al. (2016), which is a particular form of the design-by-treatment interaction model (Higgins et al., 2012). In addition to proposing a new model for multivariate network meta-analysis, we also develop a corresponding new estimation method. This estimation method is based on the method of moments and extends a wide variety of related methods. In particular, we extend the estimation method described by DerSimonian and Laird (1986) by directly extending the matrix based extension of DerSimonian and Laird’s estimation method for multivariate meta-analysis (Chen et al., 2012; Jackson et al., 2013). We adopt the usual two-stage approach to meta-analysis, where the estimated study-specific treatment effects (including the within-study covariance matrices) are computed in the first stage. We give some information about how this first stage is performed but our focus is the second stage, where the meta-analysis model is fitted.

The paper is set out as follows. In section 2, we briefly describe the univariate model for network meta-analysis to motivate our new multivariate network meta-analysis model in section 3. We present our new estimation method in section 4 and we apply our methods to a real dataset in section 5. We conclude with a short discussion in section 6.

2 A univariate network meta-analysis model

Here we describe our univariate modelling framework for network meta-analysis (Jackson et al., 2016; Law et al., 2016). Without loss of generality, we take treatment A as the reference treatment for the network meta-analysis. The other treatments are B, C, etc. We take the design $d$ as referring only to the set of treatments compared in a study. For example, if the first design compares treatments A and B only, then $d=1$ refers to two-arm studies that compare these two treatments. We define $t$ to be the total number of treatments included in the network, and $t_{d}$ to be the number of treatments included in design $d$. We define $D$ to be the number of different designs, $N_{d}$ to be the number of studies of design $d$, and $N=\sum\limits_{d=1}^{D}N_{d}$ to be the total number of studies. We will use the word ‘contrast’ to refer to a particular treatment comparison or effect in a particular study, for example the ‘AB contrast’ in the first study.

We model the estimated relative treatment effects, rather than the average outcomes in each arm, and so perform contrast based analyses. We define ${\bm{Y}}_{di}$ to be the $c_{d}\times 1$ column vector of estimated relative treatment effects from the $i$th study of design $d$, where $c_{d}=t_{d}-1$. We define $n_{d}=N_{d}c_{d}$ to be the total number of estimated treatment effects that design $d$ contributes to the analysis, and $n=\sum\limits_{d=1}^{D}n_{d}$ to be the total number of estimated treatment effects that contribute to the analysis. To specify the outcome data ${\bm{Y}}_{di}$, we choose a baseline treatment in each design $d$. The entries of ${\bm{Y}}_{di}$ are then the estimated effects of the other $c_{d}$ treatments included in design $d$ relative to this baseline treatment. For example, if we take $d=2$ to indicate the ‘CDE design’ then $c_{2}=2$. Taking C as the baseline treatment for this design, the two entries of the ${\bm{Y}}_{2i}$ vectors are the estimated relative effects of treatment D compared to C and of treatment E compared to C. For example, the entries of the ${\bm{Y}}_{di}$ could be estimated log-odds ratios or mean differences.

We use normal approximations for the within-study distributions. We define ${\bm{S}}_{di}$ to be the $c_{d}\times c_{d}$ within-study covariance matrix corresponding to ${\bm{Y}}_{di}$. We treat all ${\bm{S}}_{di}$ as fixed and known in analysis. Ignoring the uncertainty in the ${\bm{S}}_{di}$ is acceptable provided that the studies are reasonably large and is conventional in meta-analysis, but this approximation is motivated by pragmatic considerations because this greatly simplifies the modelling. We do not impose any constraints on the form of ${\bm{S}}_{di}$ other than they must be valid covariance matrices. The lead diagonal entries of the ${\bm{S}}_{di}$ are within-study variances that can be calculated using standard methods. Assuming that the studies are composed of independent samples for each treatment, the other entries of the ${\bm{S}}_{di}$ are calculated as the variance of the average outcome (for example the log odds or the sample mean) of the baseline treatment.

We define $\delta_{1}^{AB}$, $\delta_{1}^{AC}$, $\cdots$, $\delta_{1}^{AZ}$, where Z is the final treatment in the network, to be treatment effects relative to the reference treatment A, and call them basic parameters (Lu and Ades, 2006). We use the subscript 1 when defining the basic parameters to emphasise that they are treatment effects for the first (and in this section, only) outcome. We define $c=t-1$ to be the number of basic parameters in the univariate setting. Treatment effects not involving A can be obtained as linear combinations of the basic parameters and are referred to as functional parameters (Lu and Ades, 2006). For example the average treatment effect of treatment E to treatment C, $\delta_{1}^{CE}=\delta_{1}^{AE}-\delta_{1}^{AC}$, is a functional parameter. We define the $c\times 1$ column vector ${\bm{\delta}}=(\delta_{1}^{AB},\delta_{1}^{AC},\cdots,\delta_{1}^{AZ})^{T}$ and design specific $c_{d}\times c$ design matrices ${\mathbf{Z}}_{(d)}$. We use the subscript $(d)$ in these design matrices to emphasise that they apply to each individual study of design $d$; we reserve the subscript $d$ for design matrices that describe regression models for all outcome data from this design. If the $i$th entry of the ${\mathbf{Y}}_{di}$ are estimated treatment effects of treatment J relative to the reference treatment A then the $i$th row of ${\mathbf{Z}}_{(d)}$ contains a single nonzero entry: 1 in the $(j-1)$th column, where $j$ is the position of $J$ in the alphabet. If instead the $i$th entry of the ${\mathbf{Y}}_{di}$ are estimated treatment effects of treatment $J$ relative to treatment $K$, $K\neq A$, then the $i$th row of ${\mathbf{Z}}_{(d)}$ contains two nonzero entries: 1 in the $(j-1)$th column and -1 in the $(k-1)$th column.

Our univariate model for network meta-analysis is

[TABLE]

where ${\bm{\Theta}}_{di}\sim N({\bf 0},\tau^{2}_{\beta}\mathbf{P}_{c_{d}})$, ${\bm{\Omega}}_{d}\sim N({\bf 0},\tau^{2}_{\omega}\mathbf{P}_{c_{d}})$, ${\bm{\epsilon}}_{di}\sim N({\bf 0},\mathbf{S}_{di})$, all ${\bm{\Theta}}_{di}$, ${\bm{\Omega}}_{di}$ and ${\bm{\epsilon}}_{di}$ are independent, and ${\bf P}_{c_{d}}$ is the $c_{d}\times c_{d}$ matrix with ones on the leading diagonal and halves elsewhere. We refer to $\tau^{2}_{\beta}$ and $\tau^{2}_{\omega}$ as the between-study variance, and the inconsistency variance, respectively. The term ${\bm{\Theta}}_{di}$ is a study-by-treatment interaction term that models between-study heterogeneity. The model ${\bm{\Theta}}_{di}\sim N({\bf 0},\tau^{2}_{\beta}\mathbf{P}_{c_{d}})$ implies that the heterogeneity variance is the same for all contrasts for every study regardless of whether or not the comparison is relative to the baseline treatment (Lu and Ades, 2004). Other simple choices of $\mathbf{P}_{c_{d}}$, such as allowing the off-diagonal entries to differ from 0.5, violate this symmetry between treatments. For example, in the case $d=2$ indicating the CDE design, the between-study variances for the CD and CE effects in this study are given by the two diagonal entries of $\tau^{2}_{\beta}\mathbf{P}_{c_{d}}$, which are both $\tau^{2}_{\beta}$. The between-study variance for the effect of E relative to D is $(-1,1)\tau^{2}_{\beta}\mathbf{P}_{c_{d}}(-1,1)^{T}$, which is also $\tau^{2}_{\beta}$. The ${\bm{\Omega}}_{d}$ are design-by-treatment interaction terms that model inconsistency in the network. The model ${\bm{\Omega}}_{d}\sim N({\bf 0},\tau^{2}_{\omega}\mathbf{P}_{c_{d}})$ implies that the inconsistency variance is the same for all contrasts for every design; other simple choices of $\mathbf{P}_{c_{d}}$ also violate this symmetry.

To describe all estimates from all studies, we stack the ${\bm{Y}}_{di}$ from the same design to form the $n_{d}\times 1$ column vector ${\bm{Y}}_{d}=({\bm{Y}}_{d1}^{T},\cdots,{\bm{Y}}^{T}_{dN_{d}})^{T}$, and we then stack these ${\bm{Y}}_{d}$ to form the $n\times 1$ column vector ${\bm{Y}}=({\bm{Y}}_{1}^{T},\cdots,{\bm{Y}}^{T}_{D})^{T}$. Jackson et al. (2016) then use three further matrices that we also define here because they will be required to describe the estimation procedure that follows. The matrix ${\bf M_{1}}$ is defined as a $n\times n$ square matrix where $m_{1ij}=0$ if the $i$th and $j$th entries of ${\bm{Y}}$, $i,j=1,\cdots n$, are estimates from different studies; otherwise ${m}_{1ii}=1$, and ${m}_{1ij}=1/2$ for $i\neq j$. The matrix ${\bf M_{2}}$ is defined as a $n\times n$ square matrix where ${m}_{2ij}=0$ if the $i$th and $j$th entries of ${\bm{Y}}$, $i,j=1,\cdots n$, are estimates from different designs; otherwise ${m}_{2ij}=1$ if the $i$th and $j$th entries of ${\bm{Y}}$ are estimates of the same treatment comparison (for example, treatment A compared to treatment B) and ${m}_{2ij}=1/2$ if these entries are estimates of different treatment comparisons. The supplementary materials show a concrete example showing how these two matrices are formed. Jackson et al. (2016) also define a $n\times c$ univariate design matrix ${\mathbf{Z}}$, where if the $i$th entry of ${\mathbf{Y}}$ is an estimated treatment effect of treatment J relative to the reference treatment A then the $i$th row of ${\mathbf{Z}}$ contains a single nonzero entry: 1 in the $(j-1)$th column, where $j$ is the position of $J$ in the alphabet. If instead the $i$th entry ${\mathbf{Y}}$ is an estimated treatment effect of treatment $J$ relative to treatment $K$, $K\neq A$, then the $i$th row of ${\mathbf{Z}}$ contains two nonzero entries: 1 in the $(j-1)$th column and -1 in the $(k-1)$th column. Defining ${\mathbf{S}}_{d}=\mbox{diag}({\mathbf{S}}_{d1},\cdots,{\mathbf{S}}_{dN_{d}})$, and then ${\mathbf{S}}=\mbox{diag}({\mathbf{S}}_{1},\cdots,{\mathbf{S}}_{D})$, model (1) can be presented for the entire dataset as

[TABLE]

3 A multivariate network meta-analysis model

We now explain how to extend the univariate model in section 2 to the multivariate setting to handle multiple outcomes. We define $p$ to be the number of outcomes, and so the dimension of the network meta-analysis, so that we now consider the case where $p>1$. The ${\bm{Y}}_{di}$ are now $pc_{d}\times 1$ column vectors, where the ${\bm{Y}}_{di}$ contain $c_{d}$ column vectors of length $p$. For example, in a $p=5$ dimensional meta-analysis and continuing with the example where $d=2$ indicates the CDE design, we have $c_{2}=2$. The ${\bm{Y}}_{2i}$ are then $10\times 1$ column vectors where, taking C as the baseline treatment for this design, the first five entries of the ${\bm{Y}}_{2i}$ are estimated relative treatment effect of D compared to C and the second five entries are the same estimate of E compared to C. We define the $pc\times 1$ column vector ${\bm{\delta}}=(\delta_{1}^{AB},\delta_{1}^{AC},\cdots,\delta_{1}^{AZ},\delta_{2}^{AB},\delta_{2}^{AC},\cdots,\delta_{2}^{AZ},\cdots,\delta_{p}^{AB},\delta_{p}^{AC},\cdots,\delta_{p}^{AZ})^{T}$, so that this vector contains the basic parameters for each outcome in turn. When $p=1$ the vector ${\bm{\delta}}$ reduces to its definition in the univariate setting, as given in section 2.

We define ${\bf\Sigma}_{\beta}$ and ${\bf\Sigma}_{\omega}$ to be $p\times p$ unstructured covariance matrices that are multivariate generalisations of $\tau^{2}_{\beta}$ and $\tau^{2}_{\omega}$. These two matrices contain the between-study variances and covariances, and the inconsistency variances and covariances, respectively, for all $p$ outcomes. We refer to ${\bf\Sigma}_{\beta}$ and ${\bf\Sigma}_{\omega}$ as the between-study covariance matrix, and the inconsistency covariance matrix, respectively. We continue to treat the within-study covariance matrices ${\bm{S}}_{di}$ as if fixed and known in analysis but these are now $pc_{d}\times pc_{d}$ matrices. The entries of the ${\bm{S}}_{di}$ matrices that describe the covariance of estimated treatment effects for the same outcome can be obtained as in the univariate setting. However the other entries of ${\bm{S}}_{di}$, that describe the covariance between treatment effects for different outcomes, are harder to obtain in practice. A variety of strategies for dealing with this difficulty have been proposed (Jackson et al., 2011; Wei and Higgins, 2013).

3.1 The proposed multivariate model for network meta-analysis

In the multivariate setting, to allow correlations between estimated treatment effects for different outcomes, both within studies and designs, we propose that model (1) is generalised to

[TABLE]

where ${\bm{X}}_{(d)}=(({\bm{I}}_{p}\otimes{\bm{Z}}_{(d)1})^{T},\cdots,({\bm{I}}_{p}\otimes{\bm{Z}}_{(d)c_{d}})^{T})^{T}$, ${\bm{Z}}_{(d)i}$ is the $i$th row of ${\bm{Z}}_{(d)}$, ${\bm{\Theta}}_{di}\sim N({\bf 0},\mathbf{P}_{c_{d}}\otimes{\bf\Sigma}_{\beta})$, ${\bm{\Omega}}_{d}\sim N({\bf 0},\mathbf{P}_{c_{d}}\otimes{\bf\Sigma}_{\omega})$ and ${\bm{\epsilon}}_{di}\sim N({\bf 0},{\bm{S}}_{di})$, where all ${\bm{\Theta}}_{di}$, ${\bm{\Omega}}_{d}$ and ${\bm{\epsilon}}_{di}$ are independent, and $\otimes$ is the Kronecker product. The random ${\bm{\Theta}}_{di}$ and ${\bm{\Omega}}_{d}$ continue to model between-study heterogeneity, and inconsistency, respectively. Recalling that ${\bm{\delta}}$ contains the basic parameters for each outcome in turn, the design matrices ${\bm{X}}_{(d)}$ provide the correct linear combinations of basic parameters to describe the mean of all estimated treatment effects in ${\bm{Y}}_{di}$. Model (2) reduces to model (1) in one dimension. The definition of $\mathbf{P}_{c_{d}}$ means that ${\bf\Sigma}_{\beta}$ and ${\bf\Sigma}_{\omega}$ are the between-study covariance matrix, and inconsistency covariance matrix, for all contrasts. We continue define ${\bm{Y}}$ as in the univariate setting, where ${\bm{Y}}$ contains $n$ column vectors of estimated treatment effects that are of length $p$, so that ${\bm{Y}}$ is a $np\times 1$ column vector in the multivariate setting. We define the multivariate $np\times pc$ design matrix ${\mathbf{X}}=(({\mathbf{I}}_{p}\otimes{\mathbf{Z}}_{1})^{T},\cdots,({\mathbf{I}}_{p}\otimes{\mathbf{Z}}_{n})^{T})^{T}$, where ${\mathbf{Z}}_{i}$ is the $i$th row of ${\mathbf{Z}}$. Model (2) can be presented for the entire dataset as

[TABLE]

where we continue to define ${\mathbf{S}}$ as in the univariate case. Matrices ${\bf M_{1}}$ and ${\bf M_{2}}$ are the same as in the univariate setting, and so continue to be $n\times n$ matrices. Model (3) is a linear mixed model for network meta-analysis and is conceptually similar to other models of this type (Piepho et al., 2012). If ${\bf\Sigma}_{\omega}={\bf 0}$ then all ${\bm{\Omega}}_{d}={\bm{0}}$ and there is no inconsistency; we refer to this reduced model as the ‘consistent model’. If both ${\bf\Sigma}_{\beta}={\bf 0}$ and ${\bf\Sigma}_{\omega}={\bf 0}$ then all studies estimate the same effects to within-study sampling error and we refer to this model as the ‘common-effect and consistent model’.

Missing data (unobserved entries of ${\mathbf{Y}}$) are common in applications as not all studies may provide data for all outcomes and contrasts. When there are missing outcome data, the model for the observed data is the marginal model for the observed data implied by (3), where any rows of ${\mathbf{Y}}$ that contain missing values are discarded. We will use a non-likelihood based approach for making inferences and so assume any data are missing completely at random (Seaman et al., 2013). We define the diagonal $np\times np$ missing indicator matrix ${\mathbf{R}}$, where ${\mathbf{R}}_{ii}=1$ if ${\mathbf{Y}}_{i}$ is observed, ${\mathbf{R}}_{ii}=0$ if ${\mathbf{Y}}_{i}$ is missing, and ${\mathbf{R}}_{ij}=0$ if $i\neq j$.

4 Multivariate estimation: a new method of moments

Our estimation procedure is motivated by the univariate method proposed by DerSimonian and Laird (1986). This was developed in the much simpler setting where each study provides a single estimate $Y_{i}$, and where the random-effects model $Y_{i}\sim N(\delta,\tau^{2}+S_{i})$ is assumed. This estimation method for $\tau^{2}$ uses the $Q$ statistic, where $Q=\sum S_{i}^{-1}(Y_{i}-\hat{\delta})^{2}$ and $\hat{\delta}=\sum S_{i}^{-1}Y_{i}/\sum S_{i}^{-1}$ is the pooled estimate under the common-effect model ($\tau^{2}=0$).

Now consider an alternative representation of this $Q$ statistic. Taking ${\bf Y}=(Y_{1},\cdots,Y_{n})^{T}$, ${\bf S}=\mbox{diag}(S_{1},\cdots,S_{n})$ and ${\bf W}={\bf S}^{-1}$ means that $Q=\mbox{tr}({\mathbf{W}}({\mathbf{Y}}-\hat{{\mathbf{Y}}})({\mathbf{Y}}-\hat{{\mathbf{Y}}})^{T})$, where $\hat{{\mathbf{Y}}}$ is obtained under the common-effect model. To obtain a $p\times p$ matrix generalisation of $Q$ for multivariate analyses, we replace the trace operator with the block trace operator in this expression (Jackson et al., 2013). The block trace operator is a generalisation of the trace that sums over all $n$ of the $p\times p$ matrices along the main block diagonal of an $np\times np$ matrix. This produces a $p\times p$ matrix. In the absence of missing data we can write our multivariate generalisation of the $Q$ statistic, $\mbox{btr}({\mathbf{W}}({\mathbf{Y}}-\hat{{\mathbf{Y}}})({\mathbf{Y}}-\hat{{\mathbf{Y}}})^{T})$, as a weighted sum of outer products of $p\times 1$ vectors of residuals under the common-effect and consistent model. Hence the distribution of $\mbox{btr}({\mathbf{W}}({\mathbf{Y}}-\hat{{\mathbf{Y}}})({\mathbf{Y}}-\hat{{\mathbf{Y}}})^{T})$ depends directly on the magnitudes of unknown variance components.

4.1 A Q matrix for multivariate network meta-analysis

We define a within-study precision matrix ${\mathbf{W}}$ corresponding to ${\mathbf{S}}$. If there are no missing outcome data in ${\mathbf{Y}}$ then we define ${\mathbf{W}}={\mathbf{S}}^{-1}$, where ${\mathbf{S}}$ is taken from model (3). If there are missing data in ${\mathbf{Y}}$ then the entries of ${\mathbf{W}}$ that correspond to observed data are obtained as the inverse of the corresponding entries of the within-study covariance matrix of reduced dimension (equal to that of the observed data) and the other entries of ${\mathbf{W}}$ are set to zero. For example, consider the case where ${\mathbf{Y}}$ is a $6\times 1$ vector but only the second and fifth entries are observed; this corresponds to much less outcome data than would be used in practice but provides an especially simple example. Then we define ${\mathbf{S}}_{r}$, where the subscript $r$ indicates a dimension reduction, as a $2\times 2$ matrix whose entries are the within variances and covariances of the two observed entries of ${\mathbf{Y}}$. The $6\times 6$ precision matrix ${\mathbf{W}}$ then has all zero entries in the first, third, fourth and sixth rows and columns. However the remaining entries of ${\mathbf{W}}$ are the entries of the $2\times 2$ matrix ${\mathbf{S}}_{r}^{-1}$, so that ${\mathbf{W}}_{22}=({\mathbf{S}}_{r}^{-1})_{11}$, ${\mathbf{W}}_{25}=({\mathbf{S}}_{r}^{-1})_{12}$, ${\mathbf{W}}_{52}=({\mathbf{S}}_{r}^{-1})_{21}$, and ${\mathbf{W}}_{55}=({\mathbf{S}}_{r}^{-1})_{22}$. We define $\hat{{\mathbf{Y}}}$ to be the fitted value of ${{\mathbf{Y}}}$ under the common-effect and consistent model (${\bf\Sigma}_{\beta}={\bf\Sigma}_{\omega}={\bf{0}}$), so that $\hat{{\mathbf{Y}}}={\mathbf{H}}{\mathbf{Y}}$ where ${\mathbf{H}}={\mathbf{X}}({\mathbf{X}}^{T}{\mathbf{W}}{\mathbf{X}})^{-1}{\mathbf{X}}^{T}{\mathbf{W}}$. We also define an asymmetric $np\times np$ matrix (Jackson et al., 2013)

[TABLE]

Our definitions of ${\mathbf{W}}$ and ${\mathbf{R}}$ mean that ${\mathbf{W}}{\mathbf{R}}={\mathbf{W}}$, which results in the simplified version of ${\mathbf{Q}}$ in (4). From the first form given in (4), we have that the residuals ${\mathbf{Y}}-\hat{\mathbf{Y}}$ are pre-multiplied by ${\mathbf{R}}$, so that any residuals that correspond to missing outcome data do not contribute to ${\mathbf{Q}}$. Furthermore missing outcome data do not contribute to $\hat{\mathbf{Y}}$ because they have no weight under the common-effect and consistent model. Hence we can impute missing outcome data with any finite value without changing the value of ${\mathbf{Q}}$. This is merely a convenient way to handle missing data numerically and has no implications for the statistical modelling.

4.2 Design specific Q matrices for multivariate network meta-analysis

In order to identify the full model, we will require design-specific versions of ${\mathbf{Q}}$ that only use data from a particular design. As in the univariate setting, we stack the outcome data from design $d$ to form the vector ${\mathbf{Y}}_{d}=({\mathbf{Y}}^{T}_{d1},\cdots,{\mathbf{Y}}^{T}_{dN_{d}})^{T}$. In the multivariate setting the vector ${\mathbf{Y}}_{d}$ contains $n_{d}$ estimated effects each of length $p$, so that ${\mathbf{Y}}_{d}$ is now a $pn_{d}\times 1$ column vector. We define the design specific $n_{d}\times n_{d}$ matrix ${\bf M_{1}^{d}}$, where ${m}^{d}_{1ij}=0$ if the $i$th and $j$th estimated effect (of length $p$) in ${\mathbf{Y}}_{d}$, $i,j=1,\cdots,n_{d}$, are from separate studies; otherwise ${m}^{d}_{1ii}=1$ and ${m}^{d}_{1ij}=1/2$ for $i\neq j$. We define the $pn_{d}\times pc_{d}$ design matrix ${\mathbf{X}}_{d}$ which is obtained by stacking identity matrices of dimension $pc_{d}$, where we include one such identity matrix for each study of design $d$. Hence ${\mathbf{X}}_{d}={\bf 1}_{N_{d}}\otimes{\mathbf{I}}_{pc_{d}}$, where ${\bf 1}_{N_{d}}$ is the $N_{d}\times 1$ column vector where every entry is one. We also define the $pc_{d}\times 1$ column vector ${\bm{\beta}}_{d}={\bm{X}}_{(d)}{\bm{\delta}}+{\bm{\Omega}}_{d}$.

An identifiable design-specific marginal model for outcome data from design $d$ only, that is implied by model (2), is

[TABLE]

where ${\mathbf{S}}_{d}=\mbox{diag}({\mathbf{S}}_{d1},\cdots,{\mathbf{S}}_{dN_{d}})$. We can also calculate design specific versions of (4) where we calculate all quantities, including the fitted values, using just the data from studies of design $d$. We define these $pn_{d}\times pn_{d}$ design specific matrices as

[TABLE]

where ${\mathbf{W}}_{d}$, ${\mathbf{R}}_{d}$ and $\hat{{\mathbf{Y}}}_{d}$ in (6) are defined in the same way as ${\mathbf{W}}$, ${\mathbf{R}}$ and $\hat{{\mathbf{Y}}}$ in (4) but where only data from design $d$ are used. Hence ${\mathbf{R}}_{d}$ and ${\mathbf{W}}_{d}$ are the missing indicator matrix, and the within-study precision matrix, of ${\mathbf{Y}}_{d}$, respectively. We compute $\hat{{\mathbf{Y}}}_{d}={\mathbf{H}}_{d}{\mathbf{Y}}_{d}$ where ${\mathbf{H}}_{d}={\mathbf{X}}_{d}({\mathbf{X}}_{d}^{T}{\mathbf{W}}_{d}{\mathbf{X}}_{d})^{-1}{\mathbf{X}}_{d}^{T}{\mathbf{W}}_{d}$. When computing ${\mathbf{H}}_{d}$ we take the matrix inverse to be the Moore-Penrose pseudoinverse. This is so that any design-specific regression corresponding to this hat matrix that is not fully identifiable (due to missing outcome data) can still contribute to the estimation. We use model (5) to derive the properties of ${\mathbf{Q}}_{d}$ in equation (6).

4.3 The estimating equations

We base our estimation on the two $p\times p$ matrices $\mbox{btr}({\mathbf{Q}})$ and $\sum\limits_{d=1}^{D}\mbox{btr}({\mathbf{Q}}_{d})$, where ${\mathbf{Q}}$ and ${\mathbf{Q}}_{d}$ are given in (4) and (6), respectively. Specifically, we match these quantities to their expectations to estimate the unknown variance parameters using the method of moments.

4.3.1 Evaluating $\mbox{E}[\mbox{btr}({\mathbf{Q}})]$ and deriving the first estimating equation

We define ${\mathbf{A}}=({\mathbf{I}}_{np}-{\mathbf{H}})^{T}{\mathbf{W}}$ and ${\mathbf{B}}=({\mathbf{I}}_{np}-{\mathbf{H}})^{T}{\mathbf{R}}$, which are known $np\times np$ matrices. We also divide the matrices ${\mathbf{A}}$ and ${\mathbf{B}}$ into $n^{2}$ blocks of $p\times p$ matrices, and write ${\mathbf{A}}_{i,j}$ and ${\mathbf{B}}_{i,j}$, $i,j=1,\cdots n$, to mean the $i$th by $j$th blocks of ${\mathbf{A}}$ and ${\mathbf{B}}$ respectively. Hence ${\mathbf{A}}_{i,j}$ and ${\mathbf{B}}_{i,j}$ are both $p\times p$ matrices. In the supplementary materials we show that

[TABLE]

[TABLE]

We apply the $\mbox{vec}(\cdot)$ operator to both sides of the previous equation and use the identity $\mbox{vec}({\mathbf{A}}{\mathbf{X}}{\mathbf{B}})=({\mathbf{B}}^{T}\otimes{\mathbf{A}})\mbox{vec}({\mathbf{X}})$ (see Henderson and Searle, 1981), to obtain

[TABLE]

where

[TABLE]

[TABLE]

and

[TABLE]

Upon substituting $\mbox{E}[\mbox{btr}({\mathbf{Q}})]=\mbox{btr}({\mathbf{Q}})$, ${\bf\Sigma}_{\beta}={\hat{\bf\Sigma}}_{\beta}$ and ${\bf\Sigma}_{\omega}=\hat{{\bf\Sigma}}_{\omega}$ in equation (7), the method of moments gives one estimating equation in the vectorised form of two unknown covariance matrices.

4.3.2 Evaluating $\mbox{E}[\mbox{btr}({\mathbf{Q}}_{d})]$ and deriving the second estimating equation

Model (5) depends upon one unknown covariance matrix, ${\bf\Sigma}_{\beta}$. The intuition is that, upon using all $D$ of the ${\mathbf{Q}}_{d}$ matrices in (6) and the method of moments to estimate ${\bf\Sigma}_{\beta}$, we will then be able to estimate the other unknown covariance matrix ${\bf\Sigma}_{\omega}$ using the first estimating equation. We define design specific ${\mathbf{A}}_{d}=({\mathbf{I}}_{pn_{d}}-{\mathbf{H}}_{d})^{T}{\mathbf{W}}_{d}$ and ${\mathbf{B}}_{d}=({\mathbf{I}}_{pn_{d}}-{\mathbf{H}}_{d})^{T}{\mathbf{R}}_{d}$ , where ${\mathbf{A}}_{d}$ and ${\mathbf{B}}_{d}$ are known $pn_{d}\times pn_{d}$ matrices. We also divide the matrices ${\mathbf{A}}$ and ${\mathbf{B}}$ into $n_{d}^{2}$ blocks of $p\times p$ matrices, and write ${\mathbf{A}}_{d,i,j}$ and ${\mathbf{B}}_{d,i,j}$, $i,j=1,\cdots,n_{d}$, to mean the $i$th by $j$th blocks of ${\mathbf{A}}_{d}$ and ${\mathbf{B}}_{d}$ respectively. In the supplementary materials we show that

[TABLE]

where

[TABLE]

and

[TABLE]

Upon substituting $\mbox{E}[\sum\limits_{d=1}^{D}\mbox{btr}({\mathbf{Q}}_{d})]=\sum\limits_{d=1}^{D}\mbox{btr}({\mathbf{Q}}_{d})$ and ${\bf\Sigma}_{\beta}=\hat{{\bf\Sigma}}_{\beta}$ in (8), we obtain a second estimating equation from the method of moments.

4.4 Solving the estimating equations and performing inference

We solve the estimating equation resulting from (8) for $\mbox{vec}({\bf\hat{\Sigma}}_{\beta})$ and substitute this estimate into the estimating equation resulting from (7) and solve for $\mbox{vec}({\bf\hat{\Sigma}}_{\omega})$.

4.4.1 Estimating ${\bf\Sigma}_{\beta}$ under the consistent model

Some applied analysts may prefer to assume the consistent model (${\bf\Sigma}_{\omega}={\bf 0})$. As in the univariate case (Jackson et al., 2016), we have two possible ways of estimating ${\bf\Sigma}_{\beta}$ under the consistent model: we can use the estimating equation resulting from (7) with ${\bf\Sigma}_{\omega}={\bf 0}$ or the estimating equation resulting from (8) as in the full model. Also as in the univariate case, we suggest the former option because it uses the information made by assuming consistency when estimating ${\bf\Sigma}_{\beta}$. However this first option is valid only under the consistent model.

4.4.2 ‘Truncating’ the estimates of the unknown covariance matrices so that they are symmetric and positive semi-definite

As in the univariate case, there is the problem that the point estimates of the two unknown covariance matrices are not necessarily positive semi-definite. The method of moments does not even initially enforce the constraint that the point estimates of the unknown covariance matrices are symmetrical (Chen et al., 2012; Jackson et al., 2013). We produce symmetric estimators corresponding to an estimated covariance matrix of $\hat{{\bf\Sigma}}$ as $(\hat{{\bf\Sigma}^{T}}+{\bf\hat{\Sigma}})/2$ (Chen et al., 2012; Jackson et al., 2013). This also corresponds to taking the average of estimates that result from our ${\bf Q}$ and ${\bf Q}_{d}$ matrices and their transposes (Jackson et al., 2013). We then write these symmetric estimators in terms of their spectral decomposition (Chen et al., 2012; Jackson et al., 2013) and truncate any negative eigenvalues to zero to provide the final symmetric positive semi-definite estimated covariance matrices. Specifically, we define the truncated estimate corresponding to the symmetrical $\hat{{\mathbf{\Sigma}}}$ as $\hat{{\mathbf{\Sigma}}}^{+}=\sum\limits_{i=1}^{p}\mbox{max}(0,\lambda_{i}){\mathbf{e}}_{i}{\mathbf{e}}_{i}^{T},$ where $\lambda_{i}$ is the $i$th eigenvalue of the symmetric $\hat{\mathbf{\Sigma}}$ and ${\mathbf{e}}_{i}$ is the corresponding normalised eigenvector.

4.4.3 Inference for ${\bm{\delta}}$

Inference for ${\bm{\delta}}$ then proceeds as a weighted regression where all weights are treated as fixed and known. Writing $\hat{\bf{V}}$ as the estimated variance of ${\bf Y}$ in (3), in the absence of missing outcome data we have $\hat{{\bm{\delta}}}=({\bm{X}}^{T}\hat{\bf{V}}^{-1}{\bm{X}})^{-1}{\bm{X}}^{T}\hat{\bf{V}}^{-1}{\bm{Y}}$ where $\mbox{Var}(\hat{{\bm{\delta}}})=({\bm{X}}^{T}\hat{\bf{V}}^{-1}{\bm{X}})^{-1}.$ In the presence of missing data we can, under our missing completely at random assumption, apply these standard formulae for weighted regression to the observed outcomes. Alternatively and equivalently, we can impute the missing outcome data in ${\bf Y}$ with an arbitrary value and replace $\hat{\bf V}^{-1}$ with the precision matrix corresponding to $\hat{\bf V}$, calculated in the way explained for ${\bm{S}}$ in section 4.1 (Jackson et al., 2011). Approximate confidence intervals and hypothesis tests for all basic parameters for all outcomes then immediately follow by taking $\hat{{\bm{\delta}}}$ to be approximately normally distributed. Inferences for functional parameters follow by taking appropriate linear combinations of $\hat{{\bm{\delta}}}$.

4.5 Special cases of the estimation procedure

In the supplementary materials we show that the proposed method reduces to two previous methods in special cases. If all studies are two arm studies and consistency is assumed then the proposed method reduces to the matrix based method for multivariate meta-regression (Jackson et al., 2013). The proposed multivariate method reduces to the univariate DerSimonian and Laird method for network meta-analysis (Jackson et al., 2016) when $p=1$.

4.6 Model identification

If the necessary standard matrix inversions resulting from the estimating equations from (7) and (8) cannot be performed then both unknown variance components cannot be identified using the proposed method. A minimum requirement for any multivariate modelling is that the common-effect and consistent model must be identifiable. This means that there must be some information (direct or indirect) about each basic parameter for all outcomes. Two or more studies of the same design must provide data for all possible pairs of outcomes to identify ${\bf\Sigma}_{\beta}$. Two or more studies of different designs must provide data for all possible pairs of outcomes to identify ${\bf\Sigma}_{\omega}$. If these conditions are satisfied then the model will be identifiable. In situations where our model is not identifiable we suggest that simpler models should be considered instead. Possible strategies for this include considering models of lower dimension or the consistent model. In practice it is highly desirable to have more than the minimum amount of replication required, both within and between designs, so that the model is well identified. We make some pragmatic decisions in the next section for our example to provide sufficient replication within designs, in order to estimate ${\bf\Sigma}_{\beta}$ with reasonable precision.

5 Example

The methodology developed in this paper is now applied to an illustrative example in relapsing remitting multiple sclerosis (RRMS). Multiple sclerosis (MS) is an inflammatory disease of the brain and spinal cord and RRMS is a common type of MS. The effectiveness of a new treatment is typically measured to assess its impact on relapse rate and odds of disease progression. Magnetic Resonance Imaging (MRI) allows measurement of the number of new or enlarging lesions in the brain. Three outcomes are included in our analyses, so that $p=3$ in the full three dimensional network meta-analysis. These three outcomes are: (1) the log rate ratio of new or enlarging MRI lesions; (2) the log annualised relapse rate ratio; and (3) log disability progression odds ratio. Relapse is defined as appearance of new, worsening or recurrence of neurological symptoms that can be attributable to MS, accompanied by an increase of a score on the Expanded Disability Status Scale (EDSS) and also functional-systems score(s), lasting at least 24 hours, preceded by neurologic stability for at least 30 days. Disability progression is defined as an increase in EDSS score that was sustained for 12 weeks, with an absence of relapse at the time of assessment. Negative basic parameters indicate that treatments B-F are beneficial compared to treatment A throughout.

Data in this illustrative example were obtained from ten randomised controlled trials of six treatment options (coded in the network data as treatments A to F); placebo (A), interferon beta-1b (B), interferon beta-1a (C), glatiramer (D), and two doses of fingolimod; 0.5mg (E) and 1.25mg (F). Three of the fingolimod trials were three-arm (two doses and a control) and are included as three-arm studies. Three trials of interferon beta (one 1a and two 1b) were three-arm (also two doses and a control), and these were included as separate two-arm trials (each dose against the control, with the number of participants in each control arm halved). This ignores the differences in doses of interferon beta and was a pragmatic decision to help provide an identifiable network. Briefly, in this example there is very little replication within designs, so that identifying ${\bf\Sigma}_{\beta}$ well is very difficult without making pragmatic decisions such as this. Sormani et al. (2010) also treat these particular studies as two separate studies in this way, which helps them to identify their meta-regression models. Treating these three studies as separate two-arm trials means that the data are analysed as being from thirteen studies and a summary of the resulting data structure is shown in Table 1. There are eight different designs in Table 1 and so there is relatively little replication within designs, even when including three of the three-arm studies as separate two arm studies. See Bujkiewicz et al. (2016) for further details of these data. Figure 1 provides network diagrams that show the number of comparisons between each pair of treatments on the edges. In these diagrams the three arm studies (Table 1) are taken to contribute three comparisons, for example the CEF study contributes CE, CF and EF comparisons. Two estimates of treatment effect from this study contribute to analyses however because C is taken as the baseline; the study’s estimated EF treatment effect contains no additional information once its CE and CF contrasts are included in the analysis.

[FIGURE:]

Table 2 shows the estimates of the basic parameters (treatment effects relative to the reference treatment, placebo) obtained from univariate network meta-analyses, bivariate analyses for all three combinations of pairs of outcomes and the trivariate analysis. The results are similar across all analyses, and conclusions from univariate and multivariate analyses are the same. This is disappointing because multivariate analyses have not resulted in more precise inference. The entries of ${\bf\hat{\Sigma}}_{\beta}$ and ${\bf\hat{\Sigma}}_{\omega}$ are shown in Table 3. The positive estimates obtained for the unknown variance components suggest that this example exhibits some between-study heterogeneity and inconsistency. In order to assess the impact of the unknown variance components, we also fitted the consistent model and the common-effect and consistent model (results not shown) using all three outcomes ($p=3$). On average, the standard errors of the fifteen basic parameters from the full model are 35% greater (range: 13% to 84%) than those from the consistent model, which in turn are 58% (range: 8% to 128%) greater than those from the common-effect and consistent model. Both the between-study heterogeneity and inconsistency have notable impact.

The multivariate analysis adds to the univariate analyses in two main ways. Firstly, the finding that the multivariate analysis is in good agreement with the univariate analyses is a particularly important finding for treatment effects on MRI where a substantial proportion of data were missing. It has been demonstrated by Kirkham et al. (2012) that a multivariate approach to meta-analysis can help obtain more accurate estimates in the presence of outcome reporting bias. Hence the multivariate analysis reduces concerns that this univariate analysis is affected by reporting bias. Secondly joint inferences for all three outcomes are possible under the multivariate model. For example, and as we might anticipate, in our example the estimated log annualised relapse rate ratios and log disability progression odds ratios are highly positively correlated; from $\mbox{Var}(\hat{{\bm{\delta}}})$ in our three dimensional multivariate meta-analysis, the correlations between the five pairs of estimated basic parameters for these two outcomes are all between 0.63 and 0.75. Medical decision making based jointly on these two outcomes should take this high positive correlation into account, and this is only possible by using a multivariate approach. For example, a formal decision analysis involving these two outcomes should be based on their joint distribution rather than their two marginal distributions.

6 Discussion

We have proposed a new model for dealing with both multiple treatment contrasts and multiple outcomes, to provide a framework for conducting multivariate network meta-analysis. By using a matrix-based method of moments estimator, our methodology naturally builds on previous work (such as the well-known DerSimonian and Laird approach) and is computationally very fast, relative to other potential estimation approaches such as REML or MCMC; this is especially the case in very high dimensions and so our methodology is particularly advantageous for ambitious analyses of this type. The main disadvantage is that, as a necessary consequence of its semi-parametric nature, the method of moments is not based on sufficient statistics and so is not fully efficient. The loss in efficiency relative to maximum likelihood estimation awaits investigation but we anticipate that this will be less serious for inferences about the average effects than the unknown variance components. Furthermore the within-study normal approximations used in our model are not necessarily very accurate even in moderately sized studies.

Since our analysis uses a general design matrix, the modelling may easily be extended by adding study level covariates to describe and fit multivariate network meta-regressions. In the network meta-analysis setting these regressions have the potential to explain the reasons for inconsistency and model multiple dose level responses. Our method of moments estimation can be combined with approaches that ‘inflate’ confidence intervals from a frequentist random effects meta-analysis (Hartung and Knapp, 2001; Jackson and Riley, 2014).

In conclusion, we have developed a new model and estimation method for multivariate network meta-analysis, which can describe multiple treatments and multiple correlated outcomes.

Acknowledgements

DJ, IRW and ML are (or were) employed by the UK Medical Research Council [Unit Programme number U105260558]. SB was supported by the Medical Research Council (MRC) Methodology Research Programme [New Investigator Research Grant MR/L009854/1].